# Advances in Image Processing Techniques for Drusens Detection and Quantification in Fundus Images

André Mora, Pedro Vieira, José Fonseca

► To cite this version:

André Mora, Pedro Vieira, José Fonseca. Advances in Image Processing Techniques for Drusens Detection and Quantification in Fundus Images. First IFIP WG 5.5/SOCOLNET Doctoral Conference on Computing, Electrical and Industrial Systems (DoCEIS), Feb 2010, Costa de Caparica, Portugal. pp.297-306, 10.1007/978-3-642-11628-5_32. hal-01060752

HAL Id: hal-01060752

https://inria.hal.science/hal-01060752

Submitted on 17 Nov 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Advances in Image Processing Techniques for Drusens
Detection and Quantification in Fundus Images

André Mora¹,², Pedro Vieira³ and José Fonseca¹,²

¹ Uninova, Campus da FCT-UNL, Monte da Caparica,
2829-516 Caparica, Portugal

DEE-FCT-Universidade Nova de Lisboa, Campus da FCT-UNL, Monte da Caparica,
2829-516 Caparica, Portugal

³ DF-FCT-UNL, Campus da FCT-UNL, Monte da Caparica,
2829-516 Caparica, Portugal

atm@uninova.pt; pmv@fct.unl.pt; jmf@uninova.pt

**Abstract.** Age-Related Macular Degeneration (ARMD) is considered the leading cause of irreversible blindness in developed countries. One of its risk factors is the presence of drusens, which are retina abnormalities appearing as yellowish spots in fundus images.

In this article a methodology using image processing techniques for the quantification of drusens is presented. The method uses splines combined with a contrast normalization to correct uneven illumination, followed by a drusen detection and modelling algorithm. The detection uses a gradient based segmentation algorithm that isolates drusens. They are then fitted by Gaussian functions, producing a model that is used to compute the area affected.

To validate the methodology, 22 images were marked by three ophthalmologists and compared to the automated method. The sensitivity and specificity for the automated process (0.664 and 0.963) were comparable to that obtained among the specialists (0.656 and 0.971). Also, the *Intraclass Correlation Coefficient* showed an agreement of 74.9% between the processed images and the specialists' analysis.

**Keywords:** image processing; retina imaging; automatic segmentation.

# 1 Introduction

Some ophthalmologic lesions have visible symptoms in the retina surface. Drusens, accumulations of extra-cellular materials beneath the retina surface, are one of those mentioned above. They are a risk factor for Age-Related Macular Degeneration, the leading cause of irreversible blindness in developed countries [1]. They appear in fundus images as speckles of higher intensity (see Fig. 1) and are qualitatively evaluated by the ophthalmologists by using fuzzy variables like: shape, size, location, etc. For this examination to be accurate, by quantifying the abnormalities, it requires a fastidious work from the expert that must be skilled for this task.

An accurate quantification of drusens in retina images using automatic procedures is one of the challenges for researchers who study medical imaging. This is a
<figure><img src="image_2.png" /><figcaption>Fig. 1. Examples of retina images obtained by fundus photography.</figcaption></figure>

subjective analysis and its criteria can differ among specialists and among analysis performed by one specialist at different moments. If we also consider that the design of an automatic diagnosis system is complex, then the development of such accurate procedure is a difficult task.

## 2 Contribution to Technological Innovation

Several studies have been published in the last twenty-years for drusens detection using image enhancement techniques. They included image pre-processing for illumination compensation and contrast improvement, and for the drusens segmentation have been proposed algorithms using adaptive local thresholds [2-7], global thresholds [8] or fuzzy logic thresholds [9]. Most of these techniques proved that good noise removal algorithm is required for getting reliable results. We have also contributed to this subject by introducing new image processing algorithms for illumination correction [10] and for drusens detection and modelling [11].

In this article it is proposed a methodology for automated drusens detection and quantification. The main goal is to create an algorithm which is as non-parametric as possible. Parameters should be calculated based on the image information, requiring less user intervention. The proposed methodology includes an image pre-processing step to compensate non-uniform illumination and contrast normalization; a detection algorithm that determines the number of drusens and an image modelling to quantify the area affected. The methodology was validated on a set of 22 images with the collaboration of three ophthalmologists.

## 3 Illumination Correction

In ophthalmologic imaging it is usually difficult to obtain good quality images. Illumination is not always uniform due to numerous factors: the retina surface is convex; the macula reflectance is lower than others parts of the retina; the patient eye is not stable, etc. The consequence is that the image has different contrast areas
<figure><img src="image_3.png" /><figcaption>Fig. 2. Illumination correction: (a) Original image f(x,y); (b) Estimated illumination model i(x,y); (c) Corrected image.</figcaption></figure>

increasing the complexity of the quantification methods. We propose to correct the non-uniformity of illumination and contrast before further processing.

The proposed technique is designed for gray scale images and since the retina images used are mostly in RGB format, a RGB to Gray conversion is applied to those. This conversion is obtained using the green channel information, which is a procedure widely used in other retina image works [4, 7, 12]. It is chosen by its better dynamic range, better contrast and less influence by the illumination differences, while it maintains all retina structures visible.

Using the gray scale image, and to correct the effects of non-uniform illumination, we then determine the image background intensity and correct the image accordingly. Several methods were tested (FFT, Gaussian fitting, Gaussian blurring, Splines) [10] and the smoothing spline fitting was the one to produce better contrasted images. This method fits the image background f(x,y) with a smoothing cubic spline (Fig. 2.b). This is a cubic spline S(x,y) that minimizes the function presented in equation (1), where p is a smoothing factor and D²S(x, y) is the second derivative of the cubic spline:

$$p \cdot \sum_{y=1}^{m} \sum_{x=1}^{n} |f(x, y) - S(x, y)|^2 + (1-p) \iint |D^2 S(x, y)|^2 dxdy \quad (1)$$

The smoothing factor p, defined in the interval [0..1], controls the relative weight between fitting all data points (p=1) and having a smooth spline (p=0). For p=1 the spline is an interpolating spline that crosses every data point, while p=0 it maximizes the smoothness generating a straight line. The reference value for p was determined empirically within our set of test images and obtained p = 1e-6, for an image resolution of 12.5 µm/pixel.

The corrected image is obtained by dividing the original image by the illumination pattern, generating a unitary average intensity. To recreate a grey scale image with 256 levels a constant value, 85 (1/3 scale), is multiplied by the image. However, depending on the image contrast this value can generate out of scale values or low contrast images. To overcome this problem a contrast normalization procedure was introduced. It involves the determination of the standard deviation (σ) of the image intensity and its normalization to a predefined value.
The retinal vessels are one of the structures that are present in all images and provide good contrast with the background. We propose a method that evaluates locally the contrast between the vessels and the background and normalizes it to a predefined value.

Using a sliding window-based analysis the σ value and the mean intensity are calculated for each window. A window width of 1/8 of the optic disc diameter was considered to be adequate, as it can contain small vessel branches and background. Then from the 50 darker windows, that are typically windows containing vessels over a uniform background, is determined the median of the σ values (σ<sub>median</sub>). The image intensity is finally updated using the equation (2), where σ<sub>predefined</sub> is the predefined σ value (0.15) and A the constant background value (85).

$$norm\_image(x, y) = A \cdot \left( \frac{\sigma_{predefined}}{\sigma_{median}} * (image(x, y) - 1) + 1 \right) \quad (2)$$

## 4 Drusens Detection

Our Drusens detection and quantification algorithm is based on the modelling of Drusen spots to ensure a shape consistent image segmentation. For getting a better performance from the modelling algorithm it is recommended that an image pre-processing is applied to calculate the approximate number of Drusens, their location and optionally their size.

We propose a maximum detection algorithm based on labelling the maximum gradient path. In medium resolution images there are more then just one pixel pointing to each *intensity maximum*. Therefore, it is possible to follow one or more ascending paths that inevitably reach an *intensity maximum*.

The algorithm begins by determining the image gradient using a 3x3 Sobel operator which evaluates intensity changes in the horizontal and vertical axis. The pixel gradient is a vector which is oriented to the ascending pixel. The two-

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Example of Drusens center detection algorithm. (a) Original Image and gradient; (b) Labelling 1st stage – initial label propagation – (with label propagation on the two upper rows); (c) Labelling 2nd stage – apply compatibilities - (two Drusens centered on the highlighted pixels).</figcaption></figure>
dimensional gradient analyses the 3x3 neighbourhood of each pixel what eliminates
several noise related artefacts. But when signal/noise ratio is low the algorithm can be
significantly improved by changing every pixel gradient by the mean (vectors sum) of
its neighbour's gradients. Fig. 3.a contains an example of an image after applying a
3x3 gradient mean. The original image represents two small Drusens.

In the first stage of the labelling procedure, *initial label propagation*, every pixel is examined in a top-left to bottom-right direction and is assigned a new label to all pixels that are not marked yet. Whenever a label is assigned to a pixel it is propagated to the neighbour pixel that is in the gradient azimuth direction. This labelling continues on the propagated pixels until the next pixel is already marked. Fig. 3.b presents the result of the labelling first stage. Every time the label propagation finishes on the same label that is being propagated, that pixel is marked as an *intensity maximum*. In the same situation when there is a different label the two labels are defined as being *compatibles*. This means that they belong to the same *intensity maximum*, i.e., the same Drusen spot.

The second stage of the labelling procedure is to apply the *label compatibilities*
resulting on an image with as many labels as possible Drusens. The result is a
segmented image where all pixels that contribute to a spot have the same label. It is
important to notice that the algorithm is so far a non-parameterized algorithm.

But a problem arises when flat valleys or flat hills occur. In these cases not all
gradient paths end on the same *intensity maximum* pixel since there might exist more
then one pixel for that *intensity maximum*. Consequently, the algorithm will generate
more *intensity maximums* then it really should. For solving this problem a merging
strategy was followed. Using a connected components approach, the local *intensity
maximums* are merged to other neighbour *intensity maximums* if they can be
connected by a path that doesn't goes lower than a predefined amplitude.

At this stage drusens are detected and characterized by centre coordinates and size.
This information will be used as the initial parameters for the modelling algorithm.

5 Drusens Modelling

The final phase is to analytically characterize the shape and intensity of the
drusens, providing information to evaluate them quantitatively. The elevations shown
on drusen areas on the tri-dimensional representation motivated us to create a model
of the image intensity and from it extract the drusens dimensions. From the model it is
computed the total affected area and drusens dimensions. This information can be
valuable to an alongside comparison that evaluates sequence of images during a long
term treatment.

In a three-dimensional view of the retina surface, as in Fig. 4.a, it can be seen a
drusen that has a waveform similar to a Gaussian function superimposed to a flat
background with low amplitude white noise. These observations and the possibility to
describe the drusens by means of elementary functions motivated the use of Gaussian
functions to model drusens. This approach allows the drusen modelling and its
isolation from the background.
The modified Gaussian function used is given by:

$$G(x, y) = A * e^{-\left(\frac{\frac{x^2}{s_x} + \frac{Y^2}{s_y}}{\frac{2}{d}}\right)} + Z_0. \quad (3)$$

where

$$\begin{align}
X &= (x - x_0) \cdot \cos(\theta) + (y - y_0) \cdot \sin(\theta) \nonumber \\
Y &= -(x - x_0) \cdot \sin(\theta) + (y - y_0) \cdot \cos(\theta) \tag{4} \\
s_y &= s_F \cdot s_x \nonumber
\end{align}$$

The Gaussian function used represented on equation (3) allows: translations in the xy plane ($x_0, y_0$) and in the background plane ($Z_0$), amplitude scaling ($A$), rotation ($\theta$) and shape modifications ($s_x, s_F, d$). These latter defines the width in the x-plane ($s_x$), the shape factor ($s_F$) used to obtain the width in y-plane, allowing stretched shapes in the xy plane, and the shape factor ($d$) that controls the amplitude profile between square shape, round shape and thin shape.

For searching the function parameters that best fits the image we adopted the Levenberg-Marquardt non-linear least squares optimization algorithm with pre-defined parameters constraints ([13]). Because some Drusens are overlapped, a multiple function (one or more for each spot) simultaneous optimization must be used for studying the sum of all the Gaussian functions. In Fig. 4.b there is an example of superimposed Drusens *modeled* through multiple Gaussian functions.

Finally, for highlighting Drusen areas in the image it is shown the contour or the area of each detected spot at predefined percentage of its amplitude (without background).

<figure><img src="image_5.png" /><figcaption>Fig. 4 (a) 3D original image view; (b) Drusens modeling using multiple Gaussian functions.</figcaption></figure>
# 6 Results

For evaluation of the drusen detection and quantification process a set of 22 images containing colour and greyscale fundus images was marked by three independent ophthalmologists. Our study was focused on the macula, as suggested by the Wisconsin Grading System [14] a circular *region of interest* centred in the macula with 3000µm diameter. The specialists marked this set of images using the computer assisted procedure with the software *MD3RI* [15] which were then used to compare the results produced by the proposed method (*AD3RI*) over the same images.

For the results analysis the percentage of drusen area, the *false positive* and *false negative* pixels, were computed. The results were then evaluated using the algorithm's *sensitivity* and *specificity* as well as using the *Intraclass Correlation Coefficient* (ICC) to evaluate the agreement between all the specialists and the *AD3RI* method and among each specialist.

The *sensitivity* and *specificity* indicators were obtained by comparing the processed image with each specialist marking. The average *sensitivity* and *specificity* were 0.664 and 0.963, respectively. The same indicators were evaluated among the specialists as a comparative indicator and the average *sensitivity* and *specificity* were comparable to the previous ones, 0.656 and 0.971, respectively.

As a comparative analysis method, the *ICC* with 95% confidence levels [16] evaluated the areas' values between specialists and the *AD3RI* (*ICC<sub>AD3RI</sub>*).

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Two examples of drusens markings.</figcaption></figure>
specialists (*ICC*<sub>Spec</sub>). These indicators were all within the 95% confidence interval and obtained an average agreement of 74.9% between the specialists and the *AD3RI* and an average agreement of 81.1% among the specialists.

# 7 Discussion

The development of methods to quantitatively measure drusens in a reproducible and accurate procedure will improve the quality of the follow up of this disease and potentiate epidemiologic studies and clinical trials. In this article we presented a new method to quantitatively measure drusens and its' comparison with 22 images marked by three independent ophthalmologists. The algorithm is based on the detection and modelling of drusens to automatically calculate the affected areas including an image pre-processing for non-uniform illumination correction.

The presented non-uniform illumination correction algorithm was an important step to obtain a less parameterized methodology, since it is capable to deliver an image with normalized illumination and contrast to the following steps. The detection and modelling of drusens with Gaussian functions demonstrated to be capable to detect drusens even with low contrast and to segment the image maintaining the drusens typical shape. It also provides an analytical model which allows the determination of drusen indicators, such as confluence, average size, type of drusens and others.

From the statistical analysis of the algorithm with three specialists related to the 22 images was observed that the algorithm performs as do specialists. On the one hand, the sensitivity (0.664) was slightly higher than the one observed among the specialists (0.656). This is due to the algorithm detection of a higher overall area, increasing the probability of agreement with the specialists. On the other hand, the high specificity value (0.963) demonstrates that the automated process is very accurate and is also comparable to results obtained among the specialists (0.971). With the algorithm, the value is lower due to the increase of false positives, which results of the detection of more drusens. We consider that the lower quantification of drusens by the specialists is mainly due to the marking process being fastidious.

The other measure of agreement used, the *Intraclass Correlation Coefficient (ICC)*, demonstrated that the agreement among specialists on the overall areas values is approximately 81% and that the automatically processed images have 74.9% of agreement when compared to specialists. This analysis showed a good concordance between the automated process and the specialists, verifying that the analysis is statistically significant, as all the obtained *ICC* values were within the 95% confidence levels.

From the above, the *AD3RI* method showed promising results to be tested by a wider panel of specialists and with new sets of images. It detected the most important drusen as did the specialists, having always the same level of detail independently of the specialists' attention and accuracy. Also, the qualitative feedback from the specialists confirmed that the results were accurate and that the software can play an important role in the follow-up of retina diseases.
**Acknowledgments.** This work was financed by the Fundação para a Ciência e Tecnologia, through the POCTI and POSI Research Programs (project nº POCI/SAU-ESP/57592/2004). The authors acknowledge the University of Aberdeen, Rudolfstiftung Hospital, Hospital Santa Maria and Faculdade de Ciências Médicas for supplying the retina images used in this work, for their valuable feedback on the software testing and for marking the sample images.

## References

1. World Health Organization, International Agency for the Prevention of Blindness: State of World's Sight - Vision 2020: the right to sight (1999-2005) (2005)

2. Peli, E., Lahav, M.: Drusen Measurement from Fundus Photographs Using Computer Image Analysis. Ophtalmology, vol. 93, pp. 1575--1580 (1986)

3. Kirkpatrick, J.N.P., Spencer, T., Manivannan, A., Sharp, P.F., Forrester, J.V.: Quantitative image analysis of macular drusen from fundus photographs and scanning laser ophthalmoscope images. Eye (Royal College of Ophthalmologists) 9, pp. 48--55 (1995)

4. Rapantzikos, K., Zervakis, M., Balas, K.: Detection and segmentation of drusen deposits on human retina: potential in the diagnosis of age-related macular degeneration. Med Image Anal 7, pp. 95--108 (2003)

5. Morgan, W.H., Cooper, R.L., Constable, I.J., Eikelboom, R.H.: Automated extraction and quantification of macular drusen from fundal photographs. Australian and New Zealand Journal of Ophthalmology 22, pp. 7--12 (1994)

6. Phillips, R.P., Spencer, T., Ross, P.G., Sharp, P.F., Forrester, J.V.: Quantification of diabetic maculopathy by digital imaging of the fundus. Eye 5 (Pt 1), pp. 130--137 (1991)

7. Shin, D., Javornik, N., Berger, J.: Computer-assisted, interactive fundus image processing for macular drusen quantitation. Ophthalmology 106, pp. 1119--1125 (1999)

8. Smith, R.T., Nagasaki, T., Sparrow, J.R., Barbazetto, I., Klaver, C.C., Chan, J.K.: A method of drusen measurement based on the geometry of fundus reflectance. Biomed Eng Online 2 10 (2003)

9. Thdibaoui, A., Rajn, A., Bunel, P.: A fuzzy logic approach to drusen detection in retinal angiographic images. 15th International Conference on Pattern Recognition, Vol. 4, Barcelona, Spain, pp. 748--751 (2000)

10. Mora, A., Fonseca, J., Vieira, P.: Drusen Deposits Modelling with Illumination Correction. Biomed-2005, Innsbruck, Austria (2005)

11. Mora, A., Vieira, P., Fonseca, J.: Modeling of Drusen Deposits Based on Retina Image Tridimensional Information. Second International Conference on Computacional Intelligence in Medicine and Healthcare - CIMED-2005, Costa da Caparica, Portugal (2005)

12. Soliz, P., Wilson, M.P., Nemeth, S.C., Nguyen, P.: Computer-aided methods for quantitative assessment of longitudinal changes in retinal images presenting with maculopathy. Medical Imaging 2002: Visualization, Image-Guided Procedures, and Display, Vol. 4681. SPIE, San Diego, CA, USA, pp. 159--170 (2002)

13. Marquardt, D.W.: An algorithm for least-squares estimation of non-linear parameters. Journal of the Society for Industrial and Applied Mathematics 11, pp. 431--441 (1963)

14. Klein, R., Davis, M.D., Magli, Y.L., Segal, P., Klein, B.E., Hubbard, L.: The Wisconsin age-related maculopathy grading system. Ophthalmology 98, pp. 1128--1134 (1991)

15. Mora, A., Vieira, P., Fonseca, J.: MD3RI a Tool for Computer-Aided Drusens Contour Drawing. Fourth IASTED International Conference on Biomedical Engineering - BIOMED2006. ACTA Press, Innsbruck, Austria (2006)
16. McGraw, K.O., Wong, S.P.: Forming inferences about some intraclass correlation coefficients. Psychological Methods **1**, pp. 30--46 (1996)