Abstract
Fractal geometry has been applied widely in the analysis of medical images to characterize the irregular complex tissue structures that do not lend themselves to straightforward analysis with traditional Euclidean geometry. In this study, we treat the nonfractal behaviour of medical images over large-scale ranges by considering their box-counting fractal dimension as a scale-dependent parameter rather than a single number. We describe this approach in the context of the more generalized Rényi entropy, in which we can also compute the information and correlation dimensions of images. In addition, we describe and validate a computational improvement to box-counting fractal analysis. This improvement is based on integral images, which allows the speedup of any box-counting or similar fractal analysis algorithm, including estimation of scale-dependent dimensions. Finally, we applied our technique to images of invasive breast cancer tissue from 157 patients to show a relationship between the fractal analysis of these images over certain scale ranges and pathologic tumour grade (a standard prognosticator for breast cancer). Our approach is general and can be applied to any medical imaging application in which the complexity of pathological image structures may have clinical value.
1. Introduction
Many biological phenomena exhibit chaotic or fractal-like behavior, and these features have been studied extensively [1–5]. One area of widespread use has been the application of fractal geometry in the analysis of medical images as it lends itself naturally to the pragmatic characterization of irregular non-Euclidean structures found in medical imaging [6–9]. However, care must be taken when applying fractal analysis to natural objects [10]. The quintessential requirement for an object to be a fractal is for the object to exhibit a form of self-similarity to arbitrarily small scales. As such, actual fractals do not exist in nature, since there is a fundamental natural limitation to the scaling behaviour of natural objects when their substructures approach the atomic scale. Even real renderings of mathematical fractals cannot be truly fractal because of the finite resolution of the rendering. In general, the scaling behaviour of natural objects depends on the scale at which the objects are considered. Despite this limitation, most fractal analysis techniques have focused on characterizing the fractal behaviour of natural objects by finding an interval of scales in which these objects have an approximately constant scaling behaviour [3].
Some authors have noted the danger in making this interpretation because of the possibility of illusory “fractal” behaviour where none is actually present [11], as well as the difficulty in defining the interval of scales over which the object has a consistent scaling behaviour [12]. To limit these difficulties, our approach in recent studies was to manually pick the interval of scales based on physical considerations (i.e., range of the sizes of the histological structures of interest) and the linearity of the dependence of entropy on scale [13–16].
Another approach used by others to resolve this difficulty has been to model the scale dependence of the fractal dimension in a rendering of a mathematical fractal by fitting a functional model to the dimension as a function of scale and interpreting the fit parameters as meaningful fractal dimensions for the renderings [17, 18]. The conceptual arguments behind these models are derived for renderings of mathematical fractals, and they may not necessarily hold for images of natural objects. However, the authors of these studies are aware that natural objects do not fit the mold of real mathematical fractals.
In this study, we adopt the view expressed by others that the fractal dimension of a natural object changes with scale. Assuming that physically meaningful information in an image is contained in the scaling behaviour of the image in a certain interval of spatial dimension, then a thorough and objective approach to extract this information would be to assess the clinical relevance of different scale ranges and determine an absolute scale at which the fractal dimension is relevant for a given medical imaging application. That is, the spatial scale interval chosen for the analysis of medical images would be selected on the basis of patient classification performance and patient outcome for diagnostic and prognostic applications, respectively, rather than the resemblance of the image structures to a fractal at these image scales.
In order to address the shortcomings of existing techniques, we propose an approach to analyzing the scaling behaviour as a function of scale itself, providing an estimate of the fractal dimension as a scale-dependent parameter rather than a single fixed value. We present this fractal analysis approach in the context of the more generalized Rényi entropy [19], in which we can also compute the information and correlation dimensions of images [20].
In addition, we describe a novel fractal analysis algorithm based on integral images [21] which speeds up the computation of these generalized dimensions by orders of magnitude, and we validate our algorithm and approach by applying it to increasingly complex data, producing meaningful results throughout. We also illustrate our method with the analysis of a set of histology images acquired from tissue samples of breast tumours with modified Bloom-Richardson grades of 1, 2, and 3, in which we determine the spatial scale intervals that exhibit statistically significant differences in fractal dimensions between the different tumour grades. Looking forward, this approach provides a rapid tool to determine an absolute scale at which the generalized fractal dimensions are relevant and may allow for objective interimage comparison of medical images acquired at different resolutions.
2. Fractal Dimension
Fractal dimension (FD) is a generalization of the intuitive notion of topological dimension and allows for noninteger dimensions; for example, a point, a line, and a plane have both topological and fractal dimensions of 0, 1, and 2, respectively. However, fractal objects such as the Koch snowflake and Sierpinski carpet both have a topological dimension of 1 but have fractal dimensions of 1.26 and 1.89, respectively. The fractal dimension of an object is fundamentally a reflection of its scaling behaviour. Hence, to provide an estimate of the fractal dimension of a real image, we first need to define a scale-dependent function , with pixel intensities defined for , where is the area contained in the image and is the analysis scale. For mathematical fractals, the fractal dimension is usually defined as However, for a real image it is meaningless to analyze an infinitely small scale by letting , and so at any given scale we can define the fractal dimension as The function corresponding to the box-counting dimension is commonly used in medical imaging because of its conceptual and algorithmic simplicity. In this method, the parameter is the side length of the square boxes into which the image is partitioned, and is the logarithm of the number of squares of this size in the image that contains at least one pixel that is a part of a structure of interest.
2.1. Generalized Dimensions
The box-counting function is a special case of a more generalized class of functions, the Rényi entropies [19], which also generalize Shannon entropy. The Rényi entropy of order of a probability distribution is given by The values are the set of natural measures of the probability distribution. For an image data set and a scale , the values are found by dividing the image into squares, and for each square by finding the proportion of the total image intensity that is contained in the th square: where is the domain of the image contained in the th square. Definition (4) implies because the set is entirely covered by the subsquares .
The fractal dimensions obtained by using different orders generate the multifractal spectrum [6] of an image. There are three values of that have clear physical significance [22]. When , we treat the term in (3) as the limit of as . Hence, for , , while for , . The entropy is therefore the logarithm of the number of nonzero values of in the image, which is the same as logarithm of the number of nonempty squares of size in the image, equaling the box-counting function from before. The Rényi entropy of order 0 will thus yield the box-counting dimension of the image.
If , we again treat (3) as a limit when . In this case, by using L’Hôpital’s rule for the limit as , the Rényi entropy reduces to the Shannon entropy as follows: The dimension obtained from the Shannon entropy is known as the information dimension of the image. This measure of the fractal dimension gives the rate at which information is gained about the structure of the image as the resolution of the image increases.
The third value of with a clear physical meaning is , which gives rise to , the correlation dimension of the image, addressing the number of neighbours a point of the structure has as a function of scale. That is, gives the power law that relates the number of other image pixels that are within a range of a given pixel to the value of .
In essence, larger values assign a greater weight to the brighter parts of the image being analyzed. This is particularly useful for the analysis of medical images in which both the spatial structure and relative intensity of edge structures may carry useful information about the image.
3. Methods and Materials
3.1. Scale Dependence of Fractal Dimension
As mentioned in the introduction, natural objects do not exhibit true fractal behaviour, as their scaling behaviour depends on the scale at which the objects are considered. Hence, we propose an objective approach to determine the scale or scale interval in which the different orders of fractal dimensions may have physical relevance (e.g., diagnostic or prognostic value in medical imaging).
To estimate the fractal dimension as a function of scale, we use definition (2) and measure the entropies given by (3). To differentiate the entropy with respect to the analysis scale, we use the locally weighted regression and smoothing scatterplots (LOWESS) method [23], which is widely used in situations where a good theoretical model for the observed data does not exist. A low-order polynomial is fitted to a weighted subset of the data around each data point, and then all parameters (such as the derivatives) of the fitted curve can be extracted from this polynomial. In our case, we are interested in the first derivative of the entropy with respect to scale, which immediately gives a scale-dependent estimate of the fractal dimension of the image. The advantage of this method over direct numerical differentiation is its vastly enhanced robustness to noise, while its advantage over a fit to the entropy or numerical derivative data, as used by others [17, 18], is in its greater flexibility and scale resolution, as well as a lack of assumptions that may not be valid for images of natural objects.
3.2. LOWESS Method
To estimate parameters of a scatterplot of a data set given as a function of using the LOWESS method [23], for each , we fit a polynomial to the data in such a way as to minimize the weighted sum of the squared residuals given by , where is a weighting function. The fitted curve is then approximated around by the polynomial . In particular, for , we estimate the th derivatives of and as equal for .
The function is central to the LOWESS method, since it makes the regression locally weighted. The simplest weighting function that can be used makes the fit around a point local a Gaussian with a standard deviation : Larger values of will smooth out the fit more than smaller values, producing a fit that is more resistant to noise at the cost of resolution in . In the case of scale-dependent fractal analysis, the input data and to the LOWESS model are and , from which we can obtain .
3.3. Rapid Computation of the Natural Measures
The most common operation encountered in the box-counting approach to fractal image analysis is the determination of the natural measure given in (4) of a certain rectangular (usually square) subset of the image. If the subset is bounded by and , while the entire image is bounded by and , finding the natural measure is equivalent to calculating
For example, in Figure 1, the sum of the image intensity in image over the dotted (red) rectangle bounded by and is equal to .
Direct summation of the image intensity over a rectangle has computational cost proportional to the area of the rectangle. A method to speed up the summation is to first compute the following integral image (also known as the summed area table) [21]: This computation reduces the subsequent summation in (7) to a simple arithmetic operation (e.g., see Figure 1): where we define and for all and .
The usefulness of this algorithm lies in its ability to speed up natural measure computation. The computational complexity of a naïve algorithm, in which for each scale we determine the natural measures, (4), by directly summing the image intensity is of order , where is the area of the image and is the number of different values of considered. However, by using the integral image algorithm approach, the computational complexity falls to since the cost of computing the natural measures at a scale is approximately equal to and is small. The integral image approach is particularly advantageous for large values of where a large number of summations are replaced with only a few subtractions. For a typical fractal scale analysis, the integral image approach is 10 to 100 times faster than the naïve algorithm.
3.4. Algorithm
Our method was implemented in MATLAB Version 7 (The MathWorks, Inc., Natick, MA, USA). The steps of the method are as follows.(i)Select the analysis scales and entropy order values. The choices made are based on purely physical considerations, for example, the range of sizes of the structures of interest in the image being analyzed. To be general, one can begin by analyzing a wide range of and values and select the physically meaningful subset of scales and orders. (ii)Normalize the original image by dividing every element by the sum of the pixel intensities of the entire image, giving the natural measure contained in each image pixel. (iii)Calculate the integral matrix for the normalized image . For each scale , rescale the image using the integral matrix. That is, calculate the matrix (iv) For each value of , raise the natural measure matrix to the power except for obtaining ; for , we calculate the matrix . (v)Use (3) and (5) to give for all , while for , .
Our method was tested on the following data sets: exact renderings of deterministic strictly self-similar fractals, randomized renderings of deterministic fractals, randomized renderings of statistically self-similar fractals, and a set of breast histological tissue samples.
3.5. Breast Cancer Tissue Specimens
In a previous study [16], we had retrospectively selected 408 patients with primary invasive ductal carcinoma (IDC) of the breast from Calgary regional hospitals after appropriate ethics approval from the institutional review board. The breast tissue from these patients was used to construct tissue microarray (TMA) cores (each with a diameter) stained with pan-cytokeratin to highlight the morphology of epithelial architecture. The number of cores per patient ranged from one to three. Images of the cores were acquired with an effective magnification of 6.3x using an AxioCam HR digital camera (Carl Zeiss, Inc.) mounted on an optical microscope (Zeiss Axioscope). The images were saved at the camera’s native resolution of 1300 × 1030 pixels in tagged image file format (tiff). In this previous study [16], we found the box-counting fractal dimension of the breast cancer TMA core images to be an independent and statistically significant prognosticator. However, the study did not include an explicit examination of the role of the scale range on the fractal dimensions computed from the images. Instead, fractal dimension was computed from plots of the slope of versus over a scale range of (), which was chosen based on a visual assessment of the range of the linear region of a small random sample of plots taken from the whole image set.
In this study, we selected all the cases from our previous study set of 408 patients [16] that satisfied the conditions of having exactly three evaluable TMA core images and contained pathologic grade information. This selection resulted in a set consisting of a total of 157 patients in the following tumour grade categories: 56 grade 1, 84 grade 2, and 17 grade 3 tumours. We applied our method to these cases, and the analysis was used to demonstrate the capabilities of the algorithm and to check that our previous choice for the scale range was a judicious one. The grayscale images of the tissues were converted into black-and-white outline images by thresholding. An example of the overall analysis process for a breast cancer tissue core is shown in Figure 2. To determine the optimal thresholding level for the edge detection, we varied the threshold level for each image to maximize the fractal dimension in the scale interval.
(a)
(b)
(c)
(d)
(e)
The motivation for using three cores per patient (chosen from different tumour regions) was to ensure that we had a representative sample of a heterogeneous tumour. For the final analysis, we selected the one core from each patient that had the greatest average fractal dimension in the scale range. The rationale for choosing the core with the greatest fractal dimension is that it is likely representative of the portion of a possibly heterogeneous tumour that has deviated most from normal epithelial breast morphology, and therefore it is the most probable indicator of abnormal and/or aggressive tumour growth with metastatic potential. For the analysis, we produced a single curve of fractal dimension as a function of scale for each patient. We averaged the fractal dimensions within each tumour grade category.
3.6. Statistics
We performed the statistical analysis using the MATLAB Statistics Toolbox 7.4 (The MathWorks, Inc., Natick, MA, USA). We quantified the differences between the three tumour grade categories using the nonparametric Kruskal-Wallis analysis and a multiple comparison test (MATLAB functions kruskalwallis and multcompare, resp.).
4. Results
4.1. Assessment of Scale-Dependent Algorithm on Renderings of Strictly Self-Similar Fractals
Figure 3 shows the results of our approach applied to renderings of strictly self-similar mathematical fractals. For the Koch snowflake, Pascal triangle, and Sierpinski carpet, the estimated fractal dimensions are consistently within 0.1 of the Hausdorff dimension of the mathematical fractal being rendered throughout the scale range of 10 to 250 pixels. Note that both the mean and maximum deviations of the measured fractal dimension from the Hausdorff dimension decrease with increasing fractal dimension, increasing value of in (6), and have a minimum in the range of 20 to 50 pixels, where neither the small nor large-scale granularities of the image affect the estimate of the fractal dimension. In this smaller interval of scales, the measured fractal dimensions have a root-mean-squared deviation of 0.074, 0.040, 0.086, and 0.029 for and 0.045, 0.018, 0.052, and 0.034 for from the Hausdorff dimension, respectively, for the four fractals presented.
(a) 998 × 1152 pixel rendering of the boundary of a Koch snowflake
(b) Scale dependence of entropy of (a)
(c) Scale dependence of fractal dimension of (a)
(d) 2000 × 2000 pixel rendering of the boundary of a Koch island
(e) Scale dependence of entropy of (d)
(f) Scale dependence of fractal dimension of (d)
(g) 1460 × 1460 pixel rendering of Pascal’s triangle (mod 3)
(h) Scale dependence of entropy of (g)
(i) Scale dependence of fractal dimension of (g)
(j) 2187 × 2187 pixel rendering of the Sierpinski carpet
(k) Scale dependence of entropy of (j)
(l) Scale dependence of fractal dimension of (j).
For the Koch island boundary (Figure 3(d)), a jump of fractal dimension from 1 to 1.5 (the Hausdorff dimension) can be seen at pixels. This behaviour is consistent with the real structure of the rendering, which at small scales consists of straight line segments of 4–8 pixels long. Hence, below a scale of 4 pixels, the rendering is really linear and hence has a fractal dimension of 1 at this scale.
In Figure 3(i), we compare the estimate of fractal dimensions produced by using the values of 0.3 and 0.5 for in (6). The larger value of produces a smoother dependence of fractal dimension on scale, which is close to the Hausdorff dimension of Pascal’s triangle. On the other hand, the smaller value gives an estimate which is more locally accurate, showing the nonuniform scaling behaviour of the rendering, which can be seen in the oscillations in the fractal dimension as a function of scale. In fact, this inhomogeneous scaling behaviour is seen in all four of the sample renderings. These oscillations occur because the real renderings of the fractals have discrete characteristic scales. For example, the Pascal’s triangle mod 3 rendering (Figure 3(g)) consists of black and white triangles of several discrete scales: 3, 9, 27, 81, 243, and 729 pixels. Around each of these scale values, the fractal dimension experiences a large drop. This drop occurs because when the analysis scale grows through each of these scales, the white triangles in the rendering become “invisible” to the box-counting algorithm, causing a smaller than expected drop in the entropy and consequently causing a dip in the box-counting dimension. However, between these characteristic scales, the image becomes nearly 2-dimensional, just like the plane in which the image is contained, because the analysis cannot “see” any change in the image features.
4.2. Assessment of Scale-Dependent Algorithm on Renderings of Statistical Fractals
We further tested our algorithm on two kinds of statistical fractals, which are generated by random processes, but nonetheless possess a statistical form of self-similarity.
4.2.1. Randomized Sierpinski Triangle
An approximation to a Sierpinski triangle can be generated by an iterative random process known as the “chaos game” [24]. In each step of the game, one new point is added to the rendering of the triangle. We employed this method with a variable number of iterations (from to ) to generate several Sierpinski triangle approximations, with an example for 10000 points shown in Figure 4(a). When the structure is rendered with only a few points, the structure of the randomized triangle is essentially point-like or 0-dimensional at small scales. For example, with only 1000 points in a 2000 × 2000 pixel image, the expected mean distance between points is pixels; indeed, the fractal dimension measured for begins to grow around that scale, but already has the Hausdorff dimension scaling behaviour at larger scales above 100 pixels.
(a)
(b)
In this case, we can see that our algorithm correctly identifies the scaling behaviour of the fractal renderings: the large-scale behaviour of the fractal quickly reaches the correct scaling power law, while the smaller-scale features are really 0-dimensional. However, as more and more points are added to the approximation, the fine-scale structure of the fractal is filled in, and the fractal dimension approaches its ideal value at all scales.
4.2.2. Brownian Surfaces
The fractal dimension of nine 1025 × 1025 pixel renderings of Brownian surfaces with dimensions 2.1 through 2.9 were evaluated by finding cross-sections through the images (which should have a dimension that is one dimension less than the surfaces themselves). The results of this analysis are shown in Figure 5. Several features can immediately be seen from Figure 5(c). All of the outlines are inherently 1-dimensional at small scales because of their linear structure. The Brownian surface outlines (i.e., Figure 5(b)) pack the image more densely as the scale of the image increases, causing the fractal dimensions to approach 2 for all the cross-sections. These features are not due to a flaw in the algorithm, but rather reflect the true behaviour of the curves obtained by slicing through the Brownian surfaces. The behaviour of other measures of fractal dimension, such as other components of the multifractal spectrum, is similar to the box-counting dimension shown in the figure. The fractal dimensions match the expected values most closely in the 10–40 pixel range, with a root-mean-square error of less than 0.06 in this interval. This example illustrates the great sensitivity of fractal dimension measures to the scale at which they are computed and the consequent need for a scale-dependent measure of the fractal dimension to quantitatively estimate this sensitivity and choose an appropriate scale or scale range for analysis.
(a)
(b)
(c)
4.3. Application of Scale-Dependent Method to Breast Tissue Samples
It is readily apparent that the curves of the averaged fractal dimensions within each tumour grade category are similar at small scales below but rapidly differentiate at larger scales Figure 6(a). In addition, there seems to be a larger difference between grades 1 and 2 at smaller scales and a larger difference between grades 2 and 3 at larger scales.
(a)
(b)
(c)
The Kruskal-Wallis analysis showed that statistically significant differences exist between at least 2 of the 3 tumour grade groups for fractal dimensions averaged over scale ranges of () and (). The multicomparison test showed that a statistically significant difference () exists between grades 1 and 2 and 1 and 3 in the scale range and a statistically significant difference () between grades 1 and 3 and 2 and 3 in the scale range (Figure 6(a)). Figures 6(a), 6(b), and 6(c) also show the fractal dimension distributions for the different grades and scale ranges in the form of boxplots. For these plots, the middle 50% of the data lie within the boxes, the lines within the boxes are the median values, the lines above and below the boxes show the upper and lower 25 percent of the data, respectively, and the crosses outside the boxes show outliers. The results illustrate the way scale range can affect the results and how different ranges can be useful for distinguishing different groups (e.g., grades 2 and 3 were not significantly different in the smaller scale range, but differences become more significant in the larger-scale range).
In a previous study using a scale range of 10 to , it was found that , , and , correlated to high, intermediate, and low survival from breast cancer [16]. The results of this study are consistent with the previous finding, as higher tumour grade is also correlated to poorer survival, and grades 1, 2, and 3 correspond to similar FD ranges (Figure 6). It is important to note that any comparisons with other results reported in the literature will only be meaningful if a similar scale range is used for finding the fractal dimension and if the tissues’ specimens are stained using the same stain and tissue preparation technique.
5. Conclusion
In this study, we have described two novel ideas for the application of fractal analysis to medical images: fractal dimension as a scale-dependent scaling parameter of a statistical distribution and an application of the integral image method for rapid evaluation of fractal dimensions. The notion of considering fractal dimension as a model-free scale-dependent parameter is a fundamental shift in perspective on investigations of fractal image analysis. By forgoing a specific model for how an image should behave, we allow ourselves to extract as much information as possible from the image. Hence, the true power of our method lies in its ability to determine the appropriate scale range or ranges that need to be analyzed using fractal methods for any particular application. In testing our algorithms on both real fractal structures and medical images, we showed the algorithm’s reliability in measuring fractal dimensions and in picking up subtle scale-dependent features in the fractal dimensions. More specifically, our analysis of invasive breast cancer tissue cores from 157 patients has shown that the ability to differentiate images of different grades of cancer depends on the scale at which images are analyzed. As tumour grade is a prognosticator for breast cancer survival, it is evident that the analysis scale has an impact on the prognostic value of a fractal analysis approach, a point which has not been systematically studied or appreciated previously. Future studies are needed to further validate and extend the breast cancer results using an independent set of tissue images.
Conflict of Interests
The authors of this paper do not have a financial relation with the commercial identities mentioned in this paper. As such, there is no conflict of interests to report.
Acknowledgments
This study was supported by a grant from the Alberta Cancer Foundation and Alberta Innovates-Health Solutions.