Quantification of Parapapillary Atrophy and Optic Disc

Abstract

Purpose.: A computer-aided measuring tool was devised to automatically detect and quantify both the parapapillary atrophy (PPA) and the optic disc (OD) regions in two-dimensional color fundus images of the retina.

Methods.: The OD region was segmented using the Chan-Vese model with a shape restraint. This region was then removed from the image (OD+PPA), which was cropped in a modified Chan-Vese approach, producing a first-order estimation of the PPA region. Its boundary is subsequently refined by using thresholding, a scanning filter, and multiseed region-growing methods. Dual channels (blue and red) in the red-green-blue space are used to minimize the interference effects of blood vessels and artifacts.

Results.: The software was tested on 94 randomly selected images of eyes with PPA from 66 subjects of a well-characterized cohort database. Our proposed algorithm achieved a mean accuracy level of 93.8% (SD 5.26) and 94.0% (SD 5.88) in estimating the size of the PPA and OD respectively, compared with the ground estimate defined by an ophthalmologist. In terms of correlation between the data of the ground estimate and our estimation, we obtained a correlation coefficient of 0.98 for both the PPA and the OD.

Conclusions.: This software offers a means of quantifying the size of PPA on two-dimensional fundus images for the first time. The proposed algorithm is capable of detecting and quantifying PPA and OD regions repeatedly, with a mean accuracy of >93%, and could also provide additional information, such as the transverse and conjugate diameter of OD, which may be useful in eye-screening.

Certain ophthalmic diseases (e.g., progressive glaucoma) and eye conditions (e.g., myopia) have been associated with the development of retinal pigment epithelium parapapillary atrophy (PPA). ¹ Although thinning and degeneration in retinal tissue are generally irreversible, early detection and medical intervention may offer an opportunity to reduce or limit their progression. ² More recently, there has been a growing interest in the potential role of PPA in detecting eye diseases at an early stage. ^{3 –9}  

Ophthalmoscopically, the PPA is divided into the β zone, located adjacent to the optic disc (OD) border and showing visible sclera as well as large choroidal vessels, and the α zone, located more peripherally and characterized by irregular hyper- and hypopigmentation. This classification of β- and α-PPA was first developed by Jonas et al. ¹⁰ and has now been accepted into common use, providing the motivation for studies to better understand the development of the PPA. For instance, Healey et al. ³ investigated the association between β-PPA and both environment and genes. They found that the presence of β-PPA appeared to be under strong genetic control, with only a fraction of this genetic effect shared with genes involved in myopia. Tezel et al. ⁴ studied the relationship between PPA and the OD in patients with glaucoma or visual field damage, by manually quantifying the size of the PPA and OD regions. The presence and the progression of PPA were found to be associated with subsequent OD or visual field changes in patients with ocular hypertension. Honrubia and Calonge ⁵ further found that, in patients with ocular hypertension, the presence of PPA could imply a risk of glaucomatous deterioration, as it has a significant association with the presence of defects in the retinal nerve fiber layer (RNFL). Xu et al. ⁶ investigated the association of PPA with glaucoma in a population-based study. They found that β-PPA appeared to be larger and occurred more frequently in glaucomatous eyes than in normal eyes of Chinese adults, but no significant difference was found between chronic open-angle and primary angle-closure glaucoma. ⁶  

Uhm et al. ⁷ used a commercial software package (Photoshop; Adobe Systems Inc., San Jose, CA) to measure the size of PPA and the OD regions in assessing how closely PPA is related to functional and structural glaucomatous optic nerve damage. They discovered that the severity of glaucomatous optic nerve damage and visual field defects were related to the increases in the size of PPA and concluded that PPA could be useful for the diagnosis and progression of glaucomatous nerve damage. ⁷  

There are ophthalmic techniques, such as retinal tomography (HRT; Heidelberg Engineering, Heidelberg, Germany) and optical coherence tomography (OCT), that can provide a colorless or pseudocolor 3-D visualization of the PPA and the OD. ¹¹ They require a trained technician or ophthalmologist to manually identify the OD boundary on the image before the PPA and OD variables can be estimated from the image contour based on three-dimensional depth information. ⁸ Furthermore, several automatic segmentation schemes ^{12 –23} have been developed for the assessment of the OD and parapapillary features, such as RNFL thickness. 

The HRT and the OCT have been used in several studies about the PPA. One such application was by Uchida et al., ¹³ who investigated the association between PPA and the progression of glaucoma. However, two-dimensional color fundus imaging may be preferred by many clinicians for estimating the size of the OD and PPA, because the image from HRT or OCT is generated via computer processing rather than a direct recording of the object of interest. 

An alternative tool to detect and quantify the PPA and the OD automatically from color fundus images would reduce the workload of the human assessor and could facilitate a wider investigation about the potential importance of PPA in ophthalmic disease diagnosis involving facilities where access to HRT or OCT is limited. Moreover, in large clinical studies and pathologic monitoring programs, such a computerized tool can also improve repeatability, therefore avoiding problems associated with fatigue and habituation, and tends to be more cost effective. 

A novel approach to automatically segment and quantify the OD and PPA region from color fundus images is proposed in this article. The methodology exploits both the red and blue channels of the color fundus image to maximize extraction of information on the features (PPA), while keeping interference (blood vessels) to a minimum. A combination of several techniques, including scanning filter, autothresholding, and region-growing, as well as the Chan-Vese (C-V) model with a shape restraint is used to segment and quantify the regions of the OD and PPA. 

Methods

We used color retinal images drawn from a subsample of the Lothian Birth Cohort (LBC), a 1936 study. ²⁴ The participants comprise surviving members of the Scottish Mental Survey of 1947 (n = 70,805) who were born in 1936 and currently reside in the Edinburgh area (Lothian) of Scotland. Eight hundred sixty-six of them were successfully traced and participated in a series of investigations that included having retinal photographs taken at the Wellcome Trust Clinical Research Facility, Western General Hospital, NHS Lothian, Scotland. Their mean age at the time of the photography was 72 years. The research complied with the Declaration of Helsinki and was approved by the Lothian (Scotland A) Research Ethics Committee. 

The computer program in this work was implemented in commercial software (MatLab; MathWorks Inc., Natick, MA). The images were evaluated by an ophthalmologist (AL) who was masked to the image processing findings. The intensity variation between the bright objects (i.e., the OD and the PPA) and the blood vessels of the retinal image was relatively high (Fig 1). Conversely, the blood vessels were in general at a lower intensity level with respect to the background. 

Figure 1.

 

Original fundus image with the OD and PPA.

Our proposed algorithm combined a collection of image-processing techniques (Figure 2). Fundus images were initially preprocessed in two channels of the RGB space to reduce the interference of blood vessels and to better distinguish the regions of OD and PPA. The OD region could be reliably detected in the red channel, as it appeared brighter than the rest of the image, while the blood vessels appeared least influential. ²⁰ We also found that the region consisting of both the OD and the PPA (hereafter referred to as the region of OD-plus-PPA) was most well-defined in the blue channel. Consequently, the region of OD-plus-PPA was first segmented by a modified C-V model in the blue channel. Then, a variant of the C-V model with a shape restraint was applied to segment the OD region in the red channel. In this case, the restraint was based on an ellipse reflecting the actual shape of an OD. Removing the OD region from the region of OD-plus-PPA produced the first-order estimation of the PPA region. Moving back to the blue channel, the segmented image was then equally divided into four zones automatically. 

Figure 2.

 

Flow chart shows the extraction of the PPA and the OD regions.

Based on the autoset thresholds acquired from each zone, the image was then filtered to reduce the influence of crossing vessels and artifacts. Finally, the PPA was extracted by using a multiseed region-growing method. ²⁵  

Chan-Vese Model

The C-V model is a type of active contour model or “snake” that can trace the outline of an object from an image by minimizing an energy function associated with the current object contour. ²⁶ It combines methods including curve evolution, Mumford-Shah function, ²⁷ and level sets for applications such as shape recognition, edge detection, and image segmentation. In this work, we used it to identify topological changes, corners, and cusps associated with the presence of PPA. The C-V model has been enhanced to detect objects whose boundaries are not all necessarily well-defined. The step size of the energy function was selected carefully, to ensure that the snake would stop at the desired boundary. More details on the C-V model are given in Supplementary Material S1. 

Segmentation with the Modified C-V Model

As the OD region appears to be more or less an ellipse, we adopted the model of Tang et al. ²⁸ which integrated the C-V model with an elliptical shape restraint. More details on the model are given in Supplementary Material S2 We also used it to segment the OD-plus-PPA region. However, the PPA region may sometimes appear in an irregular shape, and so we had to modify the rules slightly for the evolution of the C-V model. The set of starting C-V model points, also known as the initial mask, was arranged to be in an ellipse as per normal. However, the model was then allowed to deform freely as it edged closer to the boundary of OD-plus-PPA with each subsequent iteration. This method allowed the model to produce an enclosed but not necessarily elliptical shape, which was always bigger than the exact region of the OD-plus-PPA. 

Next, Tang's C-V model was employed to detect the OD region. To accurately segment the OD, we introduced two modifications to the model. First, Tang's equation B-3²⁸ was restored to its original form (of an ellipse) :  

 

Second, we introduced a new way to automatically detect the center of the OD for more accurate segmentation. We divided the image into four subregions. Then, we adopted the approach used in Tang's model to estimate the initial mask center (x₀, y₀). The initial function in the equation B-2²⁸ was chosen as:  

 

where R is the estimated radius of the OD and can be defined as  

 

where w and h are the width and height of the image, respectively. 

Then, our algorithm calculated automatically the offset, f_x and f _y, of the initial mask center, based on the histogram of intensity value of each four regions. The updated initial mask center is thus  

 

 

Autosetting of Thresholds and Scanning Filter

To eliminate the unwanted pixels in the oversized OD-plus-PPA region, we acquired threshold values from the histogram of intensity values in the four subregions. In this context, the threshold was set by the brightest 30% of all pixels in each region, thus producing a better-defined OD-plus-PPA region. We then subtracted the OD region and obtained the first-order estimation of the PPA region as illustrated in Figure 2f. 

The PPA region is a nonhomogeneous area divided into multiple sections by a few crossing blood vessels. We therefore proposed the use of a scanning 1 × 3 filter to create a path through the vessels for the following region-growing model to reach different sections of the PPA. 

Multi-initial Seed Region-Growing

The seeded region-growing technique was first introduced by Adams and Bischof. ²⁵ It starts with several initial seeds and then adds adjacent points as new seeds, beginning with the points of lowest priority. The priority is defined by a distance function. The distance of each pixel to a contiguous region ²² is defined by:  

 

where I(x) is the gray image value of the point x(element) δ and δ _i , the region labeled i. Above all, setting both the right initial seed and distance function is the most important step in PPA extraction. 

Our algorithm automatically placed one initial seed in each of the four subregions and set an optimal distance function for each subregion. Each seed was then allowed to grow until the regional threshold distance set by equation 5 was met. Finally, we combined the results at all four subregions to produce an integrated PPA region. By combining the techniques listed above, our methodology permitted the full use of both global and local information for PPA and OD segmentation. 

Results

A total of 94 color fundus images (including 18 poor-quality images, as determined by an independent and experienced human assessor, AL) from 66 subjects were randomly selected from the Lothian Birth Cohort (LBC) database. Without knowing the segmentation results from the proposed tool, the human assessor provided the ground estimate of the OD and the PPA region in the images. Subsequently, the area enclosed by the ground estimate was counted pixel by pixel (Photoshop; Adobe Systems Inc.) to quantify the size of each region. The count was repeated with the segmentation results from the tool. 

Figure 3 shows six samples from the segmentation results of the proposed tool. The first column depicts the results obtained from good-quality images, and the second column depicts the results from poor-quality images. The ground estimate is enclosed by the black solid line. The results of estimated PPA and OD region are enclosed by the spots and red markings, respectively. Figure 4a shows the comparison of the OD area (in arbitrary pixel unit), based on the ground estimate, and the estimated OD area, determined by the proposed tool in the 94 trials, along with the line of best fit. Figure 4b shows a similar graph but for the PPA area size estimation. 

Figure 3.

 

Segmentation results from the proposed algorithm. Shown are (a, c, e) good-quality and (b, d, f) poor-quality images. Black solid line: the ground estimate; black dots and the red markings: the estimated PPA and OD regions, respectively.

Figure 4.

 

The correlation between the ground estimate (x-axis) and the results obtained by the proposed tool (y-axis) in quantifying the size of each region, in arbitrary pixel units. (a, b) Direct estimation results from the tool; (c, d) estimation results of the OD and PPA region after calibration, such that y = x. The correlation coefficient was 0.98 in all cases.

The results suggest that the proposed algorithm or estimation model was able to detect the general boundary of OD and PPA. However, it tended to terminate the snake evolution prematurely on all good-quality images; hence, the results appeared to underestimate the actual size. This underestimation was consistent in the case of poor-quality images, in which the intensity variation/resolution in defining the boundary of OD was limited. Figure 3d shows a good example of when the model missed the mark by pushing the boundary into the scleral rim. Overall, it appears that most of the results in the 94 trials are underestimations. This is confirmed by the gradients (both <1) of the best-fit lines in Figures 4a and 4b. We therefore calibrated our estimation model by using these values as scaling factors. The final estimation results are plotted in Figures 4c and 4d. As shown, a correlation coefficient of 0.98 (max = 1) was achieved in the size estimation of both the OD and PPA regions. This result suggests that our estimation is not stochastic but is fairly consistent with the ground estimate defined by an ophthalmologist (i.e., the best-fit line is defined by the equation y = x). 

Validity of the Tool

In this section, we use three methods to validate our estimation model. First, we calculated the mean accuracy (MA) of our estimation model, which is given by:  

 

where S _a represents the actual size (ground estimate) of PPA or OD, while S _e represents the estimated size (by our model) of PPA or OD. The numerical value n is the total number of images analyzed in our experiment. Before calibration, our estimation model achieved a mean (SD) accuracy level of 90.0% (6.20%) and 88.2% (5.85%), in defining the size of OD and PPA, respectively. With the same set of color fundus images, our estimation model after calibration achieved a mean (SD) accuracy level of 94.0% (5.88%) and 93.8% (5.26%). Table 1 summarizes the estimation results. 

Table 1.

 

The Results of PPA and OD Segmentation in 94 Trials

Results	Before Calibration		After Calibration
	PPA	OD	PPA	OD
Mean accuracy, %	88.2	90.0	93.8	94.0
Standard deviation	5.85	6.20	5.26	5.88
Correlation coefficient, R 2	0.98	0.98	0.98	0.98

Second, we calculated the sensitivity of our estimation model at different tolerance levels. In this context, tolerance level refers to the percentage of estimation error being acceptable. When the estimation result by our model falls within the range of tolerance, we count it as one correct prediction, and vice versa. The sensitivity is defined as the percentage of correction prediction over the total number of images analyzed in our experiment. Our model after calibration achieved a sensitivity of 75.5% (tolerance, ±10%) and 90.4% (tolerance, ±15%) in estimating the size of OD region. As expected, the sensitivity level is lower when the tolerance is smaller and when our model is not calibrated (Fig. 5a). The scatterplots in Figure 5b show a similar trend. After calibration, the sensitivity of our model for the PPA region was 84.0% (tolerance, ±10%) and 92.6% (tolerance, ±15%), respectively. 

Figure 5.

 

The sensitivity of the proposed tool in defining the size of (a) the OD and (b) the PPA regions at different tolerance levels in the 94 trials. The slopes represent the best fit for the scatterplots.

Third and finally, we evaluated the robustness of the calibration method for our estimation model. We randomly selected 70 of the 94 images and derived the suitable scaling factor to calibrate the segmentation results, according to the aforementioned methodology. Then, we used the remaining 24 images to obtain the accuracy of the calibrated model. In this case (labeled RA), the mean (SD) accuracy was 94.0% (5.77%) and 94.3% (4.53%) in defining the size of OD and PPA, respectively, in the 24 test images. We repeated this experiment by (1) calibrating our model with images that had the largest (70/94; 74.5%) PPA and had remaining images from the dataset—that is, those with the smallest (24/94; 21.3%) PPA as test images, and (2) calibrating with the smallest (74.5%) PPA and testing with the largest (21.3%) PPA. The earlier (labeled HI) achieved a mean (SD) accuracy of 92.7% (7.53%) and 92.8% (4.93%) for size estimation of the OD and PPA regions, respectively. On the other hand, the latter (labeled LO) achieved a mean (SD) accuracy of 93.7% (6.98%) and 93.4% (6.57%) for the OD and PPA. Table 2 summarizes the results, including the previous results from calibrating with all 94 images (labeled AL). As expected, the model performed the worst when it was calibrated with the largest PPA images and tested with the smallest PPA images. However, the mean accuracy of our model is still greater than 92.7%. 

Table 2.

 

The Statistical Data of PPA and OD Segmentation Results in Different Combinations of Calibration and Test Images

	OD				PPA
	HI	RA	LO	AL	Hi	RA	LO	AL
Scaling factor	0.92	0.91	0.90	0.92	0.88	0.87	0.89	0.88
Correlation coefficient	0.97	0.98	0.99	0.98	0.97	0.97	0.99	0.98
Mean accuracy, %	92.7	94.0	93.7	94.0	92.8	94.3	93.4	93.8
Standard deviation	7.53	5.77	6.98	5.88	4.93	4.53	6.57	5.26

 All images were taken from the same set of LBC images.

Discussion

We have developed automated software to measure the size of PPA and the OD in two-dimensional fundus images. Our experimental results with a very wide variety of fundus images showed that the proposed algorithm was not only robust for automatic PPA shape detection and area quantification, but it could also provide the transverse and conjugate diameter of the OD as well as PPA-to-OD ratio, which may be useful in early detection and grading of eye conditions such as glaucoma. 

The automatic detection and the quantification of OD in fundus images are particularly important tasks in retinal image analysis for two reasons: First, the OD has attributes similar to the PPA, both in terms of contrast and brightness, making their boundary detection a difficult task. Second, the OD is often seen as a landmark that can be used for a coarse localization of the area of interest in retinal images, reducing the search area during the preprocessing stage. 

Several schemes for automatic segmentation of the OD have been reported. ^{11,14 –19,21 –23} The OD can be detected either by finding a large cluster of pixels with high intensity ^22,23 or by the highest intensity variation at the gray level ^11,21 ; however, difficulties are often encountered when exudates coexist within the retinal image, because exudates also have a higher intensity level than the surrounding regions of the OD. Walter et al. ¹⁹ developed a combined approach of watershed transformation and morphologic filtering to detect OD, but found that morphologic filtering could not completely remove the distortion caused by overlying retinal vessels. Another approach used an area-thresholding algorithm to localize the OD, ¹⁸ before detecting its boundary by means of a Hough transform (HT) (i.e., best fitting circle based on the gradient information of the image). However, this approach proved to be time consuming and relied on certain forms of the OD that were not always encountered. Principal component analysis (PCA) for automatic detection of the OD has been reported ¹⁷ and could be used, even in the presence of bright lesions on the fundus image, although this approach could also be time-consuming. Alternatively, Osareh et al. ¹⁶ used a template-matching algorithm to detect the disc boundary automatically. Although morphologic preprocessing helped to reduce the interference effects of blood vessels, it could not remove them completely. Moreover, such processing blurred the OD boundary, making the detection unreliable. The C-V method ¹⁵ and level setting methods ¹⁴ have also been applied to OD boundary segmentation. The major advantage of these algorithms is their ability to compensate discontinuities in the boundary of the image feature to be located. However, those approaches have to be carefully initialized and can achieve only good segmentation results when the region has homogenous intensity values and a well-defined boundary. 

In this work, we explored a dual-channel approach with a modified C-V model to segment the PPA and the OD individually. The proposed algorithm particularly is designed to address the aforementioned challenges on how to maximize information extraction of features (OD/PPA) while keeping interference (blood vessels) to a minimum. 

In comparison to the works by other groups, there are three main merits of the proposed tool. First, our software tool could measure the PPA region automatically in two-dimensional color fundus images. This tool is the first tool that can make that measurement without human intervention. In previous studies, investigators had to manually measure the PPA region in either 2-D ^{3 –9} or 3-D ¹³ images, which were constructed by specially written planimetry computer programs. Second, the proposed tool not only detected the OD region and estimated its size, it also provided readings of the transverse and conjugate diameter—two commonly used parameters in retinal image analysis. Using this tool, we could further derive the normalized PPA size (i.e., the ratio between the PPA and the OD size) to explore its association with different eye diseases or conditions and establish a better understanding about the significance of the PPA development. Third, this tool has been automated, which means that not only could it reduce the dependence on the human assessor and thus potentially avoid problems associated with human errors such as fatigue, it could also be more cost effective for larger scale population-based screening. 

There remain some limitations within our method, however: First, our software stopped at undesired points on encountering irregular dark pixels before the OD boundary in good-quality images and less consistently in poor-quality images. This problem results in underestimation of the actual size. Second, the proposed algorithm estimates the sizes of OD and PPA regions, providing a means to measure the extent of PPA. It would be ideal if the software could also define the absolute shapes, allowing the patterns in PPA progression to be studied in different eye conditions. A possible way to address these above limitations would be to take into consideration additional local information (e.g., texture) and explore further the image fusion from multiple channels. 

Conclusions

We have demonstrated that PPA on a retinal image can be quantified by means of computer-aided software. Our proposed algorithm or estimation model, after calibration, achieved an accuracy of 94.0% (SD 5.88%) and 93.8% (5.26%) in defining the size of the OD and PPA, respectively, compared with the gold standard experienced human assessor. Our model also showed high reliability in estimating the size, with the correlation coefficient reaching 0.98 for both cases (OD and PPA). In terms of sensitivity, our model achieved 75.5% and 84.0% (tolerance, ±10%) in size estimation of the OD and PPA region, and higher when we increased the tolerance level. We have also investigated the robustness of our calibration method and found that our model consistently achieved a mean accuracy of more than 92.7%. In addition, the proposed software could also provide ophthalmologists additional information, namely transverse and conjugate diameter of the OD as well as the ratio between the OD and PPA size, with potential application in eye-screening programs. The methods developed so far are therefore promising as the basis for a fully automated prescreening technique that will prioritize images for subsequent expert human assessment. 

Supplementary Materials

Text s01, PDF - Text s01, PDF 

Footnotes

 Supported by Lothian NHS Endowments, Ref. 07425.

Footnotes

 Disclosure: C.-K. Lu, None; T.B. Tang, None; A. Laude, None; I.J. Deary, None; B. Dhillon, None; A.F. Murray, None