An improved technique for fractal characterization called the modified blanket method is introduced that can quantify surrounding fractal structures on a pixel by pixel basis without artifacts associated with scale-dependent image features such as object size. The method interprets images as topographical maps, obtaining information regarding the local surface area as a function of image resolution. Local fractal dimension (FD) can be quantified from the power law exponent derived from the surface area and image resolution relationship. We apply this technique on simulated cell images of known FD and compared the obtained values to power spectral density (PSD) analysis. Our method is sensitive to a wider FD range (2.0–4.5), having a mean error of 1.4% compared to 6% for PSD analysis. This increased sensitivity and an ability to compute regional FD properties enabled the discrimination of the differences in radiation resistant cancer cell responses that could not be detected using PSD analysis.
© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
Fractal patterns are very common in nature and can be observed in different types of biological structures and functions [1,2]. Image analysis of these self-repeating patterns over length- or time-scales is an attractive approach to quantify changes in biological structures . For example, mitochondria undergo fission and fusion based on the metabolic needs of the cell , and it has been demonstrated that mitochondrial organization follows the statistical properties of self-similar fractals [5–8]. Previous studies have demonstrated that mitochondrial reorganization occurs in a wide variety of pathological conditions, including Parkinson’s disease [5,9], cancer [4,6,10–15], and mitochondrial diseases . However, there are several challenges in quantifying the fractal organization of organelles within cells.
Various techniques have been developed to estimate fractal dimension. One of the oldest and simplest approaches for fractal analysis is known as the box counting method, which superimposes boxes of decreasing size over a region of interest within an image . The power law relationship between the number of boxes intersected by the pattern of interest and the size of the boxes is used to estimate fractal organization [1,17]. However, box counting can only be performed on binary images and is most commonly used to characterize the border of objects. To quantify the fractal dimension of grayscale images, Fourier-based approaches have been used, particularly in mitochondrial clustering analysis. Through radial sampling of two-dimensional power spectral density (PSD) maps, fractal dimension (FD) can be determined from images by measuring the exponent, β, from the power law relationship between PSD and spatial frequency as shown in Eq. (1) .10] calculated the power exponent β from NADH autofluorescence images as a measure of mitochondrial clustering and identified depth-dependent changes between normal and precancerous epithelial tissue. Pouli et al.  also utilized PSD analysis to evaluate mitochondrial dynamics from NADH intensity images in different epidermal layers of human skin. Similarly, Liu et al.  used PSD analysis of NADH images to measure mitochondrial clustering in relation to structural metabolic changes caused by controlled perturbations such as extrinsic and intrinsic mitochondrial uncoupling, glycolysis, and fatty acid oxidation. However, this PSD approach is susceptible to errors when cell segmentation is inaccurate and scale-dependent features remain visible after pre-processing. Furthermore, the PSD approach is also only able to obtain the overall FD of an image and cannot provide region-specific information, such as the FD of individual cells. To overcome these issues, MacDonald et al. , utilized an autocorrelation technique capable of analyzing patterns within individual cells. However, this autocorrelation technique is particularly time consuming and not suitable for mapping FD in large image sets.
The goal of our study was to develop a rapid and automated method for quantifying fractal patterns within images on a pixel-by-pixel basis. We modified a previously developed technique known as the Blanket method [11,19,20]. The Blanket method interprets images as a three-dimensional surface topography, where the intensity at each pixel corresponds to the height along the z axis. The surface area (SA) of the topographic map is computed as a function of scale (pixel size) [11,21], and the power law relationship between SA and pixel size is used to determine FD. Caldwell et al.  and Byng et al.  employed this approach for the fractal analysis of mammographic parenchymal morphologies, which are strongly related to the development of breast cancer. Chappard et al.  also used the Blanket method to characterize the texture of trabecular bone in X-Rays. In the current study, we propose a Modified Blanket Method (MBM), which utilizes convolution to rapidly compute local measurements of SA, and thus produce maps of FD within images. By acquiring local FD values at each pixel, the analysis of specific structures or individual cells is possible. Here, we present a series of studies in which we analyze simulated cell images containing fractal patterns of known FD to show that the MBM is sensitive to a wide range of FD values. We evaluate the effects of different convolution kernels on the accuracy of the FD measurements and compare them to PSD analysis. Moreover, we demonstrate the biomedical application of the MBM on NADH autofluorescence images obtained from parental and radiation resistant A549 lung cancer cells.
2.1 Modified Blanket Method
The MBM was developed as an automated technique to convolve gradient images at different resolutions and compute FD. The image analysis algorithm was coded and executed in MATLAB. The length scale over which the local fractal dimension is analyzed is based on a convolution kernel defined by a binary disk with a user-defined radius. The algorithm begins by resampling both the image of interest and convolution kernel (Fig. 1(a)) proportionally. The number times the image is resampled at different image resolutions is determined by the size of the kernel. Once the image has been resampled, the horizontal (X) and vertical (Y) gradients of the image are computed using X- and Y- gradient kernels [1,-1,0] and [0,-1,1]T respectively (Fig. 1(b)). The absolute value of the gradient in the X- and Y- directions are then summed by convolving with a binary disk kernel (Fig. 1(c)). The summed horizontal and vertical gradients are added together to compute a local surface area (SA) map (Fig. 1(d)). This SA map is subsequently resampled back to the original dimensions of the image and stored (Fig. 1(e)). These steps (Fig. 1(a-e) are repeated as the convolution kernel is reduced in size by increments of 1 pixel until a 1x1 kernel remains, ultimately producing a set of SA maps derived from different image resolutions (Fig. 1(f)). The power law exponent (β) relating local surface area measurements and pixel sizes at each pixel location can be derived from the linear change in log-transformed values between SA maps (Fig. 1(g)). Fractal dimension at each pixel location (Fig. 1(g)) can then be determined through the equation:Fig. 1(h)) can be generated by assigning the corresponding FD value to each individual pixel.
2.2 Validation: simulated cell images
To assess the accuracy of the MBM to approximate fractal values, we created simulated cells containing intensity clouds of known fractal dimension as in previous studies . The fractal clouds were obtained by filtering white noise in the frequency domain and transforming it back to the spatial domain . These fractal clouds were overlaid using masks of circular cells with varying diameter and cellular nuclei. We evaluated these simulated cell images over a range of fractal dimension values (FD = 2.0 to FD = 4.5) and cell diameters (d = 80, d = 120, and d = 160 pixels) commonly observed in biomedical images. FD values were obtained from five repetitions of each combination of FD and cell diameter. The fractal clouds and cell positions were randomized during each repetition of the simulation, and FD was derived using our MBM method. This process was then repeated for three different cell densities (10, 20, and 30% of the total image area) to measure differences caused by increasing cell quantities. Mean error and standard deviation was computed using the measured FD values obtained from the simulation repetitions.
In addition to the varying fractal clouds and cell diameters, we analyzed the simulated cell images over a range of different convolution radii (CR = 5, CR = 10, CR = 15, and CR = 20 pixels) to obtain local fractal dimension values. Errors in FD measurements associated with abrupt intensity changes at cell borders have been previously described  and at larger CR values it is possible to visualize these edge effects. In order to estimate the differences in the average FD caused by such artifacts, we examined the effect of cell mask erosion to isolate the regions of interest (i.e. the center of simulated cells) from edge effects. We ran simulations with different combinations of convolution radii (CR) and erosion radii (ER) and computed the average FD (Fig. 2).
2.3 Comparison to the power spectral density approach
To compare the results obtained with the MBM, we also employed PSD analysis as in previous work [4,8,10,18]. The PSD approach computes a radial average of the 2D Fourier transform of the image, resulting in a 1D PSD curve . To reduce abnormalities caused by regions that do not contain intercellular features, a clone stamping pre-processing step was used, which fills the background with copies of the intracellular image features . The 1D PSD curve was fit to a specific range of frequency values, corresponding to the most linear portion of the plot as in previous studies [4,8]. The power law exponent (β) commonly reported in other studies [1,8] can be related to fractal dimension using Eq. (3).
The accuracy of both methods (MBM and PSD) over a wide range of FD values (FD = 2.0 to FD = 4.5) was evaluated using the simulated cell images. The computational time of each method was also measured using five repetitions at different image resolutions (128 x 128, 256 x 256, 512 x 512, 1024 x 1024, 2048 x 2048, and 4096 x 4096 pixels). For the MBM, we evaluated three different CR values (CR = 5, CR = 10, and CR = 25) to assess their effect on the speed of the fractal analysis. The computational times recorded at each repetition were then used to compute mean times and their standard deviation.
2.4 The biomedical application to understanding cancer metabolism
To demonstrate the sensitivity of this MBM technique using experimental data, we compared FD values obtained from NADH autofluorescence images of cell mitochondria from parental and radiation-resistant A549 lung cancer cells. In order to generate radiation resistant A549 cells, the parental cell line was subjected to multiple doses of 2 Gy radiation . Images from both cell lines were then acquired at different post-radiation time points (t = 0, t = 1, t = 12, and t = 24 hours), with a total of nine images per time period. Two-photon excited fluorescence from NADH was collected using a 20x water-immersion objective (NA = 0.75) with a 460/40 nm bandpass filter and a Ti:Sapphire laser source tuned to 755 nm. Image resolution was set to 512 x 512 pixels (130 µm x 130 µm) with 16-bit depth. FD values were obtained using both the MBM and PSD method.
2.5 Statistical analysis
An ANOVA design with image fields nested within cell culture dishes was used to identify differences in FD between cell lines (parental and radiation-resistant) and time points. Significance was based on an α = 0.05.
3.1 Method validation
Utilizing simulated microscopy images of cells with intracellular fractal clouds of known FD, we demonstrated that the MBM can compute FD values that successfully estimate true FDs ranging from 2.00 to 4.50 (Fig. 2(a), (b)). For 120 pixel diameter simulated cells, a CR of 5 produced a mean error of 1.4 ± 1.0%. However, increasing the convolution radius (CR) over which the FD is computed slightly decreased accuracy (Fig. 2(a)), with mean errors of 2.1 ± 1.9%, 3.6 ± 3.3% and 5.0 ± 3.7% for radii of 10, 15 and 20, respectively. The FD values obtained from the MBM were insensitive to differences in cell size (d = 80, d = 120, and d = 160 pixels), cell density (10, 20, and 30% cell coverage per image), or overall image intensity. Using cell images with FD = 2.5, the MBM analysis resulted in a mean error of 0.72 ± 0.49%, 1.08 ± 0.34%, and 1.76 ± 0.76% for cell diameters of 80, 120, and 160 pixels respectively. A similar analysis was performed with varied cell densities, resulting in a mean error of 1.97 ± 0.88%, 2.07 ± 0.78%, and 2.14 ± 0.54% for 10, 20, and 30% cell coverage respectively.
Mean error in average FD measurements from simulated cell images also significantly increased if the cell borders were not removed from analysis through the morphological operation of erosion (Fig. 2(c)). When using a CR of 5, the lack of an erosion step prior to FD averaging resulted in a mean error of 3.6 ± 2.6%. This error was even greater for CR of 10, which increased error from 2.1 ± 1.9% to 7.1 ± 4.5% (Fig. 2(d)). Pixelwise FD maps indicate that inaccurate measurements of FD approaching a value of 2 are observed along the cell-background border (Fig. 2(e)) in simulated images. An erosion radius equivalent to the convolution radius is necessary to remove these artifacts produced when abrupt, non-fractal changes in image intensity are present.
Unlike the MBM, PSD analysis was not accurate at FDs below 2.75 (Fig. 3(a)). While the MBM had a mean error of 1.4 ± 1.0% and maximum error of 4.26% over FDs from 2.0 to 4.5, PSD yielded a mean error of 6.0 ± 10.9% and maximum error of 35.50% (Fig. 3(a), (b)). The MBM was also able to compute FD with significantly shorter computational times than the PSD approach, particularly at image resolutions of 1024x1024 and higher (Fig. 3). MBM analysis with CR = 5 had computational times on a standard desktop computer that were up to 42-times faster on 4096x4096 images and 3-times faster on 128x128 images than PSD analysis. Increasing the CR resulted in longer computational times, but CR = 10 was still 17-times faster than PSD analysis of 4096x4096 images and 1.5-times faster on 128x128 images (Fig. 4). For both the PSD and MBM, image resolution and computational time scaled according to a power law.
3.2 Lung cancer cells
The MBM revealed a significant change in the mitochondrial FD of radiation-resistant lung cancer cells between 12 and 24 hours after radiation (p<0.0001) (Fig. 5). Previous work demonstrated that the optical redox ratio of FAD/(NADH + FAD) decreases in the radiation-resistant cell line at 24 hours post-radiation (p = 0.0141), indicating a shift towards glycolytic metabolism , which matches the changes in FD observed here with the MBM. Interestingly, erosion did not produce significant differences in the FD values, despite providing a noticeable improvement in the accuracy of simulated images. While the MBM results matched biochemical changes in cell metabolism observed over time, a change in FD between post-radiation time points was not detectable through the conventional PSD analysis (p = 0.9991). Unlike the MBM and previous redox ratio measurements , PSD analysis indicated significant differences between parental (control) and radiation-resistant cancer cell lines at every time point (Fig. 5(c)).
Our study demonstrates that the MBM is capable of rapidly and accurately quantifying FD within individual cells. By computing local changes in intensity gradients at different image resolutions, the MBM can obtain fractal dimension values from the power law exponent relating surface area and pixel size . Cell simulations demonstrated that the MBM yields accurate FD measurements over a wide range of FD values (2.00-4.50). In addition, the improved accuracy of the MPM relative to PSD analysis, allowed for the detection of structural changes in parental and radiation-resistant A549 cell mitochondria over time, which was supported by previous optical redox ratio measurements that suggested a shift in cellular metabolism [23,24] (Fig. 5(c)).
Compared to PSD analysis, the MBM offers accuracy over a greater range of FD in biomedical images. Cell simulations using the PSD method resulted in inaccurate FD approximations for FD<2.75 (Fig. 3(a)) with a high maximum error of 19.9% (Fig. 3(b)). This reduced accuracy for the PSD approach is caused by the susceptibility of the method to scale-dependent image features (i.e. cell boundaries). We employed a technique proposed by Xylas et. al.  to minimize the effect of scale-dependent features by replicating cellular patterns in the background of the image through clone stamping. However, Xylas et al. and the current study demonstrate that this pre-processing step cannot recover FD values below 2.5 when using PSD analysis. At lower FD values, the image intensity patterns vary over lengths that can approach or exceed the size of the cell, resulting in cell boundaries visible in the clone-stamped images (Fig. 3(d)). These cell boundaries were also visible in experimentally-derived clone-stamped images (Fig. 5(b)) as well, suggesting possibly inaccurate measures of fractal dimension through the PSD. The MBM does not demonstrate this same sensitivity to cell boundaries in experimentally-derived images (Fig. 5(a)), regardless of whether pixels at the edge of the cell are removed through erosion prior to calculating average values.
The effects of cell or object boundaries can be still observed in the MBM-derived FD maps of the simulated images. However, compared to the PSD, the MBM’s edge-effects are limited to the relatively small number of pixels at the boundary of cell objects (Fig. 2(e)). Increasing the size on the convolution disk can increase the number of pixels affected by these large (non-fractal) changes in intensity (e.g. cell or nuclei borders). To overcome the consequences of these artifacts, erosion of the cell mask using the same kernel size as that used during convolution yielded more accurate average FD values in simulated images (Fig. 2(e)). The erosion procedure thins the cell mask by a specific number of pixels and enables computation of an average FD value from only unaffected pixels within the cell or region of interest. In general, the simulations revealed that using the MBM with equal CR and ER size yields a better FD approximation. However, erosion did not have a large impact in the average FD value when it was performed on the A549 cell images (Fig. 5). This might be related to the more gradual intensity change between cell and background regions within actual biomedical images. It is important to highlight that although cell boundaries can influence the MBM in some cases, they are not as sensitive to the final FD measurement as they are when using the PSD approach.
The MBM offers improvements in accuracy and speed relative to PSD analysis of mitochondrial clustering, but the primary advantage of the MBM is the ability to report FD on a pixel-by-pixel basis. Given that single cell assessments of mitochondrial function are possible using time-resolved and intensity-based NADH and FAD autofluorescence imaging [14,25], the ability to also assess rapidly mitochondrial structure on a cell-by-cell basis can enable complementary information to assess disease and repair . Fractal analysis has been employed in a number of biomedical applications beyond mitochondrial clustering, including the characterization of electroencephalogram (EEG) signals , the morphology of the brain in MRI images , and the anisotropy degree of bone structures . The MBM may be particularly useful in detecting changes in vasculature surrounding tumors produced by aberrant angiogenesis . Changes in extracellular matrix organization in different disease or repair states may also offer new applications for this method capable of rapidly measuring regional differences in structural organization. In summary, the modified blanket method for computing fractal dimension offers a unique ability to rapidly map out changes in local image patterns for a variety of biomedical applications beyond mitochondrial organization.
National Institutes of Health (NIH R00EB017723, NIH R01AG056560); Arkansas Biosciences Institute; NIGMS P20GM103625; and the Jean Ostermeier Memorial Cancer Research Endowed Award.
NIH R00EB017723 (KPQ), NIH R01AG056560 (KPQ), Arkansas Biosciences Institute (KPQ, NR), NIGMS P20GM103625 (RPMD), and the Jean Ostermeier Memorial Cancer Research Endowed Award (IV).
The authors declare that there are no conflicts of interest related to this article.
1. J. Li, Q. Du, and C. Sun, “An improved box-counting method for image fractal dimension estimation,” Pattern Recognit. 42(11), 2460–2469 (2009). [CrossRef]
3. J. H. Brown, V. K. Gupta, B. L. Li, B. T. Milne, C. Restrepo, and G. B. West, “The fractal nature of nature: power laws, ecological complexity and biodiversity,” Philos. Trans. R. Soc. Lond. B Biol. Sci. 357(1421), 619–626 (2002). [CrossRef] [PubMed]
4. D. Pouli, M. Balu, C. A. Alonzo, Z. Liu, K. P. Quinn, F. Rius-Diaz, R. M. Harris, K. M. Kelly, B. J. Tromberg, and I. Georgakoudi, “Imaging mitochondrial dynamics in human skin reveals depth-dependent hypoxia and malignant potential for diagnosis,” Sci. Transl. Med. 8(367), 367ra169 (2016). [CrossRef] [PubMed]
5. A. A. Dukes, V. S. Van Laar, M. Cascio, and T. G. Hastings, “Changes in endoplasmic reticulum stress proteins and aldolase A in cells exposed to dopamine,” J. Neurochem. 106(1), 333–346 (2008). [CrossRef] [PubMed]
6. K. P. Quinn, G. V. Sridharan, R. S. Hayden, D. L. Kaplan, K. Lee, and I. Georgakoudi, “Quantitative metabolic imaging using endogenous fluorescence to detect stem cell differentiation,” Sci. Rep. 3(3432), 3432 (2013). [CrossRef] [PubMed]
9. S. B. Berman and T. G. Hastings, “Dopamine oxidation alters mitochondrial respiration and induces permeability transition in brain mitochondria: implications for Parkinson’s disease,” J. Neurochem. 73(3), 1127–1137 (1999). [CrossRef] [PubMed]
10. J. Xylas, A. Varone, K. P. Quinn, D. Pouli, M. E. McLaughlin-Drubin, H. T. Thieu, M. L. Garcia-Moliner, M. House, M. Hunter, K. Munger, and I. Georgakoudi, “Noninvasive assessment of mitochondrial organization in three-dimensional tissues reveals changes associated with cancer development,” Int. J. Cancer 136(2), 322–332 (2015). [CrossRef] [PubMed]
11. C. B. Caldwell, S. J. Stapleton, D. W. Holdsworth, R. A. Jong, W. J. Weiser, G. Cooke, and M. J. Yaffe, “Characterisation of mammographic parenchymal pattern by fractal dimension,” Phys. Med. Biol. 35(2), 235–247 (1990). [CrossRef] [PubMed]
12. T. M. Cabral and R. M. Rangayyan, Fractal Analysis of Breast Masses in Mammograms (Morgan & Claypool Publishers, 2012).
15. D. Pendin, R. Filadi, and P. Pizzo, “The Concerted Action of Mitochondrial Dynamics and Positioning: New Characters in Cancer Onset and Progression,” Front. Oncol. 7(102), 102 (2017). [CrossRef] [PubMed]
16. F. Bartolomé and A. Y. Abramov, “Measurement of mitochondrial NADH and FAD Autofluorescence in Live Cells,” in Mitochondrial Medicine Methods in Molecular Biology, V. Weissig and M. Edeas, eds. (Humana Press, 2015).
17. K. Foroutan-pour, P. Dutilleul, and D. L. Smith, “Advances in the implementation of the box-counting method of fractal dimension estimation,” Appl. Math. Comput. 105(2), 195–210 (1999). [CrossRef]
18. Z. Liu, D. Pouli, C. A. Alonzo, A. Varone, S. Karaliota, K. P. Quinn, K. Münger, K. P. Karalis, and I. Georgakoudi, “Mapping metabolic changes by noninvasive, multiparametric, high-resolution imaging using endogenous contrast,” Sci. Adv. 4(3), eap9302 (2018). [CrossRef] [PubMed]
20. M. P. Paskaš, I. S. Reljin, and B. D. Reljin, “Multifractal Framework Based on Blanket Method,” The Scientific World Journal 2015, 894546 (2015)
22. D. Chappard, A. Chennebault, M. Moreau, E. Legrand, M. Audran, and M. F. Basle, “Texture analysis of X-ray radiographs is a more reliable descriptor of bone loss than mineral content in a rat model of localized disuse induced by the Clostridium botulinum toxin,” Bone 28(1), 72–79 (2001). [CrossRef] [PubMed]
23. K. Alhallak, S. V. Jenkins, D. E. Lee, N. P. Greene, K. P. Quinn, R. J. Griffin, R. P. M. Dings, and N. Rajaram, “Optical imaging of radiation-induced metabolic changes in radiation-sensitive and resistant cancer cells,” J. Biomed. Opt. 22(6), 060502 (2017). [CrossRef] [PubMed]
24. D. E. Lee, K. Alhallak, S. V. Jenkins, I. Vargas, N. P. Greene, K. P. Quinn, R. J. Griffin, R. P. M. Dings, and N. Rajaram, “A Radiosensitizing Inhibitor of HIF-1 alters the Optical Redox State of Human Lung Cancer Cells In Vitro,” Sci. Rep. 8(1), 8815 (2018). [CrossRef] [PubMed]
26. W. Y. Hsu, C. C. Lin, M. S. Ju, and Y. N. Sun, “Wavelet-based fractal features with active segment selection: application to single-trial EEG data,” J. Neurosci. Methods 163(1), 145–160 (2007). [CrossRef] [PubMed]
28. R. Jennane, R. Harba, G. Lemineur, S. Bretteil, A. Estrade, and C. L. Benhamou, “Estimation of the 3D self-similarity parameter of trabecular bone from its 2D projection,” Med. Image Anal. 11(1), 91–98 (2007). [CrossRef] [PubMed]
29. C. E. Priebe, J. L. Solka, R. A. Lorey, G. W. Rogers, W. L. Poston, M. Kallergi, W. Qian, L. P. Clarke, and R. A. Clark, “The application of fractal analysis to mammographic tissue classification,” Cancer Lett. 77(2-3), 183–189 (1994). [CrossRef] [PubMed]