An efficient optic cup segmentation method decreasing the influences of blood vessels

Background Optic cup is an important structure in ophthalmologic diagnosis such as glaucoma. Automatic optic cup segmentation is also a key issue in computer aided diagnosis based on digital fundus image. However, current methods didn’t effectively solve the problem of edge blurring caused by blood vessels around the optic cup. Methods In this study, an improved Bertalmio–Sapiro–Caselles–Ballester (BSCB) model was proposed to eliminate the noising induced by blood vessel. First, morphological operations were performed to get the enhanced green channel image. Then blood vessels were extracted and filled by improved BSCB model. Finally, Local Chart-Vest model was used to segment the optic cup. A total of 94 samples which included 32 glaucoma fundus images and 62 normal fundus images were experimented. Results The evaluation parameters of F-score and the boundary distance achieved by the proposed method against the results from experts were 0.7955 ± 0.0724 and 11.42 ± 3.61, respectively. Average vertical optic cup-to-disc ratio values of the normal and glaucoma samples achieved by the proposed method were 0.4369 ± 0.1193 and 0.7156 ± 0.0698, which were also close to those by experts. In addition, 39 glaucoma images from the public dataset RIM-ONE were also used for methodology evaluation. Conclusions The results showed that our proposed method could overcome the influence of blood vessels in some degree and was competitive to other current optic cup segmentation algorithms. This novel methodology will be expected to use in clinic in the field of glaucoma early detection.

image has been widely used in many hospitals. Thus, it is possible to develop the CAD system used for ophthalmologic diagnosis based on fundus image processing technique [8][9][10][11].
The OD includes two distinct parts, namely, a circumjacent zone called the rim and a central bright region called the optic cup [12]. The optic cup can be divided into nasal and temporal regions. The former is generally occluded by the main blood vessels. Physiologically, the loss in optic nerve fibers (ONF) leads to the enlargement of optic cup which is called large cupping, and atrophy of neuroretinal rim which is called rim loss. Large cupping and rim loss are two important indicators of glaucoma. They could be detected by measuring the vertical optic cup-to disc ratio (CDR) [13]. Once more ONF disappear, the optic cup will become larger with respect to the optic disc, which corresponds to CDR increasing. In generally, an abnormal vertical CDR indicate a high risk of glaucoma. Therefore, automatic segmentation of the optic cup is crucial for CAD system for glaucoma.
To date, there are few studies on optic cup segmentation algorithms, which include thresholding-based methods [14][15][16], region growing [17,18], model-based methods (active contour models or snakes [19,20], level sets [21,22], and elliptical shape model [23]), anatomical evidence-based methods [24], and superpixel classification [25][26][27]. Optic cup segmentation is challenged because the intensity has a sudden change in areas which blood vessels pass across the cup-disc boundary. However, most studies haven't solved this problem effectively. Thus, inpainting of the blood vessels is an important step in optic cup segmentation.
This study aims to develop an efficient method to overcome the aforementioned problem in optic cup segmentation. Morphological operations were first performed to get the enhanced green channel image. Then, blood vessels were extracted and filled by using an improved Bertalmio-Sapiro-Caselles-Ballester (BSCB) model. Finally, local chart-vest (LCV) model was used to segment the optic cup.
The remainder of this paper is organized as follows: first, we introduced the image data used in this study and presented the proposed methodology, an efficient optic cup segmentation method decreasing the influences of blood vessels. The experimental results are then provided, followed by discussions and conclusions.

Data acquisition
Digital fundus images were acquired from local ophthalmologic hospital. The 24-bit color images were captured using digital fundus camera (Canon CR-DGi) with the array size of 1440 × 960 pixels. The photographic angle of the fundus camera was set to 60°, and the optic disc was adjusted at image center. The average optic cup boundary identified by two experts was taken as the ground truth for the following evaluation. All the experimented 94 images were acquired by the same instrument and the subjects were collected randomly to take the ocular disease screening. A total of 94 images including 32 patients with glaucoma (18 male, 14 female) and 62 healthy subjects (34 male, 28 female) were included in our experiment. For the 32 patients with glaucoma, the ages were ranged from 35 to 58 years old (45.23 ± 3.31 years). None of the left 62 normal subjects had history of hypertension, nor cardiovascular disease and diabetes. Their ages were ranged from 31 to 63 years old (43.51 ± 5.27 years). In addition, 39 glaucoma images from the public dataset RIM-ONE (An Open Retinal Image Database for Optic Nerve Evaluation) were also used for evaluation [28].

Optic cup segmentation
The edge of optic cup is more difficult to identify compared with that of optic disc, primarily because the image is blurred where the blood vessels pass across the optic cup. After preprocessing such as image enhancement, blood vessels were extracted and inpainting. Then the optic cup was segmented by using LCV model.

Preprocessing
Among the image components of color which are red (R), green (G) and blue (B), the G channel shows the optimal image contrast for the optic cup (see Fig. 1). Therefore, the G channel image U G was selected for the subsequent processing.
In the proposed algorithm, both top-hat and bottom-hat transformations were applied to enhance the image contrast [27]. Top-hat transformation refers to the subtraction of the opening operation result from the image itself. By contrast, bottomhat transformation refers to the subtraction of the image from the result of closing operation. Both top-hat and bottom-hat transformations are based on a predefined neighborhood or structuring element (SE). The above two transformations are illustrated as Eqs. (1) and (2), respectively.
Here, a structuring element with a size of 5 × 5 pixels was used (more details of size selection were provided in "Results" section). The above two transformations could be represented as Eq. (3), the boundary of optic cup become clearer after the operations (see Fig. 2). (1)

Blood vessel extraction
Since the intensities of blood vessels were lower than those of background and other structures, the intensities of blood vessels will become higher after median filtering. Therefore, the blood vessels could be possibly identified based on the intensity difference. In order to capture the change of image intensity, the contrast enhanced image U was subtracted from the median filtered image U med . Then the intensity differential image could be acquired which was represented by Eq. (4).
The binary image U D was generated according to the value of U sub [29]. Only the intensities of pixels with corresponding U sub > 0 were set to 1 according to Eq. (5). These pixels were considered belonging to blood vessels (see Fig. 3d).
In this study, the size of median filtering was set to 9 × 9 pixels (more details of size selection were provided in "Results" section).

Blood vessels inpainting
Bertalmio et al. [30] proposed the novel BSCB model for image inpainting in 2000 based on theory of partial differential equation (PDE). The BSCB model uses Laplace operator to measure the neighborhood information for image inpainting. It smoothly propagates the information to the same region along the direction of isophote. At the same time, an anisotropic diffusion function is adopted to prevent the prolongation lines from crossing one to another. Generally, the model comprises two steps: inpainting and diffusion. In this study, the improved BSCB model was developed, that is, the neighborhood intensities were used as the propagation information, instead of using only one single pixel.
After the above processing, the intensity will become more uniform within the area of optic cup. The influence caused by the blood vessels will be decreased. The following were the detailed descriptions about the processing procedure.
Assuming Ω is the region to be inpainted and ∂D is the boundary of Ω. The BSCB model can be described by Eqs. of (6) and (7) [31,32].
Equation (6) represents inpainting, where ∇L is the propagation information and T is the isophote direction. Equation (7) is used for diffusion, where k is the Euclidean curvature of the isophote and Ω ɛ is the dilation of Ω with a balling radius of ε, and g ɛ is a smoothing function of Ω ɛ . To be easier understanding, Eqs. of (6) and (7) can be comprehensively described as follows. Here is the output of the algorithm, and Δt is the improvement rate.
In the traditional BSCB model, propagation information L n (i, j) is substituted by the discrete Laplace operator, which is shown in Eq. (10).
Although texture could be the transfer information for image inpainting, it is not necessary in optic cup segmentation. In contrast, local image information is more effective. Therefore, using neighborhood mean value as transfer information may eliminate the influence of noise. According to the above deduction, u n xx i, j and u n yy i, j were replaced by u n xx i, j and u n yy i, j in the following vascular inpainting, represented by Eqs. (11) and (12) respectively. In this study, the size of neighborhood was set to 3 × 3 pixels.

Optic cup boundary identification
Wang et al. [33] proposed a LCV model which included local statistical information in level set based segmentation framework. In the algorithm, extended structure tensor (EST) was combined which intensity inhomogeneity could be decreased effectively. In this study, the above-mentioned LCV model was used for the following optic cup segmentation.
First, the centroid (x c , y c ) was selected from the region with relative high intensity, which was acquired according to the following criteria represented as Eq. (13) Here N is the total number of pixels within the target region. In this step, two parameters need dynamic adjusted in the LCV model, i.e., α and μ. Commonly, α was set to 0.1 or 1, depending on whether the image has intensity inhomogeneity. For μ, two corresponding values are adopted: 0.01 × 255 2 and 0.1 × 255 2 . If several targets need to be detected, μ should be small, and vice versa. Because the optic cup area has intensity inhomogeneity and only the optic cup be the target, the values of α and μ were then set to 0.1 and 0.1 × 255 2 respectively.

Evaluation
The experimental results were evaluated by three statistical criteria, namely, F-score (area-based), distance (curve-based) and vertical CDR, which were explained in "F-score (F)", "Distance (D)", and "Vertical CDR" sections.
Besides the proposed algorithm, the experiments were also performed by two other known methods, proposed by Joshi et al. [17] and Liu et al. [21]. There are two reasons why the above methods were selected to be compared. First, since the optic cup segmentation could be classified into framework of region based, edge based and hybrid, the proposed algorithm and another two selected methods all belong to region based method. Second, the proposed method aimed at decreasing the influence of blood vessels, which was similar to the other selected methods putting forward to the solutions.
In details, Joshi et al. [17] imposed the expected symmetry of optic cup region by setting threshold to recover the under-segmentation areas located in the blood vessels. Liu et al. [21] used ellipse to redraw the optic cup boundary after employing a combinative algorithm with level set and threshold setting, which was to weaken the noising by blood vessels.
For each method, the evaluation of F-score and distance were compared. The vertical CDR values were presented against manual segmentation results achieved by experts.

F-score (F)
The pixel-wise precision and recall values were computed to assess the overlap area between the computed region and the ground truth. These values were defined as Eqs. (14) and (15), where TP, FP, and FN represented the number of true positive, false positive, and false negative pixels respectively. The harmonic mean of the precision and recall values, called F-score (F), was computed to better appreciate the results. The F-score was expressed in Eq. (16), (14) precision = TP TP + FP (15) recall = TP TP + FN Given that both recall and precision are evenly weighted, the F-score value lies between 0 and 1, and the recall, precision, F-score should be all ideally close to 1.

Distance (D)
To assess the accuracy of the boundary, a curve-based evaluation is performed. Let C e be the boundary identified by the ophthalmologist and C m be the boundary achieved by the proposed algorithm. The distance (D) which was computed in pixels between two curves was expressed as Eq. (17), Here, n is the number of angles, the distance from centroid of C e to the points on C e in direction of θ was defined as d θ e , similar with the definition of d θ m . D should be close to 0 for an accurate algorithm.

Vertical CDR
To estimate the vertical CDR, the optic disc must be segmented ahead. Compared to optic cup, optic disc segmentation is relatively easier. In this study, to ensure the accuracy of the CDR and evaluate the effectiveness of the proposed algorithm, the optic disc was manually delineated by experts. The vertical CDR was calculated by Eq. (18) proposed by Gloster et al. [34], where C V represented the vertical diameter of optic cup, and D V represented the vertical diameter of optic disc. C V and D V were determined by the distances from the top to the bottom of these diameters, as shown in Fig. 4. In clinic, the risk of suffering glaucoma increases with the value of CDR.

Results
A total of 94 digital fundus images, including 32 glaucoma images, were experimented using the proposed algorithm. The algorithm was individually evaluated against manual segmentation results by expert-1 and expert-2, and also against the average expert marking called expert-X. For comparison, the results obtained by two other known methods of Joshi et al. [17] and Liu et al. [21], were also presented in evaluation of F-score and distance.
The optic cup segmentation results were shown in Fig. 5. In order to compare with the other above mentioned methods, more examples were shown in Fig. 6. The evaluation results of F-score and distance were shown in Tables 1 and 2. The values of CDR achieved by experts and the proposed method were also compared shown in Fig. 7. The statistical results of CDR for normal and glaucoma data were listed in Table 3.
The selections of parameters used in the algorithm were also evaluated. For the SE used in morphometric operations of preprocessing, the size from 2 × 2 pixels to 15 × 15 pixels have been experimented. Representative results were shown in Fig. 8 and   Table 4. Representative results of median filtering (window size is from 3 × 3 pixels to 15 × 15 pixels, interval of 2 × 2 pixels) were shown in Fig. 9 and Table 5.
From the results of Figs. 5 and 6, the proposed algorithm achieved the satisfied segmentation results which were better than the methods proposed by Joshi et al. [17] and Liu et al. [21]. Although both of them also considered the influence caused by blood vessels and proposed the corresponding ways trying to decrease the influence, our proposed method showed competitive to overcome the difficulty in areas with a lot of blood       9 Comparison results with different window size of median filtering. From left to right are optic cup segmentation results acquired with size of 5 × 5 pixels, 7 × 7 pixels, 9 × 9 pixels, 11 × 11 pixels and by the experts vessels. In addition, both the maximum F-score and minimum distance of our proposed method were acquired. That is, the accuracy of the proposed method is better than the other two. Furthermore, the standard deviations of F-score and distance were also lower which showed the robustness of the proposed method. From results of Fig. 8 and Table 4, the size of SE was suggested to be selected as 5 × 5. From results of Fig. 9 and Table 5, the window size of median filtering was suggested to be selected as 9 × 9.
In order to evaluate whether our proposed method is feasible and comparable to the most recent published optic cup segmentation algorithms [35][36][37], also 39 glaucoma images from the public database RIM-ONE were experimented. F-score and CDR acquired by our proposed method and these referenced algorithms were listed in Table 6. From the results, F-score acquired by our proposed method is acceptable and the corresponding CDR is more close to the manual results by experts.

Discussion
The proposed algorithm and the other two above mentioned methods used for comparison all put forward the solutions aimed at reducing optic cup segmentation errors caused by blood vessels. Given the figures listed in results, it can be concluded that the proposed method will give rise to higher F and lower D, which was hence better in optic cup segmentation. Joshi et al. [17] imposed the expected symmetry of cup region in nasal and temporal side after threshold processing to recover the vessel errors. A vertical axis of symmetry passing through optic disc center was considered and nasal region was obtained by mirroring the temporal region. Therefore, the segmentation results were decided completely by the segmentation quality of temporal region. However, most of the cups are not symmetrical about the vertical axis passing through optic disc center. What is more, as the acquisition conditions like photographic angle changes, the shapes of both sides are various which also causes the uncertainty of the relative size on both sides, leading to under-segmentation or oversegmentation of nasal side after the operation of symmetry. Liu et al. [21] applied ellipse fitting to redraw the cup boundary for weakening vessel influence in the segmentation process. Nevertheless, the ellipse fitting was implemented after employing combinative algorithm of level set and color intensity thresholds which didn't take  the vessel into account because of its low intensity. Thus only the minor remission of the errors caused by vessels could be reached after fitting, serious under-segmentation still cannot be avoided. By contrast, the proposed scheme firstly fills the blood vessels by an improved BSCB model, making the whole cup area more uniform for the following segmentation,then the optic cup is determined through the LCV model approach which can capture the unsharp boundary of cup, reducing under-segmentation or over-segmentation to a great extent compared with other methods. For the CDR evaluation in Table 3, the CDR values by the proposed algorithm are little smaller than those by expert, which indicates that the cup region determined by the proposed method was little smaller toward upper cup edge and lower cup edge than that determined by expert. The CDR values are little too small caused the misdiagnosis of four patients with glaucoma as normal in the experiment. The reason is that, for the dense vascular area near the cup edge, the grayscale is close to that of the neuroretinal rim after inpainting, leading to little local insufficient segmentation. But the difference of CDR between the proposed method and expert is small, thus the CDR was in the range of normal.
To sum up, we have proposed a novel automatic optic cup segmentation algorithm based on inpainting of blood vessels in this study. Compared with other techniques, the proposed methodology has several advantages.
1. The proposed algorithm repaired the vascular regions by using an improved BSCB model, reducing the errors caused by blood vessels to a large degree and improving the segmentation accuracy; 2. In LCV model, the centroid point located in area with higher intensity was selected as the starting point, which realized the automatic selection of the seed point. In addition, the loose selection of circular radius make the initial contour easily be automatically selected. 3. The proposed algorithm realized fully automatic segmentation with high accuracy, which reduced the human intervention of the traditional semi-automatic method.
The proposed method also has some limitations. First, the region of optic cup segmented by the proposed method was slightly smaller than that identified by experts, especially in the dense vascular area. It may be because parts of the image intensities near the edge of optic cup are close to these of the neuroretinal rim after inpainting. Accordingly, local insufficient segmentation will be happened. Second, the robustness of the algorithm is expected to improve. Since the used experimental data were collected from the same type of camera, more data from different type of camera should be tested and the robustness of the algorithm will be improved.

Conclusion
This study proposed a novel optic cup segmentation method based on inpainting of blood vessels using an improved BSCB model. The proposed algorithm realized fully automatic segmentation, which captured the optic cup boundary more accurately compared with other previous methods also dedicating to reduce the influence of blood vessels. Future work is needed to experimented more samples acquired from different type of camera and also developed more robust algorithms.