INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

THE FRACTAL DIMENSION (D) of a self-similar spatial structure is often estimated by calculating spatial heterogeneity as a function of the scale of spatial resolution. The present study proposes a new method for estimating D directly from a functional image. In general, spatial heterogeneity of an image of a measured quantity X, which represents a physiological function, e.g., pulmonary blood flow, is a power law function of the spatial resolution scale R used to make the measurement (2, 3, 7, 9, 19, 25). If the spatial distribution of the measured quantity is fractal, then a relationship exists between X(R) and some power of R used to measure that variable, i.e.

(1)

where R₀ is an arbitrary resolution scale used as a reference point. A parameter used commonly to quantify spatial heterogeneity is the coefficient of variation (cov), defined as the standard deviation of the data () normalized by its mean or cov = /

Pulmonary perfusion heterogeneity has been measured from the distribution of intravenously injected radioactive-labeled or fluorescent microspheres in the pulmonary circulation (8, 10). Microspheres (15 µm diameter) are thought to deposit in capillaries of the lung in direct proportion to local blood flow. The animal is then killed, and the lungs are excised, inflated to total lung capacity, dried, and sliced into cubic pieces. The coefficient of variation is calculated for increasingly larger resolution scales by grouping adjacent pieces into progressively larger cubes. A limitation of this method is that increasing cube size requires the exclusion of peripheral data, which would incorporate pieces outside the lung, or the inclusion of cubes of smaller size, which biases the measurements. Furthermore, the evaluation of D via this method may be subject to misregistration artifacts between the measuring grid and the physiological structure: the grid and the organ do not have the same geometric shape, may not be an integer multiple size of the fractal pattern, or may have an offset relative to the fractal pattern.

We investigated an alternative method to assess the D of functional heterogeneity of tomographic data sets. We theorized that two-dimensional low-pass filtering a functional image could be used to assess its D, given that the corner frequency (f_c) of the low-pass filter should be equivalent to the inverse of the resolution scale. The f_c of an ideal low-pass filter is the frequency above which all energy is filtered out of the signal.¹ An advantage of this method would be that no data would need to be excluded from the analysis. In this study the low-pass-filtering method was evaluated in three types of images: random noise images in the shape of a tomographic lung cross section, synthetically generated fractal images, and positron emission tomographic (PET) images of pulmonary blood flow.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Measurement of D

(2)

(3)

(4)

Generation of Two-Dimensional Data Sets

Random noise images. Using a pseudorandom number generator, noise images were created with nonzero voxels forming the shape and size of a typical dog lung cross section. Voxel elements were chosen from a normal distribution with a mean and variance equal to 1.0. For each noise image, a log-log plot of cov(f_c) vs. 1/f_c was generated, and D was determined from the slope of the least squares regression line (Eq. 4). One hundred different random noise images were generated and analyzed with the low-pass-filtering method.

Synthetic fractal images. Synthetic fractal images were used to compare the results obtained using the conventional method of estimating D with results obtained using the two-dimensional low-pass-filtering method. To create synthetic fractal images, we first defined a 2 × 2-similarity matrix made of four nonequal positive elements summing to unity. Fractal images were generated by dividing a square matrix of 128 × 128 unity values into four identical squares and multiplying the voxels in each of these by the corresponding element of the 2 × 2-similarity matrix. Each of these four squares was further subdivided into four, and each was multiplied again by the corresponding elements in the 2 × 2-similarity matrix. The process was repeated recurrently until the smallest subdivision was the size of the similarity matrix (Fig. 1). The conventional method of estimating D was simulated by successively grouping data into adjacent squares of equal size and calculating the coefficient of variation of the corresponding averages. The process was repeated for larger and larger grouping squares until the fractal image consisted of four squares. The resolution scale was taken as the side length in number of voxels of the grouping square.

View larger version (180K): [in this window] [in a new window]   Fig. 1.   Square synthetic fractal image (128 × 128 voxels) used to compare fractal dimension (D) by use of conventional and low-pass-filtering methods.

PET images of pulmonary perfusion. A PET camera was used to image pulmonary perfusion in six anesthetized and mechanically ventilated supine dogs with normal lungs, as approved by the Massachusetts General Hospital Subcommittee on Animal Care. The method used to assess local pulmonary perfusion has been described in detail elsewhere (11, 12). Briefly, the method consists of injecting a bolus of ¹³NN-labeled saline into a jugular vein catheter at the initiation of a period of apnea. Infusion time ranged from 10 to 15 s, depending on the specific activity of the infusate (100-200 µCi/ml) to produce images with consistent levels of radioactivity per voxel. Simultaneously with the initiation of injection, a series of six sequential PET images (each 10 s in duration) was initiated. Because of the low solubility of ¹³NN gas in blood and tissues (partition coefficient = 0.018) in normal aerated lungs, virtually all radioactivity is transferred to alveolar gas at first pass, and regional tracer content is proportional to local perfusion (16, 18). At the end of the 60 s of apnea, mechanical ventilation was restarted, and the tracer was washed out.

Image Processing and Analysis

A significant portion of the image heterogeneity is caused by a gradient in the vertical direction, which is superimposed to the expected fractal pattern. To remove this gradient from the images, the best-fit vertical gradient plane was obtained by linear regression and subtracted from the original image. The resulting images represented residual spatial heterogeneity and were also analyzed with the low-pass-filtering method. Surface plots of a representative perfusion image, its vertical gradient plane, and the residual image are shown in Fig. 2.

View larger version (26K): [in this window] [in a new window]   Fig. 2.   Cross-sectional positron emission tomographic (PET) image of pulmonary perfusion presented as an illuminated surface plot, where height above x-y plane represents relative magnitude of regional perfusion (A). Gradient plane (B) was subtracted from original image, resulting in residual image (C) used in fractal analysis.

Filter Design

1

1

Edge Effects

View larger version (27K): [in this window] [in a new window]   Fig. 3.   Lung mask (A), low-pass-filtered lung mask (B), low-pass-filtered PET image in Fig. 2 after low-pass filtering (C), and edge-corrected filtered image resulting from normalizing C by B (D).

	RESULTS

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

Random Noise Images

View larger version (11K): [in this window] [in a new window]   Fig. 4.   Fractal plot for low-pass-filtering method of natural logarithm of coefficient of variation (cov) vs. natural logarithm of inverse corner frequency (fc) of filter normalized by reference frequency (fc,0). Values are means ± SD of 100 random noise images. Linear regression fit through mean data yields D = 1.5045 with R2 = 0.9998.

Synthetic Fractal Images

View larger version (74K): [in this window] [in a new window]   Fig. 5.   Synthetic fractal image in Fig. 1 processed with low-pass-filtering method (A) and conventional method (B) with increasing resolution scale (top to bottom). For low-pass-filtering method, image was filtered with fc of 0.492 fN, 0.2265 fN, and 0.0498 fN, respectively, where fN is Nyquist frequency. For conventional method, side lengths of grouping squares in number of voxels were 4, 8, and 32 voxels, respectively.

View larger version (14K): [in this window] [in a new window]   Fig. 6.   Simulated conventional method fractal plots of natural logarithm of coefficient of variation expressed as percentage [ln(cov%)] vs. natural logarithm of square root of number of voxels per square piece (nv) applied to synthetic fractal image in Fig. 1. Each plot corresponds to analysis conducted with proper grid, misregistered grid, and misaligned grid.

View larger version (12K): [in this window] [in a new window]   Fig. 7.   Low-pass-filtering method fractal plots of natural logarithm of coefficient of variation expressed as ln(cov%) vs. natural logarithm of inverse fc of filter normalized by a reference frequency {ln[(1/fc)/(1/fc,0)]}. Method was applied to synthetic fractal image in Fig. 1 and to same image after periphery voxels were padded with zeros. Best-fit lines to both data sets coincide exactly.

Analysis of PET Images of Pulmonary Perfusion

View larger version (16K): [in this window] [in a new window]   Fig. 8.   Low-pass-filtering method fractal plots of ln(cov%) vs. ln[(1/fc)/(1/fc,0)]. Method was applied to PET images of pulmonary perfusion of 6 normal supine dogs. A plateau is reached at same spatial frequency for all dogs, which corresponds to intrinsic spatial resolution of PET camera. For reference, corresponding length scales are presented below x-axis.

The ln(cov%) vs. {ln[(1/f_c)/(1/f_c,0)]} for the same perfusion images not detrended for vertical gradients followed a linear relationship for 0.28 f_N > f_c > 0.02 f_N with average R² = 0.96. D ranged from 1.11 to 1.14 (1.12 ± 0.01). The fact that the vertical gradient contributes low-frequency energy across the frequency spectrum studied is reflected by higher cov% at every resolution length for the uncorrected images than for the corresponding vertical gradient detrended images (Fig. 9).

View larger version (16K): [in this window] [in a new window]   Fig. 9.   Low-pass filtering method fractal plots of ln(cov%) vs. ln[(1/fc)/(1/fc,0)]. Method was applied to a PET image from dog 3 in Fig. 8 with (R2 = 0.9869) and without (R2 = 0.9435) detrending for vertical gradient in perfusion. For reference, corresponding length scales are presented below x-axis.

	DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS DISCUSSION REFERENCES

The most important findings of this study are as follows. 1) Successive low-pass filtering of two-dimensional functional images (e.g., PET images of blood flow) can be used to assess their D. 2) The edge-smearing artifact inherently created by low-pass-filtering images can be virtually eliminated by normalization with equivalently filtered unity masks. 3) Fractal analysis of residual spatial heterogeneity in pulmonary perfusion after removal of the vertical gradient yielded linear log-log plots (average R² = 0.99). 4) In contrast, analysis of the same perfusion images without removal of the vertical gradients reduced the goodness of fit of the log-log plots to the linear model and systematically lowered their estimated D.

To our knowledge, this is the first application of low-pass filtering to derive a D of one- or two-dimensional data sets. Different types of D analysis have been used to provide a measure of texture of two-dimensional radiographic images of bone (4, 14), lungs (21), and breast (15), of brain single photon emission tomographic images (13), and of electron micrograph images (6). These methods have not been applied to characterize the regional distribution of blood flow images. There are a number of approaches to define a D from images, including Fourier decomposition, average intensity difference of pixel pairs, and the reticular cell-counting approach (6). The Fourier decomposition method is applied to one-dimensional samples of the image (21) or to the two-dimensional data set in polar coordinates yielding power spectrum plots (14). The Fourier method is, therefore, equivalent to the low-pass-filtering method, in that it provides a measure of the heterogeneity content of the image at the different spatial frequencies. Artifacts produced by edge effects and power introduced in the frequency domain by the shape of the object, however, limit the use of Fourier analysis to small areas of the organ away from the edges.

Applicability of the low-pass-filtering method to functional image data is based on the inverse proportionality between the resolution scale and f_c. The conventional technique of assessing fractal characteristics varies the resolution scale by changing the number of adjacent cubic tissue pieces that are averaged before estimation of the heterogeneity. With the low-pass-filtering method, a change in the resolution scale is achieved by changing f_c of the low-pass filter. In general, the low-pass filters operate by removing high-frequency heterogeneity from the images (spatial frequency components above f_c), thereby creating a blurring effect. This blurring effect is analogous to the effect created by averaging adjacent groups of pieces in the conventional method (Fig. 5).

The low-pass-filtering method overcomes two inherent disadvantages associated with the conventional method: 1) it is insensitive to misregistration artifacts, and 2) it includes all the image data in the analysis at all resolution scales. These advantages were demonstrated in the analysis of synthetic fractal images. These images can be thought of as idealized two-dimensional images of flow distribution in a branching network where each parent vessel distributes flow to four daughter branches in proportions determined by the self-similarity matrix. In that sense, these images are equivalent to the dichotomous branching model proposed by Glenny and Robertson (8) to simulate the fractal properties of pulmonary blood flow distribution. We used four instead of two subdivisions at each parent-to-daughter generation to simplify the computation in generating and analyzing the two-dimensional image.

Analysis of these images demonstrated that the low-pass-filtering method is insensitive to misregistration effects, yielding values of D that were not changed when zero padding was added to the images. The conventional method yielded log-log plots with reduced goodness of fit to the linear model and D that changed with zero padding. This is because the conventional method imposes an orthogonal coordinate system with an arbitrarily selected location for the origin to divide the image into cubic pieces. As demonstrated, if the imposed coordinate axes are not in alignment with the geometrical pattern governing the fractal distribution, this misalignment causes the grouping of erroneous adjacent units. Glenny and Robertson (8) quantified changes in D, inasmuch as this type of misalignment was artificially forced on data sets created using a branching dichotomous model. They found a maximum 2.5% error in the measured D compared with the true D when the partitioning grid was offset by one-half the length of the perfusion unit. The true D refers to that calculated when the imposed partitioning boundaries corresponded to the boundaries of the fractal partitioning. The small error, however, stems in part from the assumed square shape. Greater errors could be expected in the estimation of D with misalignment when noncubic peripheral pieces are considered. We tested this hypothesis in the synthetic fractal images by introducing a misalignment of one voxel in each topological dimension between the measuring grid and the image by padding one row and one column with zeros. Groups of voxels including zeros were not included in the fractal analysis. This meant that, as the resolution scale was increased, a greater fraction of voxels was excluded. Calculated D from this analysis was ~5.5% lower than the true D. Thus combined misalignment and exclusion of peripheral data increase the error in estimation of D. In contrast, the low-pass-filtering method yielded equal values of D, despite zero padding, demonstrating the insensitivity and robustness of the low-pass-filtering method. A second type of misregistration, where the geometry of the grid does not correspond to that of the fractal pattern, was simulated by defining the size of the initial measuring grid, or 3 × 3 voxels, instead of the fractal pattern (2 × 2 voxels) used to generate the image. This case resulted in an estimation error for D of 1%. Because the low-pass-filtering method does not require the definition of an initial grid, the method is independent of the geometry, size, or shape of the fundamental fractal pattern.

Another advantage of the low-pass-filtering method over the conventional method is the fact that there is no need to exclude peripheral data in the calculation of heterogeneity. In the conventional method, as data are grouped in larger cubes, peripheral data have to be excluded to avoid including in these cubes voxels outside the organ boundary. In fact, even at the smallest cube size, peripheral pieces are typically smaller than the rest of the cubic pieces that comprise the lung. To correct for errors in the blood flow measurement caused by nonuniform piece size, the data are normalized by piece dry weight, and spatial heterogeneity is then calculated from data representing blood flow per unit of dry tissue weight. However, because of the noncubic shape of the lung, peripheral pieces tend to be smaller than internal pieces, and thus they contribute with the higher heterogeneity expected at smaller resolution scales. If the distribution of blood flow is fractal, then smaller resolution scales should contribute with increased spatial heterogeneity and result in overestimation of the total heterogeneity. Also, peripheral pieces are given the same statistical "weight" as the cubic pieces in the calculation of spatial heterogeneity, giving additional bias in the calculation.

The conventional method may also require excluding data from the fractal plot when fitting it with linear regression, whereas the low-pass-filtering method does not. Increasing piece sizes causes the coefficient of variation to decrease more rapidly than predicted by the fractal model of Eq. 3. Caruthers and Harris (5) experienced the need to exclude fractal plot data because of digression from linearity. The point at which their fractal plots digressed was chosen by inspection and typically occurred for resolution scales corresponding to 10 pieces of lung. Van Beek (22) theorized that, as too few pieces are used in the calculation of heterogeneity, the fractal plot diverges from the theoretical line because the flow heterogeneity (and self-similarity) is obscured as flow is considered in very large amounts. An alternative explanation for that observation has been proposed on the basis of the argument that at large piece sizes the fractal algorithm combines pieces that are further away from each other (8), and because blood flow to those pieces tends to be negatively correlated, by conservation of blood flow, large pieces include regions of high and low flow, thus decreasing heterogeneity more than predicted by the fractal model. However, application of the conventional method to our synthetic fractal images, which have no restrictions of conservation of flow (Fig. 1), also showed a departure from linearity at the large resolution scales. In contrast, for the same image analyzed with the low-pass-filtering method, fractal plot data maintained linearity up to the highest resolution scales (Fig. 7).

Of critical importance to the low-pass-filtering method is correction for edge effects. Without this correction, an artifact is created when images are low-pass filtered. That is, perfusion values at the edges of the lung are artifactually reduced. This effect extends deeper into the lung as the filter's f_c is lowered. The effect of edge smearing in the fractal analysis is to increase the calculated heterogeneity, resulting in overestimation of the coefficient of variation for a given resolution. The artifactual smearing of the edges in the low-pass-filtering method could be thought of as being equivalent to the conventional method, in that it would require that peripheral data be excluded from the analysis at increasing resolution scale.

The rationale for our correction method follows an observation noted in a study of the noise characteristics of PET images (23). In that study, it was noted that artifacts from image reconstruction affected edges of images of the same lung in equal proportion and that the ratio of two images similarly processed had edge effects canceled. In our analysis of perfusion images, edge effects induced by the low-pass-filtering process are proportionally equivalent on the image and mask. Thus normalization of the former by the latter eliminates edge smearing and returns their original values. The effectiveness of the method is best exemplified in an image of random noise (Fig. 10). With correction for edge smearing, there was an excellent fit to the linear model (R² = 0.9997), and appropriate D for noise (D = 1.53) was obtained from the regression analysis on the fractal plots. Without correction, heterogeneity increases with increasing 1/f_c, resulting in considerable deviation from linearity (R² = 0.90) in the fractal plot and an erroneous D = 1.25. 

View larger version (13K): [in this window] [in a new window]   Fig. 10.   Low-pass-filtering method fractal plots of ln(cov%) vs. ln[(1/fc)/(1/fc,0)]. Method was applied to a pseudorandom noise image with and without normalization by an equivalently filtered unity mask to correct for filtering edge effects.

Before conducting the fractal analysis, we removed any vertical gradient present in the image. The effects of systematic gradients superimposed onto one-dimensional data sets have been recognized to cause artifacts in the fractal analysis. As in our image analysis method, these artifacts are customarily eliminated by "detrending" the data, a process that consists of subtracting a first- or second-order polynomial from the data before the fractal analysis is conducted (21).

D values, obtained for perfusion PET images of normal supine dogs with the low-pass-filtering method in this study (D = 1.31 ± 0.04), are somewhat higher than those reported in the literature. D values range from 1.09 ± 0.02 for supine dogs (8) to 1.14 ± 0.09 for supine sheep (5). None of the literature mentions detrending the data for vertical gradients in the lung perfusion data sets. Glenny and Robertson (8) observed dorsoventral perfusion gradients. Barman et al. (1) studied the effects of pulmonary blood flow rates on D and observed cephalocaudal perfusion gradients. According to their observations, the lung lobe under investigation was unevenly perfused irrespective of cardiac output, in that the caudal regions received the highest blood flow, whereas the cephalic regions near the hilum exhibited the lowest blood flow. Parker et al. (17) reported significant lung perfusion gradients in the unanesthetized canine left lung: 4.7%/cm decreasing gradient running in the gravity-dependent direction, 7.2%/cm gradient running radially outward from the lung midpoint toward its periphery, and 2.5%/cm gradient running from the base to the apex.

Reanalysis of our six supine images with the dorsoventral gradients left intact yielded D = 1.12 ± 0.01 with R² = 0.96. These values are comparable to all those in aforementioned literature (Table 1), but the goodness of fit of the fractal plot data to the linear model was reduced. As expected, at all resolution lengths, the (cov)² values of the uncorrected images were higher than those of the detrended images because of the low-frequency energy introduced by the vertical gradient. Furthermore, as shown in Fig. 9, the fractal analysis conducted in images without removal of the vertical gradients resulted in fractal plots that deviated from the linear model more than from the corresponding detrended images. In contrast, analysis of the same images after removal of the vertical gradients resulted in plots with less deviation from the linear model.

Understanding the effect of vertical gradients on the calculation of D can be facilitated by a physical interpretation of the parameter. From Eq. 3, it is apparent that D is related to the exponential reduction in the (cov)² as length scale is increased and thus reflects the relative contribution of the different length scales to the spatial heterogeneity of the variable. A vertical gradient in perfusion adds (cov)² at length scales larger than the dimension of the lung. Because their contribution to (cov)² will be constant in all low-pass-filtered images, they should weight D toward 1.0, i.e., the value of D for perfect spatial homogeneity or perfectly uniform flow. The deviation from the linear model of fractal plots from images including vertical gradients caused a reduction in R² from 0.99 to 0.96 and suggests that vertical gradients in perfusion (and possibly other large-scale gradients) 1) are not fractal in nature, 2) degrade the goodness of fit to the fractal model, and 3) may be superimposed to a fractal distribution pattern of blood flow.

It is not entirely clear why D of a two-dimensional coronal slice should be the same as that of a volume. In fact, in the presence of substantial craniocaudal gradients, a coronal slice through a fractal volume could be more homogenous than the volumetric set. We tested this hypothesis in each of the dogs studied by comparing the (cov)² of the basal slice analyzed by the low-pass-filtering method with that of the aggregate of the 15 slices simultaneously acquired with the PET camera covering ~70% of the dog lungs. From this test, we found that the (cov)² of the slice analyzed was only 15.8 ± 4.5% lower than that of the volumetric data for four of the six images studied and that the (cov)² of the volumetric data was only 1.4 and 2.0% lower than that of the slice analyzed for the other two dogs studied. This suggests that the (cov)² of the slice analyzed is a reasonable estimate of the (cov)² of the lung as a whole, thus explaining the similarity of D between the microsphere data and the low-pass-filter-analyzed PET data without detrending for vertical gradients.

Finally, a plateau was observed in our fractal plot for all six dogs at a low-pass f_c corresponding to the intrinsic spatial resolution of the PET camera (7 mm). In none of the reported microsphere data was a plateau observed in their fractal plots. Fractal analysis has been proposed as a means to identify the functional unit of perfusion by finding the resolution scale where the coefficient of variation inflects toward a plateau. Theoretically, the relationship between spatial heterogeneity and the resolution scale will remain linear on a log-log scale as long as the heterogeneity of blood flow distribution is fractal. As the anatomic limit of the capillary is reached, flow heterogeneity will cease to increase, causing the fractal plot to bend toward a plateau and stabilize. Glenny and Robertson (8) did not observe this bend or stabilization for tissue pieces as small as 24 mm³, corresponding to a side length of ~2.9 mm. This suggests that the resolution scale of the functional unit of perfusion could be <2.9 mm. The spatial resolution in PET is limited by the camera's geometry and the range of the positron before annihilation with an electron (~3-4 mm for positrons emitted from ¹³NN in lung tissue of density 0.5-0.7 g/cm³). Thus the point of plateau in the fractal plot corresponds to the resolution limit of the imaging method.

In summary, we have devised an alternative method to characterize fractal properties of functional images. Successive two-dimensional low-pass filtering of the images effectively changes spatial resolution and thus their heterogeneity. The resolution scale is inversely proportional to the f_c of the low-pass filter. D can be computed from the slope of the log-log plot of the coefficient of variation vs. the inverse of the low-pass-filtered f_c. The method does not require the exclusion of peripheral perfusion data with increasing resolution scale and, therefore, ensures a robust measure of D. Moreover, the method is insensitive to the effects of misregistration and lung geometry.