CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This applications claims priority to U.S. Provisional Application Serial No. 61/810,516, filed April 10, 2013, titled "AUTOMATED MEASUREMENT OF HETEROGENEITY IN CT IMAGES OF HEALTHY AND DISEASED LUNGS USING VARIOGRAM ANALYSIS OF AN OCTREE DECOMPOSITION," and U.S. Serial No. 14/166,601 , filed January 28, 2014, titled METHOD OF ASSESSING HETEROGENEITY IN IMAGES".

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH OR DEVELOPMENT

|0002| The invention was made with Government support under Contract DE-AC05- 76RLOI 830, awarded by the U.S. Department of Energy. The Government has certain rights in the invention.

TECHNICAL FIELD

[0003] This invention relates to three-dimensional image analysis. More specifically, this invention relates to assessing heterogeneity in three-dimensional images using comparative analysis of decomposed images.

BACKGROUND OF THE INVENTION

[0004] Early detection of disease is generally associated with improved health outcomes. For example, in emphysema the Global initiative for Chronic Obstructive Lung Disease (GOLD) staging is a well-accepted approach for categorizing disease as mild, moderate, severe, or very severe (Pauwels et al. 2001 ). However, GOLD scores are global with no regional information, in addition, they are based only on spirometry and clinical

presentation. Thus, GOLD is a helpful non-invasive starting point, but uncertainty about its diagnosis may warrant additional diagnostic procedures such as CT (computed tomography) imaging (Pauwels et al. 2001 ). Furthermore many diseases, such as emphysema, are heterogeneous diseases, and signatures of tissue heterogeneity may facilitate more accurate diagnoses. An image-based metric would be useful for detecting such heterogeneities.

SUMMARY OF THE INVENTION

[0005] The present invention is directed to methods of assessing heterogeneity in images. In one embodiment, a method of assessing heterogeneity in images is disclosed. The method includes acquiring 3D images of an object. The method also includes performing decomposition on the acquired images. The method further includes performing comparative analysis of the decomposed images using distance.

[0006] The method may further include masking and filtering the acquired images.

[00071 In one embodiment, performing decomposition on the acquired images includes performing iterative decomposition on the acquired images. Performing iterative decomposition on the masked images ma include performing octree decomposition on the masked images, and performing octree decomposition on the masked images may include reducing the images to obtain image subdivisions that are relatively homogeneous. Reducing the images to obtain image subdivisions that are relatively homogeneous may include dividing boxes of the images according to a criterion. The criterion for dividing boxes of the images may include a measurement of standard deviation within each box and comparison of the standard deviation to a threshold standard. Each box may include location and mean value information. [0008j I ^R one embodiment, performing comparative analysis of the decomposed images using distance includes performing variogram analysis of the decomposed images.

Performing variogram analysis may include calculating a square of the difference of the means between each box. The box may be, but is not limited to, an 8x8x8 voxel box.

Calculating a square of the difference of the means between each box may include determining a spatial relationship between regions of the images that are relatively homogeneous.

[0009] The 3D image may be, but is not limited to, biological tissue.

[0010] In one embodiment, the images are reduced to homogeneous blocks of 2x2x2, 4x4x4, and 8x8x8 voxels.

[0011] In one embodiment, the threshold standard includes a Hounsfie!d Unit (HU) threshold.

[0012] The image are, but not limited to, at least one of the following: CT images, multidimensional images, ultrasound images, MRI images, data arrays, data matrix, data sets, pictures, multi-spectral images, and voxels of data generated using histograms.

[0013] In one embodiment, performing comparative analysis of the decomposed images using distance includes performing correlogram analysis of the decomposed images.

[0014] In another embodiment of the present invention, a method of assessing heterogeneity in images is disclosed. The method includes acquiring 3D images of an object; filtering the acquired images; masking the images; performing iterative decomposition on the masked images to obtain image subdivisions that are relatively homogeneous; and performing variogram anaiysis or correlogram analysis of the octree decomposed images to determine spatial relationships between the regions of the images that are relatively homogeneous. BRIEF DESCRIPTION OF THE DRAWINGS

[0015 J Figure 1 shows an example octree decomposition for reducing images to obtain image subdivisions.

[0016] Figure 2 shows unfiltered coronal slices from 3D images of (A) a control rat, (B) a rat dosed to the full lung with elastase to induce mild emphysema, and (C) a rat dosed to the distal portion of the left lobe with elastase to induce local emphysema.

[0017] Figure 3 shows an example of an image that was masked and is also (A) an unfiltered image, (B) the image with a five voxel diameter 3D median filter, and (C) the same image downsampled to 1/8 resolution, or 64 x 64 x 64 pixels.

10018} Figure 4 shows octree decomposition results from two rats. The top row is for a control rat and the bottom row is a single-lobe-dose rat (the same as those shown in Figure 2). Columns A, B, C, and D show the resulting 2 x 2 x2, 4 x 4 x 4, 8 x 8 x 8, and 16 x 16 x 16 boxes, respectively.

[0019] Figure 5 shows graphs of octree decomposition histograms showing the relationship between box size. HU. and percentage of the lung by each, box size at each HU level for (A) control, (B) full-lung dose, and (C) single-lobe dose.

100201 Figure 6 shows variogram s averaged for each dose group. The dashed line shows the extent of d _max , the characteristic lung distance to which variograms were analyzed for heterogeneous variations.

[0021 ] Figure 7 shows bar graphs of (A) percentage of image pixels w ith HU values below two standard deviations from the mean of the control group, (B) coefficient of variation (CoV), and (C) the average distance of each individual rat's variogram from the mean of the control group. The x-axis labels are animal numbers: 100s are controls, 200s are full-lung-dose, and 300s are single-lobe-dose. 301 and 302 received the dose in the right caudal lobe, and 303 received the dose in the distal region of the left lobe (Figure 2C). [0022] Figure 8 is a table that summarizes by measurement type and dose group the pulmonar function tests (PFT) results and other physiologic measurements.

[00231 Figure 9 shows variograms from 0- and 15-Gy rais showing the average HU variance of 8 x 8 x 8 cubes versus distance of separation (A) from computed tomography images at functional residual capacity and (B) volume maps.

[0024] Figure 10 shows reconstructed variograms based on the mean parameters of each dose group for the CT images at (A) functional residual capacity and (B) volume maps.

[0025] Figure 1 1 shows box plots of variogram parameters versus dose of (A) the exponent a from the computed tomography images at functional residual capacity and (B) range (R) and still (S) from volume maps.

[0026] Figure 12 shows variograms constructed from all the image sets, out to a distance of 100 mm.

1 027] Figure 13 is a box plot showing the group medians of the exponent a.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

(0028] The present invention is directed to methods of assessing heterogeneity in images. The present invention combines methods of 3D image analysis and applies them to the problem of discerning disease signatures in organs, also known as image-based biomarkers. Octree decomposition, which iteratively subdivides an image - with each division producing eight evenly sized octants (e.g., an initial 512 x 512 x 512 region decomposed into eight 256 x 256 x 256 octant regions, and so on) - reduces the image to boxes of maximum size with minimal heterogeneity w ithin each box. Octree decomposition isolates regions of tissue that are relatively homogeneous and calculates their spatial relationship.

[0029] The present invention also includes comparative analysis, such as variogram or correlogram analysis, in which the spatial relationships between the boxes are examined in, for example, a variogram - a common geostatistical tool that compares variance to distance. [0030] The present invention also discloses a heterogeneity metric that relates the relationship of the variograms of diseased tissue to healthy tissue.

[0031] The present invention also eliminates portions of the tissue that do not contribute to the disease. For example, in the lung, the present invention generally eliminates boundaries, airways, vasculature, etc. from consideration and focuses on the diseased tissue.

[0032] The present invention also provides a single metric that can be used to distinguish disease from healthy tissue.

Experimental Section

[0033] The following examples serve to illustrate embodiments and aspects of the present invention and are not to be construed as limiting the scope thereof.

EXAMPLE 1

Materials and Methods

Disease Model

[0034] An elastase-induced model of emphysema was used in this study. Nine male Sprague-Dawley rats with an average weight of 212 ± 1 1 g were orally intubated and dosed intratracheal ly with: 250 U/kg elastase dissolved in 200 μΙ_- saline to the whole lung (n=3), or 50 U kg elastase in 200 μί saline to a single lobe (n=3), or 200 ih saline as a control (n=3). Dosing levels were based on our previous work in which emphysematous changes were detected using "He diffusion MRI and histology (Jacob et al. 2008). All animal use was approved by the Institutional Animal Care and Use Committee at Paci ic Northwest National Laboratory (PNNL). CT IMAGING

[0035] Three weeks following dosing, the rats were imaged using micro-CT. At this time, rats weighed 357 ± 10 g. This imaging procedure is described in detail in Jacob et al. 201 1 . Briefly, rats were anesthetized, intubated, and mechanically ventilated at 1 Hz with 40% inhale and 60% exhale durations. Peak inhalation pressure was ~8 cmH ₂ 0, and no peak end expiratory pressure was used so that images could be acquired at functional residual capacity. Anesthesia was maintained by providing 3-4% isotlurane in air (30% 0 ₂ , balance N ₂ ). Sigh breaths were delivered periodically to maintain lung recruitment. A respiratory- gated GE eXp!ore 120 micro-CT scanner was employed with the following settings: 100 kV peak voltage, 50 mA tube current, 16 ms exposure time, and 360 projections with 1 degree angular steps. Images were reconstructed with supplied software to 150 μιτι isotropic resolution. Total imaging time was about 90 minutes due to the collection of multiple images throughout the breathing cycle; however, only the images acquired at full exhalation were used for the analyses herein. Post-mortem histology and tissue analysis were not possible due to in situ lung casting for other research purposes (i.e. supplying airway tree geometries for computational fluid dynamics models (Corley et al. 2012, Minard et al. 2012)).

Image Preparation

[0036] A lung mask image was semi-automatically generated from the above-mentioned reconstructed images using ImageJ and the ImageJ 3D Toolkit plug-in. Starting from a seed- point inside the lung, a 3D connected threshold was applied with the threshold value empirically selected to exclude major vasculature and all external tissue, so that only lung tissue would be included in analysis. Mask boundaries were smoothed using the Region Dilate and Region Erode functions in succession. Applying the mask to the reconstructed images, all background was thus assigned an intensity value of 0, and all unmasked lung tissue retained original HU values. A five pixel diameter 3D median filter was then applied in order to reduce noise while contributing minimal blurring. Finally, the image canvas size was increased to 512 x 512 x 512 by zero-filling in each direction. Zero-filling to a power of

2 served to conveniently restrict all octants in the octree decomposition to isotropic cubes. Octree Decomposition

[0037] An automated octree decomposition was performed by iteratively subdividing an image, with each division producing eight evenly sized octants (e.g., the initial 512 x 512 x 512 image was decomposed into eight 256 x 256 x 256 octant regions, and so on, as shown in Figure 1. After each division, the maximum, minimum, mean, and standard deviation of the signal intensity were calculated for each octant region. Then, iterative octant subdivision of a region occurred if either the standard deviation of that region exceeded a user-defined threshold, or if the region contained at least one voxel - but not all voxels - with intensity equal to the background signal (i.e. 0). However, regions reaching a minimum size of 2 x 2 x 2 were not subdivided. At the end of the decomposition, all boxes with a mean HU value≠ 0 and dimensions greater than 2 2 x 2 represented portions of lung with relatively uniform HU values. Following decomposition, results were binned according to HU and box size for histogram display. Boxes containing at least one 0 were not included in the binning. Octree decomposition was implemented in Python and executed on a Mac Pro model 3.1 .

[0038] There are multiple approaches to selecting the aforementioned user-defined threshold criteria. For example, boxes may be subdivided if the HU range in the box exceeds 10% of the range of values in the entire image (Subramaniam et al. 2007). Because the rats in this study had different ranges of HU values depending on disease severity standard deviation was used as the threshold criterion, as it does not require normalization.

V'ariogram Analysis

[00391 Vasograms were used in this context to determine the spatial relationships between regions of the lung that are relatively homogeneous. The octree decomposition isolated regions of relatively high homogeneity; therefore, each resulting cube was treated as a homogeneous "voxel" using the mean value of the cube as the voxel intensity. Furthermore, octree cubes of the same size were used. The largest cube sizes that resulted from the octree decompositions were mostly 8 8 x 8 with a few 16 x 16 x 16. In order to ensure a sufficient number of cubes that were distributed throughout the lung and to enhance statistical power (Keil et al. 2012), each 16 x 16 x 16 cube was further broken down into 8 x 8 8 boxes, and all 8 x 8 x 8 boxes were used in the variogram analysis. This approach assumed that emphysematous disease was manifest on a scale greater than an 8 x 8 x 8 cube, which for these images is 1.2 mm on a side - this is consistent with what observed in this disease model in previous work (Jacob et al. 2008). Excluding the smaller boxes largely eliminated vasculature, conducting airways, boundaries, and other features that had high spatial variability (Zhang et al. 2005).

[0040] The calculation of the variance (or semi-variance) y(d) f the differences in signal intensity / is described i detail in many geostatistics texts (Dua et al. 2010, Gringarten et al. 2001 ) and is given by:

[00411 where d is the distance between voxel centroids. x is the voxel location in 3D, and N is the number of voxel pairs. For an image with n voxels,

N = n(n - l)/2. [2]

[0042] The average semi-variance at each distance was calculated. Results were plotted on a distance vs. variance graph, or variogram, which graphically represents the spatial dissimilarity within the image. {0043] Variograms characterize spatial relationships between boxes of relatively high homogeneity that were dispersed throughout the lung. However, as described in previous work with brain images ( eil et al. 2012), increased distances between voxels, or octants, lead to decreased likelihood that the voxels are related in any way. This is especially true in the lung, which is composed of largely independent lobes (five in the rat). Moreover, the spatial relationship between parenchymal signal intensity in different lobes is most accurately described by a unique path up and back down the airway tree, potentially covering 20 or more orders of branching (Schulz 2000). The spatial distances are analogous to the

"regionalized variable ^" in the geological context of mineral distribution. Indeed, the regionalized variable, or the maximum distance d _max over which the variogram is expected to be reliable, is estimated to be one half of the diameter of a region (Keil et al. 2012).

Therefore, the variogram analysis in this study was limited to a region that represented half of the average dimension of the left lobe, which is ty pically the largest lobe in the rat. This assumes that the disease varies slowly compared to the box size but rapidly compared to the lobe size. Using the 3D images and a lung cast, d _max was measured to be ~ 8 mm (53 pixels) in the healthy rats, w hich was rounded up to the equivalent of seven face-bordering 8 8 8 boxes, or 56 pixels.

Octree Decomposition Tests

{0044] The following evaluations of the appropriateness and robustness of the octree and variogram approach were performed.

Threshold Range

{0045] The octree decomposition compared the standard dev iation of the lung signal within a box to a threshold range to determine whether the box should be subdivided. The optimal threshold range was determined semi-empirically. As a starting point, the typical standard deviation of the lung tissue of the three control rats was calculated. This was accomplished by producing a histogram of each masked image and then fitting the main lung peak to a Gaussian curve. The mean of the control rats ^" standard deviations (σ) was set as the initial threshold range. The effects of thresholds that were σ/3, 2σ/3, 4σ/3, and 2σ were tested by rerunning the decomposition on the same images. The effectiveness of each threshold range at grouping the control rats and at distinguishing the full-lung-dose rats from the control group in the variograms was evaluated.

Image Filtering

[0046] The following were examined: the effects of image filtering by generating variograms from un filtered images, images filtered with a five pixel diameter 3D median filter, images filtered with a 3D Gaussian blur of radius=2, and images filtered with a 3D Gaussian blur of radius=4. Filters were applied in I mage J prior to applying the mask. The number of octree boxes and differences in the variograms were compared. Adjustments to the threshold range were made for each filter test based on the standard deviation of the filtered image.

Image Translation

[0047] Rat positioning during imaging can vary from animal to animal, and the octree decomposition should be independent of this. Therefore, the effects of image translation on the variogram results were tested. For the variogram analysis, the largest unit box retained after decomposition was 8 8 x 8; therefore, a shift in the image by 8 pixels in any (or all) dimensions should result in no change to the resulting variogram. since the original image had isotropic dimensions of 2". Conversely, maximum change would occur from a four-pixel shift. To test the significance of the effects of translation, the image of one rat was shifted in X. y, and z by four pixels and by eight pixels using ImageJ, and variogram results were compared to those of the original image.

Image Rotation

[0048] Animal positioning can also have an apparent effect on image rotation, and octree decomposition should be insensitive to this. Using the same rat image as used in the translation test, an arbitrary 3D rotation was imitated by rotating the image π/4 radians about the x, y, and z axes using ImageJ. Results were compared to the original image and to the results of the translat ion test.

Image Downsampling

[0049] To demonstrate the increased value of the combined octree and threshold approach for establishing the 8 x 8 x 8 boxes versus simply downsampling the image, variograms were generated using standard image downsampling as a comparison. For this, a bilinear downsampling to the 512 5 12 x 5 12 image was applied, creating a 64 x 64 x 64 image - each voxel the equivalent size of an 8 x 8 8 octree box. Then variograms were generated utilizing the intensity value of every voxel in the 64 x 64 x 64 image (excluding the background voxels, again defined as those voxels with intensities equal to zero).

Emphysema Index and CoV

[0050] The percentage of lung below a HU threshold value, or emphysema index, was calculated for each rat to compare this conventional measurement of disease severity with the variogram results. Because there is no established HU threshold level in rat models of emphysema, the percentage of voxels with HU values below two standard deviations from the control-group mean were counted. This level was determined to be -717 HU. Calculations were made on the same masked images used for octree decomposition.

[0051] The CoV was calculated for each rat by fitting a histogram of the masked lung images to a Gaussian curve and taking the ratio of the standard deviation to the mean.

Statistical Analysis

[0052] Comparisons of variograms of the three control group rats were made using the Kruskal-Wallis (KW) rank sum test with a null hypothesis of a=0.05. This indicated whether there were significant differences within the group. For pairwise comparisons, a Mann- Whitney- Wilcoxon (MWW) rank sum test was employed, also with a null hypothesis of a=0.05.

Results

CT Images

[0053] Figure 2 shows a representative unfiltered corona! slice from a rat in each dose group: panel A is a control rat, panel B is a full-lung-dose rat, and panel C is a partial-lung- dose rat. The dose was delivered to the distal portion of the left lobe. Overall, it was difficult to discern qualitative differences between rats in the control group and the full-lung-dose group: the mean HU value (± standard deviation) of the former was -544 ± 58 and the mean of the latter was -594 ± 38. in spite of marginally lower HU values that would be expected from emphysematous disease, an analysis of variance showed that the two groups were indeed not statistically different (p=0.26). On the other hand, the partial-lung-dose rats showed a bimodal distribution o HU values because the diseased regions of the lung were distinct and considerably different than the healthy regions. An example of this can be seen in the lower portion of the left lobe of the partial-lung-dosed rat (panel C) where the signal intensity is substantially lower than the rest of the lung by -200 HU - indicative of severe tissue destruction and air trapping, which are characteristics of emphysematous lungs (Spencer 1985). Based on the overall HU measurements, it is presumed that the full-lung- dose rats developed a mild - and difficult to distinguish - emphysematous disease while the single-iobe-dose rats developed a more severe, albeit localized, disease.

Octree Decomposition Tests

Threshold

[0054] The starting octree decomposition threshold level determined from the standard deviations of the control rats (Section 2.6.1 ) was 60 ± 1 HU: therefore, 60 HU was used as the initial threshold level. Thus, the other threshold levels tested were 20, 40, 80, and 120 HU. Octree decompositions were performed using the different threshold levels on the three control rats, and variograms were compared using a KW rank sum test. The threshold level that showed the least difference among the control rats was 40 HU (ρ=Ό.12), with 60 HU also showing no significant differences (p=0.09). The 20, 80, and 120 HU thresholds did have significant differences among the controls (p=0.02, pO.0001 , and pO.0001 , respectively). Based on this result, the threshold level of 40 11 U was used for all subsequent octree decompositions (unless otherwise noted).

Image Filtering

[0055] Results of applying the different filters showed that the variograms from the un filtered image and the image with the median filter were indistinct in an MWW test (p=0.62). However, the relatively noisy unfiltered image resulted in about 20% fewer 8 x 8 x 8 boxes than the filtered image. This was in spite of a higher threshold that was used, 48 HU, based on the unfiltered image ^' s standard deviation. The radius=2 and radius=4 Gaussian filter results also did not differ significantly from the unfiltered image (p=0.41 and p=0.48, respectively) or from the median-filtered image (p=0.71 and p=0.84, respectively). A KW test showed that the radius=2 filter resulted in control group variograms that were not distinct (p=0.16), but the radius=4 variograms were (p=0.QQGl). However, because of the blurring caused by Gaussian filters, the number of 8 8 x 8 boxes that resulted from the Gaussian- filtered images was about 30% higher than the median filtered image, presumably because vascular structures and edges were not well preserved. Threshold levels used for the radius=2 and radius=4 images were 37 1-1 LI and 39 HU, respectively, based on the post-filtering standard deviations.

Effect of Downsampling

[0056] Figure 3 shows an example of an original image (panel A), the image with the five voxel diameter 3D median filter (panel B), and the same image downsampled to 1 /8 resolution, or 64 x 64 x 64 pixels (panel C). The variograms made from the downsampled images of the three control rats were statistically different in a KW rank sum test (pO.OOO l ). On the other hand, the variograms made from only the 8 8 8 boxes that came out of the octree decomposition were not (p=0.12).

Translation and Rotation

[0057] Results of image translation and rotation were compared for one rat. For the eight-pixel translation, variogram results were exactly identical to thai; of the original image, as expected. The four-pixel shift was not statistically distinct from the original image as confirmed by a M WW test (p=0.96). In addition, the result of the rotation showed no significant change from the original (p=0.56) or from the four-pixel shift (p=0.59). Thus, it is confirmed that shifts or rotations to the image (i.e. alternative positions of the lung during imaging) do not result in significant changes to the resulting variograms.

Octree Decomposition

[0058] Figure 4 shows the results of the octree decomposition on a control rat and on one with severe disease in the lower left lobe (see Figure 2C). Column A shows the 2 x 2 x 2 boxes, and column B shows the 4 x 4 x 4 boxes. These box sizes largely define the fine structures and edge details, including the conducting airways and vasculature. For this reason, we ignored the 2 x 2 x 2 and 4 4 x 4 boxes for the variogram analysis. Column C shows the 8 8 8 boxes and their relatively uniform distribution throughout the lung. Conversely, the less uniform distribution of the considerably fewer 16 x 16 x 16 boxes - the largest that resulted from octree decomposition - is shown in Column D. The lower left lobe of the treated rat (column D, bottom) has a noticeably higher density of 16 x 16 x 16 boxes than that of the control rat. indicating homogeneous and prevalent, albeit localized, CT signal intensity.

[0059] Histograms of the octree decomposition for three rats are shown in Figure 5. The marker size corresponds to the box size, in the control and mild disease rats (panels A and B), the different sized boxes, other than the 2 x 2 x 2's, are approximately Gaussian distributed about the mean HU value of the lung, whereas the single-!obe-dose rat (panel C) shows a bimodal distribution, indicative of at least one large region of the lung with considerably lower HU values. These histograms show what percentage of the lung at each HU intensity is defined by the different box sizes, but they do not convey any information about the spatial relationships of the boxes.

Variograms

[0060] In order to visualize the spatial relationships between the 8 x 8 x 8 boxes (and decomposed 16 x 16 x 16 boxes), variograms were constructed. Figure 6 shows the average variograms from the three different dose groups. The dashed line in Figure 6 denotes the range of d, _mx (see section 2.5). For d < d _max , A KW test showed that the control rats were statistically similar (p=0.12), and a MWW test showed that the dose groups were each statistically distinct from the control group (p<0-0GGl ) and from one another (pO.OOO i ). Emphysema Index, CoV, and Heterogeneity Score

[0061 ] Results of the emphysema index are shown in the Figure 7 A. Although the number of rats is too small to reliably calculate sensitivity and specificity, the graph shows that there is considerable overlap between subjects in each dose group, likely indicating poor sensitivity and specificity . The CoV for each rat is shown in Figure 7B. Similarly, there is no distinction between the control and full-lung dose groups. The part-lung dose group had considerably higher CoV values, because the standard deviation was enlarged due to a bimodal (and non-Gaussian) distribution of HU values. [0062] For comparison, the average distance Δ of each rat's variance from that of the control group average was calculated, as shown in Figure 7C. The lowest full-lung dose rat has a Δ that is 2x higher than the highest control rat. There is no overlap between dose groups, implying better sensitivity and specificity than the other metrics.

Discussion

[0063 J This work combined two different data analysis tools, octree decomposition and variograms, to study tissue heterogeneity in lung disease. This merged approach was better able to differentiate rats with mild emphysematous disease from the healthy control group than methods that relied on absolute HU values. The main criterion for octree decomposition was based on the standard deviation of HU values within an octree box. An advantage to this approach is that it avoids thresholding according to HU values; rather, it focuses on heterogeneity-based signatures that may characterize disease (Mets et al. 2012). It is proposed that a heterogeneity score Δ, the average distance of a rat's variance from that of the control group average, may be useful to classify disease severity. Furthermore, to visualize the regions of the lung with the greatest heterogeneity, one could determine which boxes had the highest semi-variance within d _max and map them back to the original image. This would provide 3D information about the spatial distribution of lung tissue heterogeneity and, potentially, disease distribution.

[0064] Another approach to a disease metric might be that of fitting the data to an established variogram model, most of which describe an asymptotic rise in variance (i.e. variance becomes independent of distance indicating that spatial relationships become random (Clark et al. 2000)). This is seen to some degree in the data of Figure 6. However, within the range of d<d _max , it was found that a power law model generally fit the data better than other models. This model typically describes fractal behavior (Bohling 2007). nitial investigations into this model showed potential promise at using fit parameters to distinguish dose groups; however, results lacked statistical significance, and the full-lung dose group tended to not fit this model as well as the other two groups.

[0065] There was no single variogram model that satisfactorily fit the all the data over the entire range of distances, because the complex geometries found in the lung result in some problems for variograms. In particular, the direct linear path between two regions of the lung separated by large distances often crosses non-lung tissue, such as the heart, which are essentially treated as holes or voids in the geometry. Furthermore, neighboring lobes generally do not interact physiologically except through the vascular and airway trees, which may only connect regions through many orders of branching. To limit the problem in this study, the distance o variogram analysis was constrained to approximately half the characteristic diameter of the largest lobe. However, to better understand the inter- and intra- lobe variance relationships, the lung could be segmented into lobes (if the image is of sufficient resolution to discern lobar boundaries), and the decomposiiion/variogram process repeated on the segmented images. Though, by ignoring inter-lobar variances, this approach would likely not capture information about disease that was confined to a single lobe, particularly if the entire lobe was affected homogeneously.

[0066] Octree decomposition was employed prior to generating the variograms. In doing this, it was assumed that the emphysematous disease is generally slowly varying over space, and that the disease causes changes to homogeneity on the order of or less than the lobar length scale but greater tha the octree box length scale (Litmanovich et al. 2009, Spencer 1985). This is consistent with what was previously observed using 3 He MRI in the same disease model (Jacob et al. 2008). Without this assumption, variograms would have to be made directly from the raw images. This is possible but impractical, because the semi- variance computation time (and resulting file size) is proportional to the number of voxel pairs, which rises approximately as n2 (see Eq. 2). The variogram computation time was measured for 1078 octree-decomposed boxes from one rat to be 6.6 seconds (on a MacPro model 3.1), and it was verified experimentally that the computation time indeed rose in proportion to the square of the number of voxels. Therefore, to create a variogram on the entire masked 3D image, which consisted of 1 .48x 106 voxels (of lung tissue only), it would take -1.25x 107 seconds, or about 5 months - with a resulting file size on the order of 20 GB. Therefore, octree decomposition dramatically reduces the computation time while focusing non-objectively on regions of the lung that are of greatest interest. The octree decomposition itself was performed in ~2 minutes.

[0067] An alternative to octree decomposition is downsampiing, which is a quick and straightforward approach to reducing image size. However, it was verified that the octree decomposition approach performed much better at separating dose groups than simply downsampiing the image and then calculating the semi- variance using every voxel. The octree decomposition assures that only the homogeneous regions of the lungs are singled out for comparison, whereas downsampiing the image blurs together proximal voxels, including vasculature, airways, and lung boundaries irrespective of signal intensity or tissue type. Thus, the downsampiing approach apparently causes a loss of information. The radius=4 Gaussian filter had a result similar to downsampiing.

EXAMPLE 2

[0068] This second example demonstrates the 3D imaging and analysis approach using variogram s of octree-decomposed images, as described above, to detect subtle injury in radiation-exposed rat lungs. It has been shown that lung injury in radiation-exposed rat lungs includes acute inflammation in the weeks after exposure and chronic fibrosis in the months after exposure while providing a relatively uniform injury over the dosed area (Lehnert et al. 1989, Ward et al. 2004, and Ward et al. 1993). It is shown that the results of the variogram analy sis on octree-decomposed CT images and 4DCT-based volume maps of irradiated lungs correlates with radiation dose better than physiologic measurements, conventional pulmonary functions tests, or CT density measurements.

MATERIALS AND METHODS

[0069] All animal use followed a protocol approved by the Institutional Animal Care and Use Committee at PNNL. Twenty-five female Sprague Daw ley rats weighing 216 ± 9 grams were used.

[0070] Lung irradiations employed a -6400 Ci Co-60 gamma source with a -30 cm thick lead collimator. The collimator was trapezoidal to match the basic outline of the lungs, w ith dimensions based on a 3D image of a weight-matched control rat: 15 mm wide at the top, 30 mm wide at the bottom, and 26 mm high. The bottom was convexly tapered by 5 mm in the center to minimize exposure to the liver. The resulting beam was then calibrated in terms of absorbed dose to tissue (Gy/min) using a tissue-equivalent ionization chamber (Exradin model A 12, Standard Imaging, Middleton, WI) connected to an electrometer (model 617, Keith!ey Instruments, Cleveland, OH) to collect the resulting current. Appropriate corrections were applied to convert from exposure in air to absorbed dose in tissue (Report 30. 1979). The measured absorbed dose rate for the estimated location of the lung center (-1.5 cm from collimator face) was 3.74 Gy/min. It was determined that 1 .5 cm of tissue results in approximately 5% reduction in absorbed dose rate, resulting in an estimated 3.5 Gy/min at lung center. One of five calculated doses of 0.0, 5.9, 8.8, 1 1.8, and 14.7 Gy (hereafter referred to as 0, 6, 9, 12, and 15 Gy) was delivered to the thorax; this dose range has been shown to cause significant injury in other rat strains (Ward et al. 1993, Ghosh et ai. 2009, and Novakova-Jiresova et al. 2007). The dose rate to the body was measured to be <0.3% of that at the lung center. The dosed region was confirmed in a weight-matched rat using an x-ray source and Polaroid radiographic film.

[0071] Anesthetized rats were placed in a custom-made, contoured holder to facilitate reproducible positioning of the rat thorax directly in front of the collimator. Rats were randomly assigned a radiation dose, with 5 rats per group. Irradiations were blind to the staff performing other measurements in order to reduce bias.

[0072] Following irradiation, rats were returned to the animal facility where they were individually housed, provided food and water ad iibitium, and observed daily for general well-being. One rat from the 6 Gy group developed a ~3 cm growth on its back and was eliminated from the study. Otherwise, no mortality or outward signs of poor health were observed. Two additional rats were kept in the same room as health sentinels, and at the end of the study they were confirmed to be seronegative for common rat pathogens.

[0073] At 6 months post-irradiation, rats were subjected to pulmonary function tests (PFT) and microCT imaging. First, a whole body plethysmography (WBP; Buxco Research Systems, Wilmington. NC) was used to measure breathing rate, tidal volume, and minute volume of unanesthetized, unrestrained rats for ~5 minutes. Rats were acclimated to the chamber for -~ 10 minutes per day for several days prior to the test.

[0074] Next, rats were imaged using 4DCT (multi-time-point 3D imaging). Details of animal preparation, ventilation, and 4DCT imaging closely follow those described in Jacob et al 201 1. Rats were anesthetized with 4% isoflurane in oxygen, orally intubated with a 14 gauge catheter tube, and connected to a computer-controlled mechanical ventilator (model 830/AP, CWE Inc., Ardmore, PA). Rats were maintained on isoflurane and ventilated with 30% 0 ₂ (balance 2) at 54 breaths per minute, with a 500 ms inhale duration and no breath hold. Periodic sighs were delivered every 100 breaths to maintain lung recruitment. The ventilator recorded tracheal pressure, inspiratory volume, and expiratory volume. Peak inspiratory volume (PIV) was -2.1 mL. No positive end expiratory pressure was used so that images could be acquired at full passive exhalation when the lung volume was at functional residual capacity (FRC). A microCT scanner (eXplore 120, GE Healthcare, Waukesha, WI) with ventilatory gating was used to acquire eleven images throughout the breathing cycle in 26 minutes with 100 ms temporal resolution (It should be noted that on Iy the images at the breathing cycle extremes - FRC and PIV - were analyzed for this study). Gating was tested by comparing the ventilator gate signal to a signal sent by the scanner upon firing x-rays; the delay between the ventilator trigger and x-ray firing was found to be < 250 με. CT imaging parameters were: 80 kVp, 32 mA, 16 ms exposure time, and 360 projections with 1° angular separation. The estimated radiation dose from all images was 940 mGy. Images were reconstructed to 200 μι isotropic resolution using supplied software, it was empirically found that this resolution provided sufficient detail for later analysis.

[0075] Following imaging, PFT were performed. Rats were anesthetized with an intraperitoneal injection of Ketam i ne/Xy laz i ne and surgically intubated for measurement of inspiratory volumes, forced expiratory volumes, and quasistatic chord compliance using a Forced Maneuvers system (Buxco Research Systems, Wilmington, NC). PFT measurements took ~5 minutes. Immediately after the pulmonary function tests, animals were euthanized via CC asphyxiation, and lungs were excised to obtain wet and dr weights.

10076} Acquired images were processed with a 5 pixel diameter 3D median filter to remove noise while maintaining feature boundaries. These filtered images were then masked to assign intensity value 0 to non-lung regions; filtered images were multiplied with a binary image created using the 3D connected threshold tool in the 3D Pl gins Toolkit of Image J to delineate lung from non-lung based on intensity . From these masked images, the mean and standard deviation σ of the distribution of Hounsfield units (HU) in the lung were determined. The coefficient of variation (CoV) was then calculated by taking the ratio of σ to the mean.

[0077] Octree subdivision reduces the original volume into eight octants of equal dimension (see Figure 1). Images were first zero-filled to 256x256x256 so that octree subdivision would result in isotropic cubes. Iteratively, each octant was further subdivided if it contained any mask voxels or if its σ exceeded a threshold value t, unless the octant reached a minimum size of 2x2x2. All octants containing the mask value of 0 were discarded from further analysis. Remaining octants larger than 2x2x2 represented relatively homogeneous sections of the lung, while the 2x2x2 cubes defined boundaries and regions with high spatial variability. In developing this approach, the optimal was empirically determined to be approximately two-thirds the mean of the control group ^' s whole-lung σ; t was 38 HU in this study. Octree decomposition was executed in Python on a Mac Pro 3.1 in ~2 minutes.

[0078] A ter octree decomposition, variogram analysis was performed to determine spatial relationships between octants. As mentioned above, a variogram is a distance vs. variance (or semi-variance) graph, whereby variance is the square of the difference in mean signal intensity between two given octants, and distance is between octant centroids. Variograms were constructed using only 8x8x8 octants - the largest that were generated from the decomposition. This selection assumes that any radiation-induced abnormalities would generally be manifest on a scale larger than an 8x8x8 cube, which had dimensions of 1.6 mm x 1.6 mm x 1.6 mm. Furthermore, this eliminated the smaller cubes that tended to define features of high spatial variability. Because of the complex geometry of the Sung - generally independent lobes often separated large distances by non-lung tissue - the variogram analysis was limited to a maximum distance d _mca (Keil et al. 2012), which we defined as approximately one half the characteristic diameter of the left lobe. From the control group images, d _max was measured to be about 12 mm. Variogram calculation was executed in Python in ~7 seconds. The octree decomposition improved the speed of the 3D variogram calculation ~700-fold as compared to previous approaches (Keil et al. 2012).

[0079] To examine potential alterations to lung ventilation patterns, volume change maps (maps showing the volume change between PRC and PIV, also known as ventilation maps) were calculated from CT images as described in (Jacob et al. 2013) using a method closely following Yin 201 1. In brief, the PIV image was warped to the PRC image, then voxel-by- voxel volume change was calculated using the determinant of the Jacobian of the resulting warp vector field and the HU values of the original CT images. CoV measurement, octree decomposition, and variogram analysis were performed on the volume change maps as described above.

[0080] We attempted to fit standard variogram models to the resulting data. For the raw CT images, a power- law equation of the form:

v = C + mcf [3]

[0081] was least-squares fit to the data. C is a constant, m is the initial slope, d is the distance between octree cube centers, and a is an exponent that represents how rapidly the octree cubes become increasingly spatially dissimilar. Higher values of a (for a> 1 ) indicate that the variance is increasing at an increasing rate, implying that the lung becomes more heterogeneous on average with increasing distance from any given location. The ventilation maps, on the other hand, were fit to a model of exponential rise to an asymptotic value:

v = S[ l - e.\p(- 3< )1. [4]

[0082] where S is the asymptote (referred to the "sill" in geostatistics), d is the distance between octree cubes, and R (the "range") is the distance at w hich the variance reaches -95% of 5. The range characterizes the local "roughness," with a larger R indicative of smoother variation. At distance R, the image becomes random, with little spatial change in variance. The sill indicates the degree of randomness or variety across the image. To compare measured parameters of the dosed groups against those of the control group, an analy sis of variance was performed w ith a Dunnett post hoc test, using a null hypothesis of a=0.05. The significance of the correlation of dose level to different measured parameters (i.e. dose- response) was tested using a paired t-test, also with a null hypothesis of a=0.05.

Results

[0083] Figure 8 summarizes by measurement type and dose group the PFT results and other physiological measurements. "Physical Measures" refers to body eight and lung weights; "Forced Maneuvers" are PFT results for anesthetized rats, "Spontaneous Breathing" are the WBP results, "Ventilation" are parameters from mechanical ventilation during imaging, and "imaging" are results derived from the CT images at FRC and PiV. All measurements show the mean and σ from n=5 rats, except the previously noted 6 Gy group with n=4. In addition, the volume map of one rat in the 12 Gy group displayed a large anomalous ventilation defect in the right lung and was therefore eliminated from all volume map calculations. Significance levels, as compared to the control group, are indicated in the table. Of the measurements in Figure 8, those that had a significant correlation with dose level were the dry lung weight (p=0.03, r=0.45) and the volume map CoV (p=0.013, r=-0.5 l ). It should be noted that significant functional deteriorations were seen only in the 6 Gy group.

[0084] Figure 9A shows representative variograms. made from the FRC images of two typical rats, one 0 Gy and one 15 Gy. The more rapid rise of the control variogram is indicative of comparatively greater variance at a given distance, which is interpreted as greater heterogeneity. Figure 9B shows representative variograms from volume change maps.

[0085] Figure 10 shows reconstructed variograms based on the mean parameters of each dose group for the CT images at FRC (Figure 10A) and the volume maps (Figure 10B). Distinct differences between the 0 Gy and the 15 Gy groups are visible, with the other groups distributed in between. Also the presence of outliers influences the mean, resulting in some overlap between intermediate dose groups.

[0086] Figure 1 1 shows box plots of the variogram parameters. In Figure 1 1 A, a at FRC becomes increasingly more distinct from the 0 Gy group w ith increasing dose to the point that the 15 Gy group is significantly different from control. Indeed, a correlates significantly with dose (p=0.Q03, r=-0.57). However, the correlation of dose with a is not significant at PIV (p=0.09, r=-0.35). For the volume change maps, S correlated significantly with dose (p=0.008, r=-0.53). However R showed different behavior, apparently varying with dose in a quadratic manner (see Figure 1 1 B). Discussion

[0087] We have reported a rapid, objective method of lung CT image analysis that was better able to detect radiation-induced lung damage than conventional pulmonary function tests or image analysis techniques. In particular, variogram parameters from both the octree- decomposed CT images at FRC and the volume change maps correlated significantly with radiation dose. Results indicated that increased exposure to radiation was associated with reduced heterogeneity in this animal model, perhaps indicative of a widespread damage or repair response across the lung. We hesitate to speculate further about the nature of the reduced heterogeneity; however, this point may be elucidated by studies that incorporate blood and alveolar lavage fluid biomarkers as well as lung histological analysis.

[0088] Prior studies of radiation-induced lung damage in rats have focused on biomarker expression, DNA damage, non-invasive imaging, pulmonary function changes, and mechanical changes (Ward et al. 1993, Novako va-J i resova et al. 2007, Wiegman et al. 2003, Vujaskovic t al 1998, Van Erde et al. 2001, Zhang et al. 2008, Pauluhn et al. 2001, Calveley et al. 2005), often with a lack of consensus. An issue that complicates the body of research is the variety of animal models used, including rat strains, radiation doses, irradiated lung fractions, and elapsed time before data collection. Apparent discrepancies include whether there are significant changes to body weight, respiratory rate, collagen content, or CT density. For example, in this study HU did not correlate with dose, consistent with the findings of Wiegman et al. 2003, but contrary to the results of Vujaskovic et al. 1998. Both of these studies used male albino Wistar rats and employed hemithoracic or partial-lung irradiation. Pauluhn et al. 2001 measured significant differences in several physiological parameters between control and irradiated rats (e.g. body weight, lung volumes, etc.) using a hemithoracic exposure of 20 Gy in Fisher rats. In our study, we observed no significant differences between dose and control groups in most of the measurements made, indicating that the radiation damage was below the threshold of detection using these methods and group sizes.

10089] The dose-response correlation i the variogram analysis detected at PR was not detectable at PIV, suggesting that spatial variations in the lung decrease as the lung expands and fills with air. A decrease in a of the control group between PRC and PIV was found, which supports this hypothesis, since the lower exponent indicates a slower increase in spatial variation of the mean 1 1U values. Thus, inflation results in a loss of information that masks the effects of subtle disease and potentially confounds diagnosis. Indeed. Kauczor et al. 2002 showed that the mea HU values in the human lung correlated better with several lung function parameters at expiration than at inspiration, also implying a loss of information in images acquired at inspiration.

[0090] The correlation with dose of volume change map parameters CoV and S implies that volume change maps may be useful diagnostic tools, in addition to their potential to characterize ventilation defects. However, such maps require multiple computational steps, and they demand the acquisition of 4DCT images at two different known inflation levels. The CT-based variograms, on the other hand, require only a single CT image, and a correlated as well or better with dose than S. These attributes make the CT approach more amenable to a clinical setting than the 4DCT approach.

[0091] A limitation of this study was the small numbers of animals in each group. Although a significant dose-response correlation was found from the variogram parameters, a clear separation of dose groups was not possible due to the range of response within each group; however, this could be also a consequence of closely-spaced dose levels.

[0092] The results suggest that variogram analy sis of octree -decomposed CT images is a rapid, automated approach that may be highly sensitive to subtle lung damage or disease. EXAMPLE 3

10093] In this third example we studied 3D CT images from five healthy volunteers and two diagnosed with mild emphysema. Pulmonary function tests were performed to confirm diagnosis. Figure 12 shows variograms constructed from all the image sets, out to a distance of 100 mm. This is equivalent to d _max . Beyond this distance, changes in variance become random indicating no spatial relationships. The variogram data were fit to a standard power law equation of the form y=m*x ^a , where m is the initial slope and a is the exponent describing curvature. For example, a> \ indicates that the line has upward curvature. β<1 indicates downward curvature, and a=\ is a straight line. As shown in the figure 12, the normal subjects (01 -08) all have upward curvature (a> 1 ), whereas the emphysematous subjects (09-10) have a<\ .

10094] To better illustrate this point, Figure 13 is a box plot showing the group medians of the exponent a. The clear separation of the normal subjects (upward curvature) from the emphysematous subjects (flat or downward curvature) suggests that the variogram is able to identify abnormalities present in human lung due to disease. Furthermore, it should be pointed out that in the Figure 12 subject 08 appears to have a highly abnormal variogram, rising sharply with considerable separation from the other normal subjects. However, after fitting with the power fit equation, it was seen that a for this subject is consistent with the other normal subjects (evidenced in the Figure 13 by the lack of outliers in the "normal" box plot), suggesting robustness of the variogram approach across subjects.

Conclusion

[0095] Results of this pre-elinical study of elastase-treated rats suggests that automated octree decomposition and variogram analysis based on image heterogeneity may provide a non-objective and sensitive metric for characterizing emphysematous lung and other diseases, even in early disease stages. The method outperformed conventional approaches that utilize thresholding and absolute HU values. This approach may be applicable to human datasets and other diseases.

References

[0096] West JB: Pulmonary Pathophysiology, 5th edn. Philadelphia: Lippincott Williams

& Wilkins; 1998.

[0097] Mets ONI de Jong PA. van Ginneken B, Gietema HA, Lammers JW: Quantitative computed tomography in COPD: possibilities and limitations. Lung 2012. 190(2): 133-145.

[0098] Robertson HT, Buxton RB: Imaging for lung physiology: what do we wish we could measure? J Appl Physiol 2012, 1 13(2):317-327.

[0099] Pauwels RA, Buist AS, Calverley PM Jenkins CR, Hurd SS: Global strategy for the diagnosis, management, and prevention of chronic obstructive pulmonary disease. NHLBI/WHO Global initiative for Chronic Obstructive Lung Disease (GOLD) Workshop summary. Am J Respir Crit Care Med 2001 , 163(5): 1256- 1276.

[00100] Kohler D, Fischer J, Raschke F. Schonhofer B: Usefulness of GOLD classification of COPD severity. Thorax 2003. 58(9):825.

[00101] Pescarolo M, Sverzellati N, Verduri A, Chetta A, Marangio E, De Filippo M, Oiivieri D, Zompatori M: How much do GOLD stages reflect CT abnormalities in COPD patients? Radiol Med 2008, 1 13(6): 817-829.

[00102] Litmanovich D, Boiselle PM, Bankier AA: CT of pulmonary emphysema— current status, challenges, and future directions. Eur Radiol 2009, 19(3):537-55 1 .

[00103] Chong D, Brown MS, Kim HJ, van Rikxoort EM, Guzman L, Mc itt-Gray MF, Khatonabadi M, Galperin-Aizenberg M, Coy H, Yang K et al: Reproducibility of volume and densitometric measures of emphy sema on repeat computed tomography with an interval of 1 week. Eur Radiol 2012. 22(2):287-294. [00104] Uppaluri R, Mitsa T, Sonka M, Hoffman EA, McLennan G: Quantification of pulmonary emphysema from lung computed tomography images. Am J Respir Crit Care Med

1997, 156( 0:248-254.

[00105) Besir FH, Mahmutyazicioglu K, Aydin L. Altin R. Asil K, Gundogdu S: The benefit of expiratory-phase quantitative CT densitometry in the early diagnosis of chronic obstructive pulmonary disease. Diagn Interv Radiol 2012, 18(3):248-254.

[00106] Irion KL, Marchiori E, Hochhegger B, Porto Nda S, Moreira Jda S, Anselmi CE,

Holemans J A, Irion PO: CT quantification of emphysema in young subjects with no recognizable chest disease. AJR Am J Roentgenol 2009. 192(3):W90-96.

[00107] Reske AW, Busse H, Amato MB, Jaekel M, Kahn T. Schwarzkopf P, Schreiter D,

Gottschaidt U, Seiwerts M: Image reconstruction affects computer tomographic assessment of lung hyperinflation. Intensive Care Med 2008, 34( 1 1 ):2044-2053.

[00108] Yuan R, Mayo JR, Hogg JC, Pare PD, McWilliams AM, Lam S, Coxson HO: The effects o radiation dose and CT manufacturer on measurements of lung densitometry. Chest 2007, 132(2):617-623.

[00109] Chae EJ, Seo JB, Song JW, Kim N, Park BW, Lee YK, Oh YM, Lee SD, Lim SY: Slope of emphysema index: an objective descriptor of regional heterogeneity of emphysema and an independent determinant of pulmonary function. AJR Am J Roentgenol 2010. 194(3):W248-255.

[00110] Boehm HF, Fink C. Attenberger U, Becker C, Behr J, Reiser M: Automated classification of normal and pathologic pulmonary tissue by topological texture features extracted from multi-detector CT in 3D. Eur Radiol 2008. 18(12):2745-2755.

[00111 ] Xu Y. Sonka M. McLennan G, Guo J, Hoffman EA: MDCT-based 3-D texture classification of emphysema and early smoking related lung pathologies. IEEE Trans Med Imaging 2006, 25(4):464-475. [001 12] Hersh CP, Washko OR, Jacobson FL, Gill R. Estepar RS, Re illy JJ, Silverman

EK: Interobserver variability in the determination of upper lobe-predominant emphysema.

Chest 2007, 131(2):424-431.

[00113] Mascalchi M, Diciotti S, Sverzellati N, Camiciottoli G. Ciccotosto C, Faiaschi F, Zompatori M: Low agreement of visual rating for detailed quantification of pulmonary emphysema in whole-lung CT. Acta Radiol 2012, 53(1 ):53-60.

[00114] Fain SB, Panth SR, Evans MD, Wentland AL, Holmes JH, Korosec FR, O'Brien MJ, Fountaine H, Grist TM: Early emphysematous changes in asymptomatic smokers: detection with 3 He MR imaging. Radiology 2006, 239(3):875-883.

[00115] Quirk JD, Lutey BA, Gierada DS, Woods JC, Senior RM. Lefrak SS, Sukstanskii AL, Conradi MS, Yablonskiy DA: In vivo detection of acinar mierostructural changes in early emphysema with (3)He lung morphometry. Radiology 201 1, 260( 3}:866-874.

|00116] Afford SK, van Beek EJ, McLennan G, Hoffman ΕΛ: Heterogeneity of pulmonary perfusion as a mechanistic image-based phenotype in emphysema susceptible smokers. Proc Natl Acad Sci U S A 2010, 107(16):7485-7490.

[00117] Vidal Melo MF, Winkler T, Harris RS, Musch G, Greene RE, Venegas JG: Spatial heterogeneity of lung perfusion assessed with (13)N PET as a vascular biomarker in chronic obstructive pulmonary disease. J Nucl Med 2010, 51(1 ):57-65.

[00118] Vogt B, Puiletz S, Elke G, Zhao Z, Zabel P, Weiler N, Frerichs I: Spatial and temporal heterogeneity of regional lung ventilation determined by electrical impedance tomography during pulmonary function testing. J Appl Physiol 2012, 113(7): 1 154- 1 161 .

[00119] Chen CC, Daponte JS, Fox MD: Fractal feature analysis and classification in medical imaging. IEEE Trans Med Imaging 1989, 8(2): 133- 142.

[00120] Gienny RW, Robertson FIT: Fractal properties of pulmonary blood flow: characterization of spatial heterogeneity. J Appl Phy siol 1990, 69(2):532-545. [00121 ] Copley SJ, Giannarou S, Schmid VJ, Hansell DM, Wells AU, Yang GZ: Effect of

Aging on Lung Structure In Vivo: Assessment With Densitometric and Fractal Analysis of

High-resolution Computed Tomography Data. J Thorac Imaging 2012, 27(6):366-371 .

[00122] Kido S, Ikezoe J, Naito H, Tamura S, Machi S: Fractal analysis of interstitial lung abnormalities in chest radiography. Radiographics 1995. 15(6): 1457-1464.

[00123] Uppaluri R, Mitsa T, Gaivin JR: Fractal analysis of high-resolution CT images as a tool for quantification of lung disease. In: Medical Imaging 1995: Physiology and Function from Multidimensional Images. Volume 2433, edn. Edited by Hoffman EA. Beil ingham,

WA: SPIE; 1995: 133- 142.

[00124] Kido S, Sasaki S: Fractal analysis for quantitative evaluation of diffuse lung abnormalities on chest radiographs: use of sub-ROls. J Thorac Imaging 2003, 18(4):237-241 .

[00125] Subramaniam K. Hoffman EA, Tawhai MH: Quantifying Tissue Heterogeneity using Quadtree Decomposition. In: 34th Annual International Conference of the IEEE EMBS. San Diego, CA, USA; 2012: 4079-4082.

[00126] Shephard MS, Georges MK: Automatic 3-Dimensional Mesh Generation by the Finite Octree Technique. Int J Numer Meth Eng 1991. 32(4):709-749.

[00127] Zhang JH, Owen CB: Octree-based animated geometry compression. Comput Graph-Uk 2007, 3 1 (3):463-479.

[00128] Szeliski R. Lava! lee S: Matching 3-D anatomical surfaces with non-rigid deformations using octree-splines. Int J Comput Vision 1996, 1 8(2): 171-186.

[00129] Zhang Yj, Bajaj C. Sohn BS: 3D finite element meshing from imaging data. Comput Method Appl M 2005. 194(48-49):5083-5106.

[00130] Dua S, Kandiraju N, Chowriappa P: Region quad-tree decomposition based edge detection for medical images. Open Med Inform J 2010, 4:50-57.

[00131] Clark I, Harper WV: Practical Geostatistics 2000. Columbus, Ohio: Ecosse North America, LLC; 2000. [00132] Gringarten E, Deutsch CV: Variogram interpretation and modeling. Math Geol

2001 , 33(4):507-534.

[00133] Keii F, Oros-Peusquens AM, Shah NJ: Investigation of the spatial correlation in human white matter and the influence of age using 3 -dimensional variography applied to MP- RAGE data. Neuroimage 2012, 63(3): 1374- 1383.

[00134] Jacob RE, Minard KR, Laicher G, Timchalk C: 3D 3 He diffusion MR! as a local in vivo morphometric tool to evaluate emphysematous rat lungs. J Appl Physiol 2008, 105(4): 1291-1300.

[00135] Jacob RE, Lamm WJ: Stable small animal ventilation for dynamic lung imaging to support computational fluid dynamics models. PLoS One 201 1 , 6(1 l ):e27577.

[00136] Corley RA, Kabilan S, Kuprat AP, Carson JP, Minard KR, Jacob RE, Timchalk C,

Glenny R, Pipavath S, Cox T et al: Comparative computational modeling of airflows and vapor dosimetry in the respiratory tracts of rat, monkey, and human. Toxicol Sci 2012,

128(2):500-5 16.

[00137] Minard KR, Kuprat AP, Kabilan S, Jacob RE, Einstein DR, Carson J P. Corley RA: Phase-contrast MRI and CFD modeling of apparent (3)He gas flow in rat pulmonary airways. J Magn Reson 2012, 221 : 129-138.

[00138] ImageJ [http://imagej.nih.gov/ij/]

[00139] ImageJ Plugins: 3D Toolkit [http://ij-plugins.sourceforge.net/plugins/3d- toolkit/index.html]

[00140] Jackins CL, Tanimoto SL: Oct-trees and their use in representing three dimensional objects. Computer Graphics and Image Processing 1980, 14(3):249-270.

|00141 ] Schulz H, Miihle H: Respiration. In: The Laboratory Rat. edn. Edited by Krinke GJ. San Diego: Academic Press; 2000: 323-336.

[00142] Spencer H: Pathology of the Lung. In. Volume 1, 4th edn. New York: Pergamon Press; 1985: 557-594. [00143] Bohling GC: Introduction to Geostatistics. In.: Kansas Geological Survey Open

File Report no. 2007-26, 50pp.; 2007.

(00144] Hsia CC, Hyde DM. Ochs M, Weibel ER: An official research policy statement of the American Thoracic Society/European Respiratory Society: standards for quantitative assessment of lung structure. Am J Respir Crit Care Med 2010, 181 (4):394-418.

[00145] Jacob RE, Carson JP. Gideon KM, Aniidan BG, Smith CL, Lee KM: Comparison of two quantitative methods of discerning airspace enlargement in smoke-exposed mice. PLoS One 2009, 4(8):e6670.

[00146] Lynch DA, Travis WD, Muller NL, et al. Idiopathic interstitial pneumonias: CT features. Radiology 2005; 236: 10-21.

[00147] Lehnert S, el-Khatib E. The use of CT densitometry in the assessment of radiation-induced damage to the rat lung: a comparison with other endpoints. hit J Radiat Oncol Biol Phys 1989; 16: 1 1 7- 124.

[00148] Ward ER, Hedlund LW, Kurylo WC, et al. Proton and hyperpolarized helium magnetic resonance imaging of radiation-induced lung injury in rats. Int J Radial Oncol Biol Phys 2004;58: 1562-1569.

[00149] Ward HE, Kemsley L, Davies L, et al. The pulmonary response to sublethal thoracic irradiation in the rat. Radiat Res 1993; 136: 1 5-21.

[00150] Quantitative Concepts and dosimetry in radiobiology: Report 30: Internationa! Commission on Radiation Units (1CRU); 1979. pp. 13-14.

[00151] ovakova-J iresova A, van Luijk P, van Goor H, et al. Changes in expression of injury after irradiation of increasing volumes in rat lung. Int J Radiat Oncol Biol Phys 2007;67: 1510-1518.

[00152 ] Rasband WS. ImageJ U.S. National Institutes of Health, http://imagej.nih.gov/ij/.

[00153] Sacha J. Image J Plugins: 3D Toolkit http://ij-plugins.sourceforge.net/piugins/3d- tooikit/index.html. [00154] Jacob RE, Carson JP, Thomas M, et al. Dynamic Multiscale Boundary Conditions for 4D CT of Healthy and Emphysematous Rats PLoS One 2013 :8:e65874.

[00155} Yin Y. MDCT-based dynamic, subject-specific lung models via image registration for CFD-based interrogation of regional lung function. Mechanical Engineering. Vol Ph.D. Iowa City. lA: University of Iowa; 201 1 .

[00156] Wiegman EM, Meertens H, Konings AW, et al. Loco-regional differences in pulmonary function and density after partial rat lung irradiation. Radiother Oncol 2003;69:11 -19.

[00157] Vujaskovic Z, Down JD, van t' Veld A A, et al. Radiological and functional assessment of radiation-induced lung injury in the rat. Exp Lung Res 1998;24: 137-148.

[00158] van Eerde MR, ampinga HH. Szabo BG, et al. Comparison of three rat strains for development of radiation-induced lung injury after hemithoracic irradiation. Radiother Oncol 2001 :58:313-316.

[00159] Zhang R, Ghosh SN, Zhu D. et al. Structural and functional alterations in the rat lung following whole thoracic irradiation with moderate doses: injury and recovery. Int J Radial Biol 2008;84:487-497.

[00160] Pauluhn J, Baumann M, Hirth-Dietrich C, et al. Rat model of lung fibrosis: comparison of functional, biochemical, and histopathological changes 4 months after single irradiation of the right hemithorax. Toxicology 2001 : 161 : 153- 163.

[00161] Calveley VL, Khan MA. Yeung I W, et al. Partial volume rat lung irradiation: temporal fluctuations of in-field and out-of-field DNA damage and inflammatory cytokines following irradiation. Int J Radiat Biol 2005;81 :887-899.

[00162] Kauczor HU, Hast J, Heussel CP, et al. CT attenuation of paired HRCT scans obtained at full inspiratory/expiratory position: comparison with pulmonary function tests. Eur Radiol 2002;12:2757-2763. [00163] The present invention has been described in terms of spec ific embodiments incorporating details to facilitate the understanding of the principles of construction and operation of the invention. As such, references herein to specific embodiments and details thereof are not intended to limit the scope of the claims appended hereto. It will be apparent to those skilled in the art that modifications can be made in the embodiments chosen for illustration without departing from the spirit and scope of the invention.