Contrast based circular approximation for accurate and robust optic disc segmentation in retinal images

Materials and performance evaluation

Two publicly available databases of retinal images have been considered to evaluate the performance of the proposed localization and segmentation methods. The main reason for this choice is the widespread use of these datasets and also the availability of a ground truth which can be used in the same conditions for any researcher. There exist other popular databases but the ground truth is not always accessible.

The Messidor database (Decencière et al., 2014) was created in the framework of diabetic retinopathy screening and diagnosis. It consists of 1,200 images acquired in three different ophthalmology departments using a 3CCD color video camera on a Topcon TRC NW6 non-mydriatic retinograph with a 45°  FOV and three different resolutions: 1,440 × 960, 2,240 × 1,488 and 2,304 × 1,536 pixels. This dataset has been used in most of the recent works mainly because of its size, much bigger than other public datasets, and due to the availability of a ground truth with contours traced by an expert from the University of Huelva (2012).

The ONHSD database (Lowell et al., 2004) consists of 99 fundus images with a resolution of 760 × 570 pixels taken from 50 patients in the context of a diabetic retinopathy screening program. The images were captured using a Canon CR6 45MNf camera with a 45°  FOV. As suggested by Lowell et al. (2004), only a subset of 90 images is considered due to the very bad quality of the discarded images. A ground truth is also provided.

In order to evaluate the performance of the proposed methods, different measures are proposed. The OD localization is usually considered as successful if the estimated center lies within the boundary of the corresponding ground truth and that is the approach followed in this work. The error in the localization is evaluated in terms of the error distance, $D\left(c_{\exp },c_{\mathrm{real}}\right)$, between the real location and the experimental one. Since the size of the retinal images under consideration may be different, a normalized error distance, $D^{\ast }\left(c_{\exp },c_{\mathrm{real}}\right)$, becomes a more convenient measure. (1)${\displaystyle D^{\ast }\left(c_{\exp },c_{\mathrm{real}}\right)=\frac{D\left(c_{\exp },c_{\mathrm{real}}\right)}{R}}$ where R is the radius of a circle with an area equal to the size of the corresponding ground truth mask.

For the purpose of evaluating the performance of the segmentation method, the Jaccard and Dice coefficients as well as the mean average distance between the true and estimated OD boundaries have been chosen.

The Jaccard and Dice coefficients are measures of the similarity between two sets. The Jaccard coefficient, JC, is defined as the ratio between the intersection and union of the automatic segmentation of the OD, S_OD, and the corresponding ground truth, S_Truth. (2)${\displaystyle JC\left(S_{OD},S_{\mathrm{Truth}}\right)=\frac{\left|S_{OD}\cap S_{\mathrm{Truth}}\right|}{\left|S_{OD}\cup S_{\mathrm{Truth}}\right|}.}$ The Dice coefficient, DC, is defined as twice the ratio between the intersection of S_OD and S_Truth, and the sum of their sizes. (3)${\displaystyle DC\left(S_{OD},S_{\mathrm{Truth}}\right)=\frac{2\left|S_{OD}\cap S_{\mathrm{Truth}}\right|}{\left|S_{OD}\right|+\left|S_{\mathrm{Truth}}\right|}.}$ As opposed to the JC and DC coefficients which are intended for evaluating the degree of overlap between regions, the mean average distance, MAD, is used to measure the distances between the closest points on their boundaries. Let us denote by $B_{OD}=\left(a_1,a_2,\dots ,a_N\right)$ and $B_{\mathrm{Truth}}=\left(b_1,b_2,\dots ,b_M\right)$, the finite sets of points belonging to the boundaries of the segmented OD and the corresponding ground truth, respectively. The MAD between the two sets is defined as follows: (4)${\displaystyle MAD\left(B_{OD},B_{\mathrm{Truth}}\right)=\frac{1}{2}\left(\frac{\sum _{i=1}^{N}d\left(a_i,B_{\mathrm{Truth}}\right)}{N}+\frac{\sum _{j=1}^{M}d\left(b_j,B_{OD}\right)}{M}\right)}$ where $d\left(a_i,B_{\mathrm{Truth}}\right)=\stackrel{\min }{j}\left|b_j-a_i\right|$ and $d\left(b_j,B_{OD}\right)=\stackrel{\min }{i}\left|a_i-b_j\right|$.

It is important to mention that the parameters of the methods for OD localization and segmentation have been experimentally determined using 40 images from the ONHSD dataset. The performance evaluation has been carried out using the whole ONHSD dataset and the Messidor database.

Method for optic disc localization

The method for the estimation of the optic disc center consists of two main steps: creating a mask based on vascular information to shrink the search space, and filtering the image with a detector which combines vascular and brightness information. The only preprocessing which was found to have a significant impact on the results consisted of resizing the original image using cubic interpolation so that the diameter of the retina becomes equal to 540 pixels as in Marin et al. (2015). By doing so, the proposed methodology can be used without the necessity of adapting the parameters of the algorithms. In the case of using a different angle for capturing the images, the parameters should be readjusted accordingly.

Shrinking the search space using vessel information

The first step exploits the fact that the OD is the entry point for the major blood vessels that supply the retina. Consequently, the OD region presents usually a high vessel density and can be seen also as a convergence point of this structure. The complementary of the green channel in the RGB color space, G_c, is taken as the input image for this stage since it is known to provide a good vessel-background contrast. Vessel enhancement, V_E, is carried out by applying the top-hat operator T₈ to G_c using a disc of 8 pixels in radius as the structuring element. (5)${\displaystyle V_E=T_8\left(G_c\right).}$ Vessel density, V_D, is calculated as the difference between the output of two averaging filters on V_E, an averaging filter of size 80 × 40 and an averaging filter of size 80 × 120. (6)${\displaystyle V_D=A{v}_{80\times 40}\left(V_E\right)-A{v}_{80\ast 120}\left(V_E\right).}$ The result is normalized dividing by its maximum. The rationale behind this is not only to capture the high density of vessels in the OD region but also take advantage of the fact that the density on both sides of it is significantly lower. A binarized version, $V_D^b$, is obtained by thresholding V_D using a threshold of 0.3.

The convergence of the branches of the vascular tree, V_C, is estimated by finding the intersections of lines which are used to approximate the branches of this tree. The lines are the result of the application of the Hough Transform (HT) to the output of a Canny edge detector, ED_Canny, computed on the vessel-enhanced image, V_E, so that the value of each pixel of V_C corresponds to the number of Hough line intersections at that point. Due to the noisy nature of this calculation, the resulting image is filtered to obtain an averaged version using a circular averaging filter of radius 40. (7)${\displaystyle V_C=A{v}_{r=40}\left(HT\left(E{D}_{\mathrm{Canny}}\left(V_E\right)\right)\right).}$ A binarized version, $V_C^b$, is obtained by thresholding V_C using a threshold of 0.3.

The two previous masks are combined using a logical AND operator to provide the final constraint mask, $V_{DC}^b$, based on vascular information. (8)${\displaystyle V_{DC}^b={\mathrm{logical}}_{\mathrm{AND}}\left(V_D^b,V_C^b\right).}$ Figure 1 shows an example of the different operations involved in the calculation of $V_{DC}^b$. It is interesting to remark the role of $V_C^b$ in the discarding of exudates. Very few lines are detected by the Hough Transform in these unwanted regions as shown in Figs. 1E and 1F.

Building the optic disc detector

The intensity image, I, in the second step is calculated as the sum of the three RGB color channels so that the OD usually appears as the brightest region but not always due to artifacts. In fact, the vessel pixels in I which are within the OD may take low values. In order to compensate for this low intensity, a vessel mask, V_M, is superimposed on the intensity image and the values of the pixels in the mask are set to the maximum value in I to provide the input image, I_c, for the OD detector. The vessel mask is obtained by thresholding V_E to retain only those pixels with high values. (9)${\displaystyle V_M=\left\{\begin{array}{cc}\hfill max\left(I\right),\hfill & \hfill \text{if}{V}_E>{\mathrm{percentile}}_{98.5}\left(V_E\right)\hfill \\ \hfill 0,\hfill & \hfill \text{otherwise}\hfill \end{array}\right\}}$ (10)${\displaystyle I_c\left(p\right)=max\left(I\left(p\right),V_M\left(p\right)\right)}$ where p is any pixel in the image.

The OD detector, D_OD, is finally calculated as the difference between the output of two averaging filters on I_c, a circular averaging filter of radius 40 and an averaging filter of size 80 × 160. (11)${\displaystyle D_{OD}=A{v}_{r=40}\left(I_c\right)-A{v}_{80\ast 160}\left(I_c\right).}$ The x and y coordinates where D_OD reaches its maximum value after the application of the constraint mask $V_{DC}^b$ are taken as the best approximation to the optic disc center. Figure 2 shows an example of the different operations involved in this second step.

Method for optic disc segmentation

The method for OD segmentation has been implemented on a subimage of the original RGB retinography by cropping it around the estimated center. The size of the cropping window depends on the resolution of the image under consideration and has been chosen to completely contain the OD and include a certain amount of the background around. In this way, a more robust and efficient segmentation can be achieved due to reduced space for search and consequently reduced number of artifacts present in the whole image.

The intensity image I defined in the previous section is taken as the basis for the calculation of the contrast between the regions of interest. For that purpose, a set of concentric circles were considered with center $\left(x_c,y_c\right)$, coincident with the center of the window, and increasing radii in increments of 1%, ranging from 25% to 45% of the side length of the cropping window. For each reference radius r, an outer and inner ring are defined as shown in Fig. 3.

The outer ring is the region lying between the two circles with the biggest radii and the inner ring is the region lying between the two circles with the smallest radii. The thickness of the ring was set to 4% of the side length of the cropping window. A contrast measure is defined as the subtraction of the average intensity in the outer ring from the average intensity in the inner ring: (12)${\displaystyle CM\left(r\right)=av{I}_{\mathrm{inner}}\left(r\right)-av{I}_{\mathrm{outer}}\left(r\right).}$ The circle for which CM attains the maximum value is taken as the best approximation to the OD contour under these assumptions.

We have found that the pure circular approximation obtained in this way may become too rigid in some cases. For this reason, an alternative solution is presented which retains the original simplicity but allowing some more flexibility in the form of considering four fractions of a circle instead as shown in Fig. 4. This is the approach which has been followed for OD segmentation.

The same procedure as described above is applied to each of the circle sections so that four different radii $r_u^{\max },r_d^{\max },r_l^{\max },r_r^{\max }$ are computed for which the maximum contrast is reached in each section. The final OD boundary is computed as a circle whose center $\left(x_c^{\max },y_c^{\max }\right)$ and radius r^max are calculated as follows: (13)${\displaystyle x_c^{\max }=2x_c+r_r^{\max }-r_l^{\max },y_c^{\max }=2y_c+r_d^{\max }-r_u^{\max }}$ (14)${\displaystyle r^{\max }=\frac{r_u^{\max }+r_d^{\max }+r_l^{\max }+r_r^{\max }}{4}.}$

Optic disc localization

The method proposed for optic disc localization achieves a 100% success rate in the ONHSD database and a 99.58% success rate (5 failures) in the Messidor database. Figure 5 shows the values of $D^{\ast }\left(c_{\exp },c_{\mathrm{real}}\right)$ obtained for the Messidor and ONHSD datasets with mean values of 0.089 and 0.099, respectively. Table 1 shows a comparison with other methods for both datasets.

Figures 6 and 7 show some examples of OD localization in the Messidor and ONHSD datasets.

Optic disc segmentation

Table 2 shows the performance of our method for optic disc segmentation in terms of the measures described in the materials and performance evaluation section. As it can be seen, two different scenarios have been considered, a completely automatic procedure where the OD center is estimated by means of the localization method explained in this paper and another ideal situation in which the OD center is obtained from the provided ground truth, i.e., manual localization. We have done so to remark the influence of OD localization in the proposed segmentation method and to set a bound on the maximum reachable performance. The average computation time is 1.28 s per image for the automatic version and 0.07 s. for the manual version so it becomes clear that most of the time is consumed in the localization of the OD. A comparison with other methods has also been included in the table. It should be noted that not all the measures are available for the different methods considered.

The percentages of images for different intervals of the overlapping coefficient JC are shown in Table 3 as well as a comparison with other methods for the Messidor dataset due to the unavailability of data for ONHSD.

Figure 8 shows some examples of correct and incorrect segmentations in the Messidor dataset and Fig. 9 shows some examples of correct and incorrect segmentations in the ONHSD dataset using the proposed method.