Automated detection of proliferative retinopathy in clinical practice.

Timely intervention for diabetic retinopathy (DR) lessens the possibility of blindness and can save considerable costs to health systems. To ensure that interventions are timely and effective requires methods of screening and monitoring pathological changes, including assessing outcomes. Fractal analysis, one method that has been studied for assessing DR, is potentially relevant in today's world of telemedicine because it provides objective indices from digital images of complex patterns such as are seen in retinal vasculature, which is affected in DR. We introduce here a protocol to distinguish between nonproliferative (NPDR) and proliferative (PDR) changes in retinal vasculature using a fractal analysis method known as local connected dimension (D(conn)) analysis. The major finding is that compared to other fractal analysis methods, D(conn) analysis better differentiates NPDR from PDR (p = 0.05). In addition, we are the first to show that fractal analysis can be used to differentiate between NPDR and PDR using automated vessel identification. Overall, our results suggest this protocol can complement existing methods by including an automated and objective measure obtainable at a lower level of expertise that experts can then use in screening for and monitoring DR.

Introduction

Diabetes is a significant and growing global public health problem, significant not just in itself but also for its associated complications (Taylor and Keeffe 2001). One common major complication of the disease is diabetic retinopathy (DR). Currently, between 22% and 36% of people with diabetes have some form of DR, and in one-third of them, the disease has progressed to a severe, vision-threatening stage (NHMRC 1997b).

From a clinical perspective, DR is generally considered to progress through two broad stages that correspond to respectively lesser and greater levels of visual impairment. From a pathological perspective, DR can be considered in terms of changes in retinal vascularization, where the two stages are classified with respect to the absence or presence of abnormal neovascularization as either (1) nonproliferative (NPDR) or (2) proliferative (PDR). These categories are usually further graded by severity and the characteristics of lesions present within each. (ETDRS 1991; NHMRC 1997a) During the earlier NPDR stages, microaneurysms and hemorrhages typically occur, the macula may swell as small retinal vessels leak fluid, and cotton wool spots may also appear. Over time, retinal vessels may eventually close, leading to further vascular abnormalities, more hemorrhaging, and the formation of exudate. During the PDR stages, the pathology progresses to include abnormal neovascularization. In essence, as the disease continues, the number of non-perfused capillaries and resultant ischemia increase and this eventually leads to new vessel growth, where the new vessels are abundant but dysfunctional and, rather than resupplying a starving retina, introduce further problems largely because they are fragile and often misplaced, with severe consequences. Neovascularization may extend into the vitreous, for example, where it can cause hemorrhaging and contribute to retinal traction and detachment, leading to irreversible blindness. (Osborne et al 2004) Retinal holes may also develop near the proliferation. (Kanski 1989; Yam and Kwok 2007)

We examined ways to characterize changes in retinal vasculature in DR, which is why this paper naturally focuses on that aspect of the pathology. It is important to emphasize, however, that DR is more than simply disordered retinal vascularization, and that the precise mechanisms at every stage have not yet been identified. Indeed, the knowledge base about underlying pathophysiologic changes in DR is currently changing dramatically and researchers are being challenged to approach the disease with new paradigms and perspectives. (Olson et al 1997; Krady et al 2005; Antonetti et al 2006)

As the knowledge base continues to grow, investigators are identifying potential leads that may eventually yield treatments (eg, changes in growth factors affecting neovascularization), but currently there is no cure for DR. Despite that it is incurable, it is not currently unalterable. There are excellent opportunities for intervention at the incipient stages (such as improving glycemic control) that have the potential to reverse or prevent further progress of the disease. (De La Cruz et al 2004) Moreover, there are limited opportunities for intervention even at advanced stages. If disease is detected early enough, for instance, laser treatments can diminish visual loss (albeit in a trade off which may introduce other visual problems). (Doft and Blankenship 1984; Ferris 1993; Icks et al 1997; Lövestam-Adrian and Agardh; Lee et al 2001a; Bek and Erlandsen 2006)

Diabetic retinopathy assessment

Although early detection and ongoing monitoring of DR have great potential to save vision, paradoxically, much evidence has amassed indicating that patients often do not get the necessary screening and monitoring to optimally apply the interventions that are available. The problem exists globally, and many barriers to screening affecting different populations to different degrees have been identified, including factors like inadequate access to care, patient misconceptions about the value of regular eye examinations, costs of visits to specialists, distance required to travel, and cultural reasons (eg, indigenous people may remain in their communities rather than seeking health advice in larger urban centers). (Mukamel et al 1999; McGlynn et al 2003; Paz et al 2006)

To address this serious problem, goals have been set nationally and internationally to increase the proportion of people being screened for DR (the Australian national target, for example, is 80%). (NHMRC 2001) Many programs with the goal of increasing awareness in the community have been implemented accordingly, but to solve the underlying problem – that is, to globally reduce the numbers of individuals seriously affected by PDR – another important goal is to develop and implement new or complementary models for effective assessment, treatment, and monitoring of DR. (Lee et al 2001)

One particular step towards ensuring that early and regular eye examinations become routine for people with diabetes regardless of geographic, economic, cultural, etc., considerations, is to simplify current procedures in order to hand over responsibility to a more accessible level of technical training (Taylor and Keeffe 2001). The rising use of telemedicine has started to shift part of the burden of DR screening to trained but non-expert health personnel around the world and has stimulated opportunities to develop the technology around digital retinal imaging (eg, incorporating digital images in an electronic patient record can help in ensuring diagnostic accuracy and monitoring and measuring outcomes). In addition, new methods of acquiring and assessing information about retinal pathology are being developed at a rapid pace (Hutchinson et al 2000; Neubauer et al 2003; Conlin et al 2006a; Conlin et al 2006b; Taylor et al 2007).

Whereas some studies so far having non-specialist health professionals identify PDR have found detection rates to be no better than 50%, with the success rate decreasing with features such as eye obstruction and early stages of proliferation, other studies suggest problems with the identification of PDR can be overcome using digital retinal photography and automated computer-based processes (Frame et al 1998; Maberley et al 2003; Olson et al 2006; Soares et al 2006).

Computer-based diabetic retinopathy detection

Methods for objectively identifying and quantitating patterns in retinal vasculature involve two main steps. The first is identifying blood vessels within images or segmentation. Several segmentation methods that identify retinal vessels in either fluorescein or color digital images of the fundus have been reported in the literature, including the use of mathematical morphology, matched filters, threshold probing, supervised classification, deformable models, and tracking (Englmeier 1997; Gao et al 2000; Hoover et al 2000; Zana and Klein 2000; Jiang and Mojon 2003; Kirbas and Quek 2004; Nain et al 2004; Staal et al 2004). In the work reported here, we identified retinal vasculature patterns in digital images of fluorescein angiograms using two-dimensional continuous wavelets. Recent advances using wavelet-based techniques have shown particular promise as these methods remove noise from images, are accurate, and can provide additional parameters that can be used to classify PDR and follow disease progression (Leandro et al 2003; Cornforth et al 2005; Cree et al 2005a; Cree et al 2005b; Soares et al 2006).

The second step is measuring the vessel pattern obtained. Many features that can be extracted from digital images of blood vessels (eg, density, total length, and curvature) have been used to quantitate patterns in retinal vasculature. In the present study, we assessed complex branching features of the retinal vascular tree that defy description with traditional Euclidean geometry but reflect principles of fractal geometry (Daxer 1993a; Landini 1996; Cornforth et al 2002).

Fractal geometry in proliferative retinopathy

Fractal geometry was first introduced in 1960 by Richardson to describe coastlines. He pointed out that depending on the scale used to measure a coastline, the total length arrived at would vary. Specifically, with a decrease in scale would come a nonlinear increase in length (Richardson 1960). What Richardson described was a general phenomenon applicable to many types of objects and patterns whereby as the scale or magnification at which something is observed changes, the detail measured, be that length, area, or volume, changes nonlinearly. This basic principle of fractal geometry is illustrated in Figure 1.

Figure 1

A 32-segment quadric fractal. The basic pattern shown in the upper box can be seen repeated throughout the structure at different scales. When the figure is viewed at a scale factor of 8 relative to the previous view, the pattern repeats itself with 32 pieces, each 1/8 the length of the original of which it is a part.

Richardson’s general idea of fractal geometry was popularized by Mandelbrot and developed into fractal analysis, which is in use today in many parts of pure and clinical science (Mandelbrot 1993; Losa and Nonnenmacher 1996). Fractal analysis rests on the general point that the relationship between the resolution or scale at which an object is measured and the measured outcome can be quantitatively expressed as the fractal dimension (DF) of the object (Liebovitch and Sullivan 1987; Bittner et al; McGinley et al 1994; Gitter and Czerniecki 1995; Iannaccone and Khokha 1995; Bernard et al 2001; Feltrin et al 2001; Goldberger et al 2002; Masters 2004; Naschitz et al 2004).

To elaborate, the similarity dimension (DS), is one instance of a DF. Its derivation is perhaps easiest understood starting with simple, ideal shapes. A straight line segment, for example, when measured using pieces 3−1 times its size, is found to be made of 3 such scaled pieces, in a relationship that holds through an infinite number of such scalings. Formally, a power law defines this trivial (ie, simple) but infinite self-similar scale-invariance, the DS being the exponent in that power law as in the equation below:

In this equation, N is the number of equal pieces that resemble the original form when r is the ratio or scale applied to the object to make the equal pieces. Manipulating this equation gives a solution for the DS:

The DS, thus, quantitates complexity, the detail or number of self-similar parts changes as scale changes (eg, the DS for the line described above is log 3/log 3 = 1.0 and for Figure 1 it is log 32/log 8 = 1.67).

Various measures related to the DS are readily calculated by different types of morphological image analysis software. These approximations are not always numerically identical, but they are generally relatively similar. The box counting dimension (DB), which we used in this paper, is calculated using a method that approximates the D by measuring patterns in digital images using squares of decreasing size, r, instead of actually counting new parts at decreasing scales. In this case, N is the number of squares of each size, r, required to cover the image, and the DB is the slope of the log-log regression line for these two values. More details on obtaining DFs are provided elsewhere (Feder 1988; Smith et al 1989; Peitgen et al 1992; Iannaccone and Khokha 1995).

Correlation dimension

Another numerically simple approximation for the DF that can be applied in clinical science and ophthalmology is the correlation dimension (CD) (Grassberger 1983; Family et al 1989; Daxer 1992; Masters 2004). It is a probabilistic dimension that allows the differentiation of true stochastic processes from deterministic chaos. The CD is based on the integral function C(ɛ) that defines the probability that two arbitrary points on an orbit are closer together than ɛ and is expressed as:

Local connected fractal dimension

Yet another dimension relevant to retinal vessels is the connected dimension. True fractals are scale invariant, meaning one scaling factor holds for all magnifications at which the fractal is measured. Biological objects, in contrast, are generally only statistically invariant, meaning they scale consistently over only a defined range (eg, whereas the scaling in the pattern in Figure 1 theoretically goes on infinitely, scaling in the retinal vasculature tree is limited by the sizes of the smallest and largest vessels) (Voss and Wyatt 1991; Dollinger et al 1998; Fernandez and Jelinek 2001). In addition, biological structures are generally influenced by many processes, both microscopic and macroscopic, that can create local variation in scaling (Vicsek and Vicsek 1997). The retinal vascular tree, for instance, is subject to glucose and insulin levels, nutrient availability, changing pressures and pressure gradients, the constitution of its external environment, interactions with cellular components, growth factors, genetic and developmental factors, the confines of the physical space within the eye, etc (Antonetti et al 2006). As such, whereas the vascular tree may unto itself bear a certain scaling dimension, it may also carry areas of unique subscaling within it.

One way to measure the complexity of any number of subsections within an object is to analyze local dimensions (Voss and Wyatt 1991). In contrast to global fractal analyses that generate single values for entire images, local analyses generate arrays of local dimensions that can be visualized as distributions (eg, in color-coded images and histrograms) or summarized using the mean and standard deviation of the local array. The local connected dimension (D) is a local dimension that is pertinent to DR because it takes into consideration the degree of contiguous branching within a particular area. To elaborate, for digital images, a D can be calculated for each significant pixel, P(x, y), corresponding in DR to parts of the retinal vasculature. First, the connected set of pixels within some defined area around and all traceable to P is identified, then the D is determined based on that set. This is generally done by centering on P a series of concentric sampling elements (eg, squares or circles) of decreasing size, r, and counting the pixels that fall in both the sampling element and the connected set at P. The D is then calculated as the slope of the regression line for the double logarithmic plot between r and the number of pixels counted at r (see Figure 2) (Landini et al 1995; Fernandez and Jelinek 2001; Karperien 2007).

Figure 2

Calculating the local connected fractal dimension (D). Each pixel, P(x, y), that represents part of the vascular tree in a retinal vessel image has an associated local connected set for every distance that can be defined around that pixel. As illustrated at the bottom left of the figure, the local connected set includes only the pixels that can be traced along a continuous path from P(x, y). As illustrated at the bottom right, this set is used to calculate the D at (x, y). The image shows the data gathering process in which concentric squares of decreasing size are laid on the connected set and a local dimension is determined accordingly (see text).

The D has so far only been used to differentiate between normal and occlusion angiograms in manually traced retinal vascular patterns (Landini et al 1995). Our work expanded on this to use the D to differentiate PDR from images with no proliferative changes based on variation in the complexity of patterns extracted from the retinal vasculature tree using automated segmentation procedures.

Methods

Twenty seven fluorescein labeled images from patients with varying degrees of retinopathy and eye disease ranging from no eye disease to PDR, both with and without neovascularization, were collected from patients attending the Albury Eye Clinic for eye assessment. Fluorescein was injected under supervision of an ophthalmologist and images were acquired using a Topcon digital camera (1024 × 1024 pixels) combined with the Image 2000 software package. This was a retrospective study of patients attending the clinic who required fluorescein labeling to determine the progression of their DR. Digital images were available for only one year as the clinic had only recently changed from using 35 mm slides to digital images. Two sets of skeletonised, binary images for analysis were prepared from the initial images. For one, the pattern of the retinal vasculature in each image was traced manually using a WACOM Graphire tablet to generate intermediate (nonskeletonised) images, and the resulting patterns were skeletonised using ImageJ (Rasband 1997–2007); for the other set, an automated segmentation was applied and then the result skeletonised.

Automated blood vessel segmentation

Blood vessels were segmented as described in detail elsewhere (Soares et al 2006). The continuous wavelet transform (CWT) is a powerful tool to analyze non-stationary signals (Costa and Cesar Jr 2001). Vessel-pixels responded strongly to the two-dimensional Gabor wavelet, because the wavelet is directional and could be specially tuned for vessel detection (Antoine et al 1993). The wavelet was superimposed onto each pixel of the image at various angles and scales. In order to detect the blood vessels, for each scale value chosen, the transform was calculated over the range of 0–170 degrees, at steps of 10 degrees, and the feature space was updated with the maximum value at each pixel position. The magnitude of the parameters associated with the Gabor wavelet transform were empirically determined in order to reach the best matching between wavelet and vessels. Once determined, the parameter configuration did not have to be changed from image to image. The modulus of the wavelet transform using four different wavelet scales – chosen as to span all possible vessel widths – at each pixel was taken as a pixel feature, as well as the original gray scale intensity, yielding five wavelet coefficients per pixel to compound the final feature space.

A normal transformation was applied to all features to obtain dimensionless values, which allows a comparison between features and avoids mistakes in classification. With the normal transformation, all features present zero mean and unitary standard deviation with respect to the training set.

The final segmentation was obtained by classifying the original input image pixels into two classes, namely vessel-pixels and non-vessel pixels, according to the supervised classification approach (Costa and Cesar Jr 2001). A Bayesian classifier was adopted in which class likelihoods were described using Gaussian mixture models, providing a fast classification while still allowing complex decision surfaces. The class priors were estimated by the fraction of each class’s pixels present in the training set composed of labeled samples, while the distribution parameters for each class’s Gaussian mixture model were estimated from the training set through the Expectation-Maximization algorithm (Theodoridis 1999). The training set was obtained from manually segmented fundus images, providing us with the 2 classes (ie, vessels and non-vessels). Once the normalized training set from the hand-drawn vascular trees was obtained and the classifiers parameters were estimated, the classification itself took place.

The classifier output is a binary image with pixels labeled as vessel or non-vessel.

Some misclassified pixels appeared as undesirable noise in the output and as only boundaries were classified for some vessels, post-processing was necessary. The postprocessing operations applied were ‘area open’ to eliminate small noisy components, ‘dilation’ and ‘area close’ to fill the vessels, followed by ‘skeletonisation’ to extract the vessel pattern and resulting in a 1-pixel-wide connected structure. The multiscale skeletonisation algorithm based on exact dilations was applied in this last step (Costa and Cesar Jr 2001).

Determination of fractal dimensions

The CD was calculated as previously discussed in the literature (Family et al 1989). This procedure leads to a graph C(ɛ) versus ɛ from which a log-log plot-based line fitting is able to estimate the CD. The extremities of the linear portion of the log-log slope were determined by taking the wavelet transform of the log-log plot using the 3rd derivative of the Gaussian as mother wavelet and finding its zero-crossing points (Cesar and Jelinek 2003).

The DB and D were calculated using the morphological image analysis software FracLac for ImageJ (Karperien 2007). The DB was calculated using a relative maximum box size of 50% of image size and a minimum box size of 2 pixels. The D for each pixel of the skeletonised patterns was calculated using a connected set within a square of width equal to 10% of image size, which was also the value used for the maximum box size, and a minimum box size of 2 pixels. To obtain a single value equivalent to the global dimension, the average of all D values (μD) was also calculated for each image.

Results

Receiver operating characteristics provided an indication of precision when using automated segmentation of the 27 images compared to the manual drawings used as the gold standards (Figure 3).

Figure 3

The range for area under the curve (AUC) was between 0.7185 to 0.9492, with an average and standard deviation of 0.8832 ± 0.0397. The accuracy values average and standard deviation was 0.903 ± 0.018.

Global DB and μD results for all images

Figure 4 shows the global D and the μD for all images, PDR and control combined, for each segmentation method. For the entire sample of manually segmented images, the mean global DB was 1.61 ± 0.05 and for the entire sample of automatically segmented images, it was lower, 1.52 ± 0.06. The mean μD for manually segmented images was 1.21 ± 0.04 and for automatically segmented images was 1.14 ± 0.03.

Figure 4

Results of the box-counting analysis and the local connected dimension analysis using all 27 images, for the manually segmented group and the automatically segmented group. Image names prefaced by PDR are the proliferative retinopathy images, and images with no letter in front are the control images. Box Counting: The figure shows the global box counting dimension (DB) for each image. Local Connected Dimension: The figure shows the mean local connected fractal dimension (μD) for each image, calculated by averaging the D values for each pixel in each image.

The lower values for the automatically segmented images appears to be an effect of noise that obscures essential features of the vascular tree and leaves gaps in the final image when using the automated method. The ultimate effect is that some areas that are traced by the manual method tend to be “erased” by the automated method. Further evidence for this interpretation is found in the point that the average number of foreground pixels was lower in the automatically segmented images compared to the manually segmented ones. Inspection of the images revealed that common areas of differences were areas with fine vessels at the ends of branches (these were sometimes omitted in the automated method), areas where fluorescein was unevenly distributed (see Figure 5), and areas with noteworthy hemorrhages or laser scars.

Figure 5

Noise affected automated segmentation. Shown here are two original fluorescein images (bottom) and their respective segmented and skeletonised vessel patterns (top). In image 16459__3_ on the left, fluorescein was relatively evenly distributed and the two patterns yielded similar global DBs (automated = 1.54 and manual = 1.59; see Figure 4 also). In image 20_201_ on the right, however, fluorescein was unevenly distributed, which caused some areas that were traced in the manual group to be left blank in the automated group, and the global DBs consequently differed more (automated 1.52 and manual = 1.66).

Despite these difficulties, on the whole the automated segmentation method identified the retinal vascular tree very well. An automated segmented blood vessel pattern is illustrated in Figure 6. The fluorescein image clearly shows a homogenous distribution of the dye throughout the vessel pattern and was correctly segmented using the wavelet approach.

Figure 6

Automated segmentation of a retinal blood vessel pattern. A. The original fluorescein labeled image 18_393__. B. The automatically segmented and skeletonised vascular pattern from the image. C. The automatically segmented pattern from B, dilated to a uniform diameter and shown overlain on the original image as a dark vessel pattern over the fluorescein-labeled (white) vessels to illustrate the correspondence of the automatically selected pattern with vessels in the original.

Comparison of pathological status using means

Tables 1a and 1b show the results for the four feature parameters we analysed for the automated and manually segmented images. As shown there, the only measure that was significantly different for PDR and control images in both the manually and the automatically segmented images was the μD.

Table 1a

Results of t-tests for automatically segmented images

	Control (mean ± stdev) N = 11	PDR* (mean ± stdev) N = 16	P value
DB	1.53 ± 0.05	1.52 ± 0.07	0.40
CDg	1.483 ± 0.04	1.49 ± 0.05	0.12
CDm	1.38 ± 0.1	1.44 ± 0.08	0.68
μDconn	1.10 ± 0.04	1.14 ± 0.06	0.05

Table 1b

Results of t-tests for manually segmented Images

	Control (mean ± stdev) N = 11	PDR (mean ± stdev) N = 16	P value
DB	1.62 ± 0.04	1.61 ± 0.05	0.23
CDg	1.59	1.58	0.51
CDm	1.7	1.66	0.15
μDconn	1.19 ± 0.03	1.22 ± 0.04	0.03
μDconnnonskeletonized	1.40 ± 0.03	1.43 ± 0.05	0.02

PDR = proliferative retinopathy; D = box counting dimension; CD = global correlation dimension; CD = median correlation dimension; μD = mean local connected fractal dimension; μD = mean local connected fractal dimension for intermediate, nonskeletonised manually segmented images.

The distribution of the μD comparing control versus PDR for all of the automatically segmented images is shown in Figure 7. As shown in the tables above and in the figure, overall the μD for PDR was significantly higher compared to the control images.

Figure 7

The distribution of the μD per image for all automatically segmented images. The group of PDR images, which were more severely pathological, had a significantly higher μD (see Table 1a).

Comparison of pathological status using the distribution of the D

To further investigate the significant results reported above, we looked in more depth at local variation over individual images. Figure 8 below shows the distribution of individual D values per pixel for two images that were significantly different from each other in the μD (p = 0.02). Not shown here, we also investigated the distributions shown in the figure using color-coded maps that color each pixel according to its D. This analysis further suggests that D analysis is a suitable method able to differentiate PDR from other patterns in retinal vasculature in automatically segmented images.

Figure 8

Distributions of the D per pixel for two automatically segmented images. The μD values for the two patterns were significantly different (image 20__99_ = 1.10 ± 0.14 ; image 20__215 = 1.14 ± 0.15; p = 0.02); in addition, the distributions highlight specific areas where the D per pixel varied.

Discussion

Global increases in the age of the population and the prevalence of diabetes are bringing a concomitant increase in diabetic retinopathy. On the heels of that increase is a mounting need for accessible, early, and effective identification of proliferative changes and monitoring of disease progress and treatment outcomes. Practical and effective quantitative indices of DR are therefore essential, and may be found in fractal analysis. In general, our results showed that the D, an objective and readily obtainable measure, distinguishes between images of more severe and less severe pathology of the retina, where the μD is higher for pathological as compared to nonpathological retinal branching patterns in skeletonised automated images and both skeletonised and nonskeletonised manually drawn images.

Automated procedures should aim to be comparable to the results obtained by ophthalmologists. Using an automated method that can detect changes in the eye associated with diabetes with a minimum sensitivity of 60% is therefore a useful advancement as it would lessen the burden on ophthalmologists during initial population screening. (Daxer 1993a) For clinical use, however, the methods should exceed 80% accuracy. (BDA 1994; NHMRC 2001) Overall, our work suggests that CWT-based segmentation and D analysis can provide the required levels of sensitivity and accuracy.

These findings are encouraging, but further work (eg, using a larger sample) is required, and several points about our results need to be discussed, including the point that none of the other fractal analysis methods we assessed yielded significant results.

The global DF of the retinal vasculature has been studied for some time (Family et al 1989; Masters 1989). Several research groups have demonstrated that the normal retinal vasculature has fractal-like properties with a global DF generally falling between 1.60 and 1.88 (Daxer 1992; Landini et al 1993; Cesar and Jelinek 2003; Jelinek et al 2005). Our results suggest that the manually segmented images we analyzed were a valid representation inasmuch as the global DB was within this usually reported range (Table 1), but the same conclusion cannot be drawn for the automatically segmented images, because the mean global DB for that group fell below the range reported in the literature. Reasons we identified for this include that the entire sample was pushed to a consistently lower DB by an overall relative insensitivity to the finest level of branching in the segmenting method, and that certain images were further affected by factors such as uneven signaling and scarring.

The μD values were also lower than the reported range for the global DB, in this case for both the manually and automatically segmented sets of images. This does not reflect the limitations of the segmentation method described above, however. Rather, it highlights differences between global and local connected analyses. Whereas global dimensions sample entire images at once, the μD for an image is calculated from a sampling at each pixel of a pattern (ie, it is comparable to an average of averages). Moreover, each pixel’s D considers scaling only over a contiguous area (refer to Figure 2), and in our analysis, we further limited those areas to relatively small parts of the image in order to better capture local variation in scaling.

The first study that reported automated segmentation of blood vessels combined with fractal analysis used 30 mm SLR photographs of the retinal posterior pole that were scanned into a computer and segmented by the CWT. (Cesar and Jelinek 2003) Automated segmentation of fluorescein-labeled retinal vessels and analysis using the CD has also been reported in a study of digital images of the posterior pole with automated vessel segmentation followed by using 8 feature patterns including the CD, which was able to differentiate proliferative changes. (Jelinek et al 2005) However the CDs did not contribute significantly to our outcome. This is reflected again in Table 1.

The CWT is a powerful and versatile tool that has been applied in many different image processing problems, including shape analysis (Costa and Cesar Jr 2001). We found in pilot work for our study (results not shown) that selecting only images that convey what are considered the essential patterns associated with neovascularization and comparing these with control images that have no visible retinopathy improves the discerning ability of the D and indeed of any feature parameter. To be practical, however, DR screening methods will have to be able to address all the possibilities that will arise in retinal imaging scenarios. Thus, our test sample was deliberately not idealized; rather, we included images with various commonly encountered complicating factors.

As was discussed above, these factors, such as laser scars, vessel drop out, microaneurysms, and hemorrhage, probably did affect the absolute results for each image and may have been a major reason why segmentation using the CWT was unsuccessful except when followed by analysis with the D. Reasons the D but not other methods were successful in distinguishing PDR may also be linked to at least some of the reasons behind the range in values reported in the literature. That range reflects a host of additional methodological issues known to influence the absolute results of fractal analyses including the type of DF, image size and resolution, feature extraction methods such as variation in the use of skeletonisation (eg, unlike skeletonised images, nonskeletonised images preserve differences in vessel diameter in the final pattern), and whether red free or fluorescein images are used (Masters 2004).

In addition to such factors, previous studies using DFs have found various types of local differences within the retinal vascular tree. Differences have been reported between DFs for arteries and veins, for example (Mainster 1990; Landini et al 1993). Different investigators have also found contradictory trends in the DF associated with an increase in pathological status. Avakian and collaborators, for example, applied fractal analysis to region-based vascular changes in non-proliferative retinopathy (Avakian et al 2002). They found that the DF was significantly higher in the normal macular region compared to the NPDR macular region, although not elsewhere in the retina. Daxer also used fractal analysis of region-based vascular changes, but applied to proliferative retinopathy, using a method which requires that neovascularization be identified in the first instance. More in keeping with the results we found, Daxer reported that the DF was significantly higher for vessel patterns with neovascularization at the disk than without (Daxer 1993b).

Such variable and sometimes contradictory results can be reconciled using the D. To elaborate, if vessels disappear from an extracted pattern for any reason, such as from being occluded or by being obscured to the camera by hemorrhage, the local DFs calculated in that area of the pattern would be expected to decrease. Conversely, if vessels increase over some area, as seen in PDR, the local DFs in that area would be expected to increase. Global fractal dimensions smooth away such variety, but the D is likely to capture and quantitate the interplay without conflict because it quantifies complexity over the entire branching pattern but also locally within it, thus, does not lose important details inherent in that local variation. This was illustrated in Figure 8 shown earlier, which compared the distribution of the D over two automatically segmented, skeletonised patterns. The patterns shown in the figures were extracted from, respectively, one image of PDR with extensive capillary closure and one of normal retinal vasculature. The two images differ significantly in the frequency distribution of the D, which, as was discussed, can provide both a single objective index and various types of objective graphic representations (histograms and color-coded maps of the original retinal images) of the relative differences in the distributions of pathology between them.

Conclusion

Screening for and monitoring DR can decrease visual impairments and nationwide healthcare costs, and are important to expanding existing knowledge and improving treatments for DR. Automated methods and objective indices are vital components of this equation. Segmentation of vascular patterns using the CWT and analysis with the D shows promise in automating DR assessment, and appears robust to the realistic confounding factors that affect retinal images. With advances in digital imaging and the development of computerized grading systems as suggested by the work reported here, automated reading and assessment of diabetic eye disease, especially relevant to underserved populations, is coming closer to being a reality, but further work is required to verify the protocol introduced here and fully investigate the utility of the CWT and D in this regard.