Retina Image Vessel Segmentation Using a Hybrid CGLI Level Set Method

As a nonintrusive method, the retina imaging provides us with a better way for the diagnosis of ophthalmologic diseases. Extracting the vessel profile automatically from the retina image is an important step in analyzing retina images. A novel hybrid active contour model is proposed to segment the fundus image automatically in this paper. It combines the signed pressure force function introduced by the Selective Binary and Gaussian Filtering Regularized Level Set (SBGFRLS) model with the local intensity property introduced by the Local Binary fitting (LBF) model to overcome the difficulty of the low contrast in segmentation process. It is more robust to the initial condition than the traditional methods and is easily implemented compared to the supervised vessel extraction methods. Proposed segmentation method was evaluated on two public datasets, DRIVE (Digital Retinal Images for Vessel Extraction) and STARE (Structured Analysis of the Retina) (the average accuracy of 0.9390 with 0.7358 sensitivity and 0.9680 specificity on DRIVE datasets and average accuracy of 0.9409 with 0.7449 sensitivity and 0.9690 specificity on STARE datasets). The experimental results show that our method is effective and our method is also robust to some kinds of pathology images compared with the traditional level set methods.


Introduction
Retina image segmentation as shown in Figure 1(a), also named as fundus image, has played an important role in diagnosing the pathologies in modern ophthalmology. The morphology of the vessel is a key indicator to detect retina disease in early stages and access the severity of the disease, such as diabetic retinopathy, age-related macular degeneration, and glaucoma. Manual screening process is a much time-consuming and tedious work and the rapidly growing medical imaging technology makes us get the medical image more easily and is cheap. Thus the automatic quantification of retina vessel is extremely needed for the early diagnosis of disease in the modern CAD (computer aided diagnosis) system. The vascular tree at Figure 1(b) is a prominent element that is required for the automatic analysis of the retina image in the CAD system. However, the illumination, resolution, and field of view (FOV) in the process of acquiring the medical image and the overlapping tissue in the retina make the segmentation a challenge. In particular, due to the reflectivity of the vessels and the complexity of the tissues in retina images, the intensity inhomogeneity often occurs in the retina image. However, the segmentation of such retina image vessels is still a tough problem.
The first group contains methods based on unsupervised method. Sopharak et al. [4] employed optimally adjusted mathematical morphology followed by morphological reconstruction to detect exudates in nondilated retinal images. A coarse level segmentation of exudates was carried out by Jaafar et al. [7] using local variation calculation of image pixels, in order to outline the candidate boundaries. The adaptive thresholding results that have been obtained using the coarse level segmentation are further refined using a morphological operation. Welfer et al. [13] used the contrast enhanced L channel of LUV color spaces to apply morphological operations as well as H-maxima transform for exudate detection. Harangi and Hajdu [14] have proposed an active contour-based region-wise method for identifying exudates. The second group contains methods based on supervised method. Most of these methods were prone to false positives near the vascular arch. Fraz et al. [11] introduced a novel method using bagged and boosted decision trees. The decision trees are used as classification model and the results of these weak learners are combined using bootstrap aggregation. Marín et al. [10] presented a supervised method for retinal vessel detection using a Neural network (NN) scheme. A multilayer feed forward neural network is adopted which consists of an input layer, three hidden layers, and an output layer. Franklin and Rajan [15] proposed the method of retinal vessel segmentation using multilayer perceptron Artificial Neural Network (ANN). The applied neural network has three layers of one input node, five hidden nodes, and one output node. A supervised hierarchical retinal blood vessel segmentation is presented by Wang et al. [16]. Convolutional neural network (CNN) performs as a trainable hierarchical feature extractor and ensemble random forest (RF) work as a trained classifier in that work. Supervised approaches learn a model to decide whether a pixel belongs to a vessel or not with the help of manual label, where the manual segmented training images are termed as the gold standard. Matched filter techniques [3] are the unsupervised method that uses a 2D kernel with Gaussian cross section to model the feature of the vessel in the retina image at some position and orientation. The matched filter response represents the existence of the vessel. Other unsupervised methods [17][18][19][20], which use the global and local average intensity information, employ deformable models for vessel segmentation.
In this paper, we propose a novel active contour method for the automatic segmentation of the blood vessel in the retina image. It is based on the Selective Binary and Gaussian Filtering Regularized Level Set (SBGFRLS) method [17] and the Local Binary fitting (LBF) term proposed by Li et al. [21]. In our method, the active contour is driven by the Selective Binary and local fitting term. Due to the local intensity term, the proposed model can avoid the numerical result trapped into the local minimum. The reinitialization method in our method is described in SBGFRLS model. This kind of reinitialization method simplified the computation during the evolution. Thus the proposed method is called combined global and local information method (CGLI).
Our method is tested using the following publicly available date sets: DRIVE [4] and STARE [3]. The STARE data sets are used to show the visual result of our method on pathology retina image and the DRIVE data sets are used to evaluate the performance of our methods compared with the golden standard provided. This paper is organized as follows. We firstly review some classic models and propose our method in Sections 2 and 3. In Section 4, the segmentation result of our model on retina image is provided. The performance of our method and a comparison with the previous model and other unsupervised vessel segment methods about the segmentation result are discussed in Section 5. Finally, the summaries are given in Section 6.

Previous Methods
In the following, some basic active contour models (ACM) are reviewed. Let Ω ⊂ 2 represent the image space, and : Ω → represents a gray image; ( ) : [0, 1] → 2 represents the contour . Image segmentation problem can be viewed as minimizing the energy functional and find a contour in Ω.

The Geodesic Active Contour (GAC).
Based on the Fermat principle in optics, the formula of the edge-based GAC model [5] is as follows: where is a smooth decreasing edge stopping function (ESF) where the ∇ * denotes the gradient of the image filtered by a Gaussian kernel , whose standard deviation is . Under the framework of the calculations of variation [6], the gradient flow [22] equation of (1) can be obtained and it can be used to achieve the numerical result of the minimization problem.
where is the curvature of the contour; → is the interior normal vector of the curve. is a constant used to control the evolution speed. Then the level set formulation is used as follows:

The C-V Model.
Chan and Vese proposed a region-based ACM model [23] based on the assumption that the image is intensity homogeneous. This method is also described as the well-known piecewise constant (PC) problem. Based on the variational level set framework [24], the CV model can be rewritten in the following form: in the equation is called Heaviside function and is the signed distance function (SDF). In practice, the is often approximated by a smooth version function by The difficulties in segmentation with retina images with CV and GAC models can be seen from Figure 2. The retina vessel in the first row in Figure 2 is the basic situation of the retina image within intensity inhomogeneity. For such images, it will not work when the thresholding method is used. In fact, no matter what threshold value is chosen, some part of the background or foreground is incorrectly identified as the foreground and the background which can be seen from the column three. At the same time, in the second row, it can be seen that the low resolution of the vessel cannot be detected in the CV model. This problem also exists in the GAC model. The example in Figure 2(c) shows the segmentation limitation of the CV and GAC models for the low resolution vessel detail in the retina image. In the meantime, the GAC model is much more sensitive to the image noise [5] which also limits its use in retina image vessel segmentation.

The Combination of the Global and Local Intensity Information Method.
The method is based on the GAC model. Unlike the traditional GAC model, a new edge-stop function which we called signed force is used to drive the contour to the ideal boundary which takes the global and the local intensity information into consideration. In the following, the signed force based GAC will be described. The description of the edge-stop function is as follows: This new ESF is used to replace at (1). The new level set formulation is obtained as follows:

The Construction of the Local Force.
In C-V model, the image is regarded as intensity homogenous. However the intensity inhomogeneity, as shown in Figures 2 and 3, is very common in retina images. So the local statistical information is required to get the segmentation result in such kind of images. In the local binary fitting (LBF) model [19], two local fitting terms 1 ( ) and 2 ( ) are used to present the intensity difference between inside and outside of the contour. The two local fitting terms make the level set method able to work even in some intensity inhomogeneity. The energy function is described as follows: Due to the prosperity of Gaussian kernel function , the contribution of the intensity ( ) to the energy at (9) decreases to zeros when moves away from the center . By minimizing (9), we can obtain 1 and 2 as follows: This local fitting term produced by Li overcomes the limitation of the region-based level set method, but it is much more sensitive to the position of the initial contour and it is also easily trapped into the local minimum during the curve evolution. So based on the analysis above, the local fitting term is regarded as the local information term in our CGLI method; therefore our method can work effectively in the intensity inhomogeneous images and can give response to the low resolution vessel detail shown in Figure 2. The local force part of our method is based on Li's method as shown in the following: The local fitter terms 1 ( ) and 2 ( ) represent the outside and inside intensity of the pixel point near the boundary. This term makes our method able to work well in very low contrast image and is proper to segment the capillary vessel in retina image, but it is not enough. A global term is required to prevent the model from trapping into the local minimum value. This global force was proposed in the SBGFRLS model [17]. Equation (12) will be explained in the following. It is assumed that the intensity outside and inside the vessel is homogeneous. It is obvious that min( ( )) ≤ 1 , 2 ≤ max( ( )). Hence, there is The intensity of the black vessel-like part of the image is the min(I(x)) and the gray part is the max( ( )). Based on (13), the signed force at the black part will be minus if the contour is not at the desirable boundaries of this vessel-like part. The average intensity of the vessel-like part will be larger than min( ( )) if the gray part could increase the average intensity inside the contour. As a result, the vessel-like part of the level set function shown in Figure 4(b) will be pulled down because of the signed force of image and the other area of the image will be pulled up. The zero level set of the will be the final segmentation contour.
All the energy functionals in the level set methods have an energy term to keep the result of the segmentation result contour smooth. Li et al. [19] add a regularization term in the energy to regularize the level set function; Kass et al. [25] used a curvature term to keep the contour smooth. In the GCV model a Gaussian filter is used to regularize the level set function after each iteration. The second term in the energy function at (8) is also the regularization term. The regularized way for Gaussian filter is chosen in our method, so the second term and the curvature term can be all omitted. The final formulation of our method is as follows: As shown above, the global force can work in the intensity homogenous images. The local force based on the local fitting term only indicates the changes at the local region described by Li et al. [19]. Both of the SBGFRLS model and the LBF model have the ability to segment different kinds of images. The CGLI method takes the global force and the local force into consideration. Thus our method has the ability to deal with intensity inhomogeneous retina images and also can get subpixel accuracy and is robust to the position of the initial contour if the weight of the force is set appropriately.

Numerical Implementation.
Here it is assumed that the target image is made of two parts, background area Ω 1 and foreground area Ω 2 . In the traditional GAC model, it is required for initial level set function to be an SDF in order to ensure the stability during the evolution, and reinitialization is also required in the evolution. To avoid this problem, Zhang et al. [17] proposed a novel method using a Gaussian filter to smooth the level set function to regularize the level set function based on the mathematic theory of scale-space [6]. Firstly, we choose the level set function and the regularized Heaviside function ( ) and its derivative ( ) in our method: The pseudocode of our CGLI level set method is as follows: :

Initialize a closed contour C in the two dimension image space Ω iterationTimes←150
For times < iterationTimes

Data Sets.
The test of proposed method is using one of the popular public sets of retina fundus images called DRIVE [7]. These public data sets are popular because they provide the ground truth images that are manually segmented by different experts that make it easy for us to measure the performance of different methods. The DRIVE database is made of 40 images, including a training set and a test set, each of which contains 20 images. The mask of the retina image of the FOV area was also provided. In the experiment, the images in the training set are used to make comparison with other methods. The performance of our proposed method is measured by comparing the automated segmentation result of the retina image with corresponding ground truth images segmented manually by experts in the training data set.

Performance Measurements.
The flowchart of the proposed method is shown in Figure 5. For each test image the final zero level set is used to obtain the vessel profile. The contour obtained using our method can divide the image pixels into two classes: vessels and nonvessels. Then, we compare the segmented binary image with their corresponding ground truth image by computing the following four performance measurements. The pixels that belong to a vessel

Parameters Description
Initialize the level set, = 1, in our method Scale parameter in Gaussian kernel = 4 (default) The weight of the local force Δ Time step Δ = 2 (default) = 3.1415 is a constant It controls the smoothness of the Heaviside function = 0.5 in our method Iteration times in the level set method The balloon force, which controls the speed during the curve evolution in GAC method in the ground truth image are classified as vessels and counted as true positives (TP); otherwise they are counted as false negative (FN). The pixels that belong to the background in the ground truth image are classified as nonvessels and counted as true negative (TN); otherwise they are counted as false positive (FP).
In order to compare with other vessel extraction algorithms, the accuracy (ACC), sensitivity (Se), and specificity (Sp) are calculated. These metrics are defined as follows: where N = TN + TP + FN + FP. As shown in Figure 1, the dark background outside the FOV area of the retina can be easily detected. In our experiment result, only the pixels of the FOV area are taken into consideration when we compute the performance value.
Several experiments are carried out using different sets of parameters for the application of the CGLI model. The proposed method has been tested with synthetic and real retina images. In addition, the proposed method is implemented in MATLAB 2014b on a 2.4 GHz PC. The parameters are described in Table 1.
The results for the X-ray images are shown in Figure 6, one X-ray image of vessels. The images are typical images with the intensity inhomogeneity. In these two images, it can be seen that parts of the vessel boundaries are quite weak. The segmentation result shows that our method has the ability BioMed Research International to deal with some kinds of images with low resolution and get vessel-like structures in the images containing different tissue. Unlike the experiment result described in Li et al. 's method [19], our method does not rely much on position of the initial contour and can get more low quite weak boundaries. Figure 6 shows the performance of the LBF model and SBGFRLS model separately in real retina image. It can be seen that both the SBGFRLS model and LBF model could not get the desirable vessel details by their own. In the SBGFRLS model, the curve is driven by the global intensity information, which is the same as the CV model. The segmentation result is shown at the first row. This model used a new way to reinitialize the level set function, which can reduce the time during the evolution. In the LBF model, Li used the local statistical information to overcome the intensity inhomogeneous problem. As shown in Figure 7, the curve does not converge at the describable location. The final segmentation result depends much on the position of the initial contour. Figure 8 shows the result of the CGLI method for the real retina images in DRIVE database. It also shows the details of the segmentation. The difficulties described in Figure 6 can be solved. Not only the desirable segmentation result can be gotten, but also our method can converge at acceptable iteration times. In LBF models [18], it is not robust enough to the position of the initial position and does not converge quickly.

Pathology Retina Image Segmentation Details.
Our method also has the ability to deal with some pathology images, as shown in Figure 9. Figure 9(a) is the original image from DRIVE and one part of its segmentation result from left to right. We only set one parameter = 0.2. Figure 9(c) shows the detail of the segmentation result. It shows that our method can get the subpixel vessel profile even if the contrast of the vessel is bad as shown in the right part of the row and it also can ignore the disease area shown in the right part of the row.

Discussion
The ACM models have widely been used in image segmentation in the decays. The active contour models have lots of advantages in image segmentation when compared with traditional methods, such as match filtering and being multiscale. Firstly, active contour models can get subpixel accuracy at the object boundary. Secondly, it can be easily formulated under the framework of the energy minimization and also allows us to add the prior knowledge. Thus it can be extended in a more natural way unlike the multiscale method or match filtering method which relies on the specified models.
The first row in Figure 10 shows the different position of the initial contour (white triangle) and the second row is its segmentation result. It can be seen that our method is robust to initial condition. It almost gets the same result. So our method is more robust to the position of the initial contour.  It is also necessary to talk about the parameters in our method listed in Table 1; all the parameters are predefined in our method except . decides the weight of the global force in the in (7). So we should set relatively small for more details and otherwise large for less details. And is the only parameter we care about in our method. Two different values of = 0.2 or 0.5 are used to show the influence of . In Figure 11, the first row = 0.2, and the second row = 0.5. Less vessel detail means larger . Generally speaking, if we want to deal with the low quality of retina images, should be set larger.

Compared with Other Unsupervised
Method. The traditional level set method is not performing well in the retina image segmentation problem. The SBGFRLS model is a less time-consuming than level set method. The LBF method is very effective in MRI image segmentation. Both of them failed in the retina image segmentation, as shown in Figure 11; the method we proposed can run fast in image processing and it is also effective in intensity inhomogeneous images. As a result, it can perform well in retina vessel extraction. We achieve on the DRIVE data sets and compare them to other supervised and unsupervised method. The performance of our method is acceptable as shown in Table 2. Our method can also achieve better segmentation results for other medical images with intensity inhomogeneity.

Conflicts of Interest
The authors declare that they have no conflicts of interest.