Fractal, multifractal and lacunarity analysis applied in retinal regions of diabetic patients with and without non- proliferative diabetic retinopathy

The aim of this study was to evaluate the blood vascular network of the retina and its regions in diabetic patients with and without signs of mild early diabetic retinopathy by using number of bifurcation points, fractal and multifractal methods. We used 33 segmented images of retinographies, 28 corresponding to the retinographies that did not show any sign of diabetic retinopathy and 5 diagnosed with non-proliferative diabetic retinopathy (NPDR). The segmented images were obtained from the DRIVE (Digital Retinal Images for Vessel Extraction) database. The segmented images of retinal blood vessels were skeletonized and fractionated into nine equal regions. Later, the skeletonized images were assessed by using the following methods: number of bifurcation points, box-counting dimension, information dimension, lacunarity parameter and multifractal analysis. The retinas and their regions of control group (diabetic without NPDR) were statistically compared to those of the NPDR by using the Z-test and the Mann-Whitney test. The multifractal analysis showed that skeletonized images of retinal vessels and its nine regions follows a multifractal behavior. In our hands, neither the fractal method nor the number of bifurcation points disclosed statistically significant differences (p>0.05). No any difference was observed in the retinal vascular network between diabetic patients with or without NPDR. Correspondence to: R.A. Nogueira, Department of Animal Morphology and Physiology, Rural Federal University of Pernambuco, Dois Irmãos, 52171-900, Recife, Pernambuco, Brazil, Tel: (+55) 81-3320 6076; E-mail: ran.pe@terra.com.br


Introduction
Among the ophthalmopathies, the diabetic retinopathy is one of the causes of vision impairment and blindness [1]. This secondary microvascular complication of diabetes mellitus occurs due to the hyperglycemia that promotes structural and functional alteration of retinal capillaries [2]. The early stage of retinopathy is termed as non-proliferative diabetic retinopathy, a disease characterized by microaneurysms, hemorrhages and capillary closure [2,3]. The proliferative phase is characterized by neovascularization, increasing of ischemic regions, hemorrhage in the vitreous cavity and tractional retinal detachment [2,4].
The retinal vascular network owns a fractal structure, the vascular branching process presents self-similarity and scaling. A fractal object or process is characterized by the following properties: self-similarity, scaling, fractal dimension [5,6]. Several works have used the fractal geometry to study the retinal vascular network [7,8]. Doubal et al. [9] have obrtained the fractal dimension of retinal vessels in the lacunar stroke and Cavallari et al. [10] in the cerebral autosomal dominant arteriopathy with subcortical infarcts and leukoencephalopathy (CADASIL). Avakian et al. [11], Kunicki et al. [12] and Ţălu et al. [13] have used fractal methods to identify non-proliferative diabetic retinopathy (NPDR), however obtaining contradictory results.
Since the fractal dimension describes how much space is filled, but does not indicate how the space is filled by the fractal structure, the lacunarity can solve this problem by making a distinction between different objects with the same fractal dimension [14,15]. The lacunarity is a parameter that indicates the distribution of gap sizes throughout the object embedded in an image, it is able to identify different fractal structures that have the same fractal dimension [14]. This parameter has been used as a tool in the characterization of retinal vascular network. Landini et al. [16] have employed the lacunarity to identify the occlusion of the artery and retinal vein, whereas Tălu et al. [17] have used it to diagnose amblyopic eyes.
Some studies have currently characterized some objects as multifractal structures (objects that have different fractal dimensions) instead of monofractal, as regarded in the previously mentioned studies. An object is considered multifractal when its different regions have different fractal properties [18]. Researchers have demonstrated the network of vessels of the retina as a multifractal object, a fact which has been proved through generalized dimensions and singularity spectrum [15,19]. Furthermore, the singularity spectrum has the capacity to evidence the disorders in the retinal vascular architecture with diseases [19].
The aim of this study was to assess the network of blood vessels of the retina and its regions with and without signs of mild early diabetic retinopathy by using the number of bifurcation points and fractal methods, such as box-counting dimension, information dimension, lacunarity parameter (complementary tool of fractal methods) and multifractal analysis: generalized dimensions and singularity spectrum.

Retinal images
The images were obtained from the DRIVE (Digital Retinal Images for Vessel Extraction) [20] database (http://www.isi.uu.nl/Research/ Databases/DRIVE/). We used 28 retinographies (control group) which do not show any sign of diabetic retinopathy ( Figure 1A) and 5 diagnosed with mild early diabetic retinopathy, non-proliferative diabetic retinopathy (NPDR) ( Figure 1B). The images were acquired from diabetic patients between 25-90 years of age, by using a Canon CR5 non-mydriatic 3CCD camera with a 45 degree field of view (FOV). Each image was captured by using 8 bits per color plane at 768 by 584 pixels. The FOV of each image is circular with a diameter of approximately 540 pixels and each image has been JPEG compressed.

Skeletonization and separation by region
The manually segmented images of retinal vessels ( Figure 1C) were also obtained from DRIVE. The segmented images of vessels (565x586 pixels) were skeletonized by the software Matlab® version 7.8 (MathWorks, Natick, Ma, U.S.A). Each image was fractionated into nine equal regions ( Figure 1D) by the software Adobe Image Ready 7.0.1. Figure 1 shows the retina in gray scale, segmented, skeletonized and divided by region (nasal superior, optic disc, nasal inferior, superior, macular, inferior, superotemporal, temporal and inferotemporal).

Counting of bifurcation points
The bifurcation point is where a blood vessel originates other vessels. The number of bifurcation points is a method that informs about the amount of branching, showing the vascular multiplication. The procedure proposed here is an optional method to measure the blood vascularization degree as another morphometric parameter to evaluate the vascular architecture (density, length, diameter of vessels). The bifurcation points were determined from skeletonized images of the retinas and their respective regions for both groups. The identification process was confirmed in the original retinographies to reduce the possibility of quantifying artifacts or overlapping of vessels.

Fractal dimension methods
Two methods were used to calculate the fractal dimension of retinal blood vessels, the box-counting dimension (D bc ) and the information dimension (D inf ) by the software Benoit 1.3 Fractal Analysis System (Trusoft, St. Petersburg, FL, USA). For the box-counting dimension (D bc ), the skeletonized image was covered with a number of boxes (N (r)) containing at least one pixel of the image. The procedure was repeated with boxes of different sizes and plotted in a double log graph of N (r) in relation to the sides of boxes r [21]. The slope of this relationship with inverted signal is the box-counting dimension: in which ε is an infinitesimal variation in the box sizes. The calculations of D bc used 19 sets of different size boxes, the length of the largest box side being 270 pixels and the reduction coefficient of the box size 1.3.
In the information dimension (D inf ), skeletonized images were also covered by boxes, but taking into account the relative probability of occupancy of the elementary boxes used to cover the fractal object: is called the Kolmogorov entropy. N is the number of boxes, m i = M i /M, and M i is the number of pixels in the n th box, M is the total number of pixels on the fractal object, r is the side of boxes and ε is a infinitesimal variation in box sizes [21]. The calculations of D inf used 8 sets of different size boxes, the length of the largest box side being 270 pixels and the reduction coefficient of box size 2.0. Figure 2 shows the box-counting dimension graph ( Figure 2A) and information dimension ( Figure 2B).

Lacunarity parameter
To evaluate the lacunarity parameter of the images of the retinal vessels, we used the software Image J (Wayne Rasband, National  Institutes of Health in Bethesda, Maryland, USA) with the FracLac plug-in (A. Karperien -Charles Sturt University, Australia). Lacunarity was obtained by measuring the gap dispersion inside an image; in other words, it was related to the pixel distribution of an object in an image. The quantification was achieved as in the box-counting method. In this case, however, different directions to the set of boxes (g) was also used.
The mean value to the lacunarity was calculated as follows: where σ is the standard deviation and μ is the mean of pixels per box with size r, in a box-counting at a direction g, n being the number of box sizes. The sum is done over all values of r and g.

Multifractal analysis
Software Image J with FracLac plug-in was also used to calculate the multifractality of retinal vasclarization. The multifractal structure is characterized by obtaining the generalized dimension D q , which is related to a value of q. Variable q is the exponent that expresses the fractal properties in different scales to an object, ranging between -∞ and + ∞. The plot of D q versus q is generally sigmoidal and decreasing to a multifractal structure. Some studies consider determined values D q that is D 0 , D 1 and D 2 , which depict the multifractality of an object when condition D 0 ≥ D 1 ≥ D 2 is satisfied. Values D q depict the multifractal object that can be compared by the single fractal dimension method, so D 0 can be considered as capacity dimension, D 1 can be related to the information dimension and D 2 to the correlation dimension [19]. In our study, values were generated from D -10 to D +10 , in other words, q values ranged between -10 and +10, where all dimensions were statistically tested. D q was calculated as follows: τ can be defined as: where P i (r)=M i /M o is density, M i is the number of pixels within the i th box and M 0 is the number of pixels for all image.
i P ∑ is the density for all boxes (i) at a determined scale r.
Another way to calculate the multifractal spectra is through the relationship between parameters f (α) versus α, where N(α) is the number of boxes for which probability P i (r) of finding a pixel within a given region i scales as is the fractal dimension to the all regions with singularity strengths between α and α + dα, where α takes on values between -∞ and + ∞.
These variables were obtained by Legendre transformation from τ into q, as described below: and ( ) ( ) Also, α (q) and f (α(q)) can be calculated by using the following equations: P i q (r) is the probability of pixels in the i th box for the exponent q. For a multifractal object, plot α versus f(α) fits a parabola with concavity turned down.
From values α, it was possible to calculate the extension of singularity length Δα and the curve asymmetry of singularity spectrum. Extension of singularity length Δα is defined as follows: α max (highest value of α) and α min (lowest value of α) indicate the fluctuation of minimum and maximum probability of pixels, respectively. The higher the Δα, the larger the probability distribution; furthermore, the multifractality is stronger and the pixel distribution of image is more complex [22,23].
The curve asymmetry of singularity spectrum (A) was calculated according to the following expression: The curve of singularity spectrum is symmetric to A=1, the curve is left-skewed if A>1 and is right-skewed if A<1. When the spectrum is left-skewed, it means that there is stronger presence of high fractal exponents and significant fluctuation; otherwise, it indicates the domain of low exponents and slight fluctuation [23].

Statistic analysis
The Shapiro-Wilk's test was used to calculate the number of bifurcations, fractal dimensions, lacunarity parameters, generalized dimensions, Δα and parameter A (curve asymmetry of the singularity spectrum). The Shapiro-Wilk's test was used to choose between a parametric test or nonparametric, in order to make a comparison between the different regions and the whole retina with and without NPDR within the two groups. The Z-test was used when the group had a normal distribution and the Mann-Whitney test for a non-normal distribution. The Z-test, here, is a method to inform if the distribution of group with NPDR had the same behavior of the distribution of group without NPDR, this test is acceptable to conditions in that one of the groups possesses a low number of samples while the other has the highest number of samples. Table 1 shows the number of bifurcation points, in which p is the value of the significance level to the Z-test or the Mann-Whitney test (with asterisks). For regions such as optic disc, nasal inferior and temporal in which the control group did not present a normal distribution, the Mann-Whitney test (nonparametric test) was used.

Number of bifurcation points
For other regions and the whole retina in which the control group presented normal distribution, the Z-test was used. The statistical tests did not show any difference in the bifurcation points between the control group and the NPDR group (p>0.05).

Box-counting dimension and information
The first step to make the statistical analysis of fractal dimensions of the two groups was to test whether the segmentations represented the original images of the retinas accurately. We selected ten segmented images of the DRIVE made by one observer, carried out the skeletonization for the control group and compared it to the skeletonization of the same ten images segmented manually by another observer. Figure 3 represents the box-counting dimensions (3A) and information (3B) to segmented and skeletonized images, respectively; t-student test showed that for both observers there were no significant difference neither between the segmented images (p=0.82 for boxcounting dimension and p=1.0 for information dimension) nor the skeletonized ones (p=1.0 for box-counting dimension and p=0.68 for information dimension). These results showed that the manual segmentation is a reliable method and can be considered as gold standard. Table 2 shows the fractal dimensions (box-counting dimension and information) of the control group and NPDR. In this case, neither the whole retina nor its different regions in the control group displayed difference to the NPDR group (p>0.05). Table 3 depicts the lacunarity values. Some regions of the two groups were submitted to the Mann-Whitney test (values with asterisks) and others to the Z-test. Statistically, there was no difference between the lacunarity parameters for the control group and NPDR (p>0.05). The Shapiro-Wilk's test was used to test the normality.

Multifractal analysis
There was a statistical difference in dimension D -1 in the inferotemporal region (p=0.04 to the Z-test), whilst other generalized dimensions of the same region corresponding to both groups were not significantly different. The generalized dimensions of other regions and the whole retina did not show any statistically significant difference. Table 4 presents the mean and standard deviations of following generalized dimensions for both groups: D -10 , D 0 , D 1 , D 2 and D 10 . Figure 4 shows two graphs that represent the mean and standard deviations of generalized dimensions (D q ) versus variable q for the whole retina and for the inferotemporal region in which there was a statistical difference in dimension D -1 . The plots followed a sigmoid fit, outlining a multifractality for blood vascular network of the retina and its regions. Each graph has two plots, one that depicts the control group and the other showing the retinopathy group. Figure 4A represents D q versus q for the whole retina in the two groups (control and NPDR), showing a behavior with less oscillation (low standard deviation) compared to Figure 4B (inferotemporal region). Figures 5A  and 5B represent the graph of f (α) versus α to the whole retina and the inferotemporal region respectively, in which each graph presents the group with and without NPDR. The graphs showed the typical behavior of multifractal structures, a parabole with concavity facing down. Table 5. Values of parameter Δα represent the mulifractality of the blood vascular network of the retina and its regions. Tests revealed no statistically significant differences for parameter Δα between the control group and the one with NPDR. In relation to the mean values with the respective standard deviations of parameter A , the retinas and all their regions for the two groups were observed to have values of A <1, or, in other words, the curve of singularity spectrum f (α) versus

Discussion
The non-proliferative diabetic retinopathy is characterized by microaneurysms, hemorrhages and capillary closure [2,3]. The rise in hemodynamic alterations and compensatory mechanisms due to the presence of such signs may promote changes in vascular architecture. In this paper, we have used methods that are capable of identifying the geometric complexity of the blood vascular network.
The fractal dimension is a statistical descriptor of the space-filling pattern and density serving as a tool capable of evaluating the vascular development and therefore it can help in the diagnosis of diseases related to the vascular disorders [7,24,25]. The retinal vascular network has a fractal dimension that ranges between 1 and 2, as being close to 2 indicates a more complex network or higher vascular density. For instance, Doubal et al. [9] have associated a low value of fractal dimension of retinal vessels to the lacunar stroke subtype. Also, according to Cavallari et al. [10], the decrease in fractal dimension is related to the reduction in complexity of retinal vessels, which reflects in the alteration in brain microvessels in patients with CADASIL.
To identify possible vascular alterations not disclosed by fractal dimensions, the parameter of lacunarity, which has the capacity to recognize different fractal structures with the same fractal dimension, was employed, describing how the pixels are organized in the full image. This is a way to measure heterogeneity, as well as the degree of invariance to change of the fractal object in an image. Whereas the fractal dimension indicates how much space is filled by vessels, this parameter (reflected in the distribution of gaps) represents how the vessels fill the space wherein they are embedded [15]. This parameter has been used to discover alterations in the arteries and retinal veins [16], as well as to diagnose retinas with amblyopia [17]. The lacunarity parameters obtained were not statistically different as regards the normal state and the disease, thus ratifying the results obtained by using the number of bifurcation points, box-counting and information dimensions.
In the multifractal analysis, only one related value D -1 of the inferotemporal region showed statistical difference between the two groups. We can say that the single exponent q=-1 is not enough to characterize the retinal vascular architecture with NPDR. All generalized dimensions to retinas and their other regions showed no any differences. The multifractal analysis ensures more information about the behavior of the network of blood vessels once it reveals several values that estimate simultaneously the level of geometric complexity or filling of space by the vessels. The graphs in figure 4 present a  The Mann-Whitney test was used in the regions with asterisks and the Z-normal test was used in the others. The Shapiro-Wilk's test was used to test the normality.  behavior that defines them as a multifractal structure, represented by the descending sigmoid curve of D q spectrum (generalized dimensions) versus q for the two groups. In case the object presents constant values for fractal dimensions, graph D q versus q should be a straight parallel to q axis, and the object must be considered a monofractal structure.

Retinal region
Singularity spectrum f (α) calculates the fractal dimension of the subset of pixels in an image described by a particular exponent, referring to the relative dominance of various fractal exponents involved in the structure [26]. Figure 5 showed that spectrum f (α) versus α also represents the multifractality of retinal vessels for both groups; in this case, the spectrum is a parabola with concavity facing down [27]. Stŏsić and Stŏsić [19] have shown that the singularity spectrum (graph f (α) versus α) of retinal vascular network follows a multifractal geometry for both normal and pathological states; they have also shown that a group of retinas with different diseases presented lower singularity compared to normal retinas. In our case, we showed a group diagnosed just with NPDR and that there was no distinct singularity spectrum in relation to the diabetic group without the NPDR ( Figure 5). In order to assess the multifractality of the vascular networks of retinas and their regions, length Δα [22,23,26] was used, which was not statistically different when comparing the two groups as observed in Table 5. Another way to analyse the singularity spectrum is through its curve asymmetry [23]. As already mentioned, this parameter was not statistically different within the two groups. Moreover, both groups showed values of low fractal exponents (A<1).
The condition D cap ≥ D inf ≥ D co has been established for each method that these exponents represent in measuring the fractal processes [28].
In terms of generalized dimensions, this condition can be described as follows: D 0 ≥ D 1 ≥ D 2 [15,19,29]. Each retina in both groups presented this condition, but when the regions were analyzed, some images from those regions lost this criterion, even for decreasing D q values with the increase in exponent q. The superior region was the one with the highest number of images (7) not following such criterion. All images of the inferior region, macular and nasal superior are within this condition. Nevertheless, we can state that the vascular network of retinal regions presents multifractality, not only the vascular network of the whole retina, but its regions can also be considered a superposition of monofractal structures [19,30].
The fractal dimensions obtained by the box-counting and information methods are not ruled by the condition stated by Grassberger and Procaccia [28] in that D cap ≥ D inf . Likewise, results obtained by Mendonça et al. [31]; Kunicki et al. [12] and Voinea and Popescu [32] were also found to be out of that condition. Considering that the mass-radius dimension of a structure is higher than the boxcounting dimension (capacity dimension) [31,32], Family et al. [7], having obtained a higher correlation dimension than the mass-radius one for retinal blood vascularization, the principle D cap ≥ D inf ≥ D co was contradicted.
Some researches have reported alterations of retinal vascularization in individuals with different states of diabetes, Daxer [33] has utilized the fractal geometry to characterize quantitatively the formation and regression of neovascularization in the diabetic retinopathy. Cheung et al. [34] observed that the increase in the fractal dimension of the retinal vasculature is associated with early NPDR signs in young individuals with type 1 diabetes. Ţălu et al. [13] also found the fractal dimension values slightly lower to the retinal images of normal patients than the values to the patients' images with NPDR. However, Lim et al. [35] have reported that the fractal dimension of retinal vasculature was not associated with incident early diabetic retinopathy in children and adolescents that possess the type 1 diabetes. In this work, the diabetic patients despite not having the signs that characterize mild early diabetic retinopathy showed vascular network geometry similar to the patients with the retinopathy, according to the results revealed by fractal methods (box-counting dimension, information dimension, generalized dimensions, singularity spectrum and lacunarity    parameter) and the number of bifurcation points in skeletonized images of the retinal vascular network. However, contradictory results have been described in other works. Avakian et al. [11] performing fractal analysis in skeletonized images from 60-degree fundus fluorescein angiography, observed that the density values of vessels in normal retina macular region was higher than the vascular density of the retinal macular region with NPDR. Our study corroborates the results obtained by Kunick et al. [12] that showed no difference between retinal images of patients with and without NPDR, utilizing fractal dimension analysis (box-couting and information dimension) in segmented images of the retinal vascular network. Diabetic patients, despite not having the signs that characterize a mild early diabetic retinopathy, show a vascular geometry similar to the patients with the disease, according to the results revealed by fractal methods (boxcounting dimension, information dimension, generalized dimensions, singularity spectrum and lacunarity parameter) and the number of bifurcation points. We can suggest the existence of vascular network alterations before to arise the NPDR, since the retinal fractal dimension in normal individuals is different than in diabetic patients [36].

Conclusion
According to our fractal analysis (box-counting dimension, information dimension, generalized dimensions, singularity spectrum and lacunarity parameter), as well as the measure of number of bifurcation points did not disclose geometrical alterations in the retinal vascular network for either group of diabetic patients, with or without non-proliferative diabetic retinopathy. In addition, we have observed that not only the vascular network of whole retina, but also its several regions follow a multifractal behavior.