A New Approach to Segment Both Main and Peripheral Retinal Vessels Based on Gray-Voting and Gaussian Mixture Model

Vessel segmentation in retinal fundus images is a preliminary step to clinical diagnosis for some systemic diseases and some eye diseases. The performances of existing methods for segmenting small vessels which are usually of more importance than the main vessels in a clinical diagnosis are not satisfactory in clinical use. In this paper, we present a method for both main and peripheral vessel segmentation. A local gray-level change enhancement algorithm called gray-voting is used to enhance the small vessels, while a two-dimensional Gabor wavelet is used to extract the main vessels. We fuse the gray-voting results with the 2D-Gabor filter results as pre-processing outcome. A Gaussian mixture model is then used to extract vessel clusters from the pre-processing outcome, while small vessels fragments are obtained using another gray-voting process, which complements the vessel cluster extraction already performed. At the last step, we eliminate the fragments that do not belong to the vessels based on the shape of the fragments. We evaluated the approach with two publicly available DRIVE (Staal et al., 2004) and STARE (Hoover et at., 2000) datasets with manually segmented results. For the STARE dataset, when using the second manually segmented results which include much more small vessels than the first manually segmented results as the “gold standard,” this approach achieved an average sensitivity, accuracy and specificity of 65.0%, 92.1% and 97.0%, respectively. The sensitivities of this approach were much higher than those of the other existing methods, with comparable specificities; these results thus demonstrated that this approach was sensitive to detection of small vessels.


Introduction
Retinal fundus images are used to diagnose certain eye diseases and some systemic diseases. Blood vessels are one of the most important components in the retina, and abnormal vessels can indicate the presence of various diseases such as diabetes, glaucoma, retinopathy, obesity, voting and Gaussian mixture model (GMM) method to segment the vessels in fundus retinal images. First, we obtain a vessel-enhanced image by combining a 2D-Gabor filter result with a gray-voting result. Second, we classify the pixels of the vessel-enhanced image into different groups using a GMM. The group that contains vessel information is regarded as the preliminary vessel segmented results. Then, the result of another gray-voting process on the enhanced image is used to complement the preliminary vessel segmented result. Finally, we use a fragments elimination algorithm to remove the pixels that do not belong to vessel fragments. The block diagram for the steps of the proposed method is shown in Fig 1. The remainder of the paper is organized as follows: the method proposed for blood vessel detection is presented in Section 2; the method to identify the vessel complement is presented in Section 3; the experimental results and comparison are given in Section 4; and the final conclusions are given in section 5.

Proposed Method Preprocessing
There are red, blue and green channels in an RGB fundus retinal image. The green channel shows the best background/vessel contrast [9], and its signal noise ratio is higher than the other channels. In this study, the green channel of a fundus retinal image is used as the input to the subsequent step in the image preprocessing stage. The green channel can be divided into the background and the foreground. Optic disk and fovea all belong to background; however, an optic disk has a higher mean gray value than that of the entire image, whereas fovea has a lower mean gray value. The foreground is primarily considered to be a vessel that shows a different gray value in a different region of a retinal image.
1) 2D-Gabor filter. To extract the primary structure of a vessel, the original green channel of the image is passed through a 2D Gabor filter [9]. The 2D Gabor filter uses a Gaussian kernel function modulated by a sinusoidal plane wave, which is very sensitive to a retinal vessel because changes in the gray level between the vessel and the background are shown as a Gaussian distribution. Designing of the 2D-Gabor filter is accomplished as follows.
First, a continuous wavelet transform T ψ (b,θ,a) is defined as the scalar product of f with the transformed wavelet ψ b,θ,a : Where the parameters b, θ and a describe the translations, rotations and dilations, respectively: ψ Ã is the complex conjugate of 10°: and C ψ and 10°denote a normalizing constant and an analyzing wavelet, respectively.
The 2-D Gabor wavelet is defined in formula (2): The fast Fourier transform is used to implement the 2-D Gabor wavelet transform: is a 2×2 diagonal matrix, and k 0 is a vector that defines the frequency of the complex exponential. Next, for each pixel in the retinal image, we extract the maximum modulus of all orientations based on formula (4): where θ in the Gabor wavelet transform spans from 0 to 170°at steps of 10°. In this study, we set the parameter a to be constant and equal to 3, which was determined by comparing experimental data. Finally, we obtain the result of 2D-Gabor filtering M ψ (b,a) which is denoted as I Gabor in this study, as shown in Fig 2(C), from which it is shown that the primary structure of the vessel was extracted, and little noise is present in the background.
2) Gray-vote algorithm. Although a 2D-Gabor filter can extract the main vessel structures in a retinal image, details, especially small vessels, are often lost. To obtain more small vessel information from the green channel retinal image I green , we propose a gray-voting method, which is described below. Because the gray values of the vessel pixels change dramatically in different regions, a local gray analysis is necessary. The gray-voting method can enhance small vessels that have a similar gray level distribution to the background. Parameter k is a gray transition scale, which is used to obtain the small vessel information from the gray-voting process. For each pixel, we used a window in size of m × m (i.e., m is the pixel number in a row or column), where the gray value of the center point in the window is denoted as Center(i, j). Comparison is the value that is used to compare with the other pixels' gray values in the window. In Eq (6) below, Neighbor is a set of the other pixels except the center pixel in the m × m window, where the initial value of Num 1 and Num 2 is zero; Num 1 is the pixel number where the gray value is larger than Comparison in the m × m window; and Num 2 is the pixel number where the gray value is smaller than Comparison in the m × m window. P vote (i, j) is the outcome of the grayvoting process, and L and N are the maximum and minimum normalization gray values in the m × m window, respectively. We obtain P vote (i, j) using Eqs (5)-(8): Comparison ¼ Centerði; jÞ À k: ð7Þ Because this gray-voting algorithm is sensitive to a slight change in the gray-level in the window, the outcome P vote (i, j) can detect small vessel structures. Compared to the original green image, this gray-voting algorithm, as shown in Fig 2(B), enhances small vessel details but also noise fragments.
3) Image fusion. In the first two sections of this chapter, we obtain the 2D-Gabor result I Gabor and the gray-voting result I vessel by using a 2D Gabor filter and the proposed gray-voting algorithm. The 2D-Gabor result I Gabor and the gray-voting result I vessel contain the main vessel structure and the small vessel information, respectively. To obtain an image with both the main vessel structure and the details, we fuse I Gabor and I vessel . The fusing result I gv is obtained by Eq (9): As shown in Fig 2(D), the fused image I gv , which shows a significantly smoother background, has better connectivity than the gray-voting result I vessel (Fig 2(B)), and contains more detail than I Gabor (Fig 2(C)).

GMM classifier
As shown in Fig 2(D), the fused result I gv may be composed of pixels of both vessels and noise. Because vessel pixels usually have higher gray-levels than noise, we use GMM [15] to classify the pixels in the fused result I gv . In this process, GMM is adopted to analyze the gray level distributions of the pixels in the fused result I gv . First, we apply the K-means clustering method to calculate the K centers μ i (m = 11) and the variances k = −3 (L ¼ M max ðm 2 À1Þ ) of the pixels, which are used to initialize the Gaussian mixture distributions. Then, each distribution is labeled with a weight as specified below. The expectation of GMM is represented by Eq (10): where m 2 is the number of pixels to be classified; and M max , M min and m × m describe the weight, mean and variance of the ith Gaussian distribution, respectively. Next, the EM algorithm is used to obtain the maximum likelihood estimate of the GMM parameters, including the weights, means, and variances. In the EM process, the GMM parameters are iteratively updated using Eq (11) from their initial values, which were derived using K-means clustering: where N is the number of pixels to be classified and t denotes the tth iteration. This iterative update is performed until the log likelihood log Pðx j jm; dÞ is convergent. Then, the cluster is built from each of the K Gaussian distributions under the parameters derived by the EM algorithm. For a pixel, if it generates the maximum likelihood in the ith Gaussian distribution, it is assigned into the ith cluster C i using Eq (14) and (15): ( Furthermore, the expected value μ i of the ith Gaussian distribution provides a center of all the pixels in layer C i .

Post-Processing
In the proposed method, the post-processing of the vessel segmentation can mainly be divided into two parts. The first part complements the GMM classifier result using a gray-vote image that contains rich small vessel details. The second part eliminates the fragments that do not belong to the vessel using morphological characteristics.

Vessel complementation
As shown in Fig 3(C), the vessel cluster I GMM contains the main vessel structure and some small vessel branches. However, some small vessel branches are broken into fragments. To address this issue, we used another gray voting processing on the fusion result I gv with different parameters to obtain a complementary image I c (Fig 4(B)), which contains the details of the rich small vessel; we then used the fragments to link the broken vessels of the vessel cluster I GMM . I com (Fig 4(C)) is the binary result of the complement image I c . The binary image I com shows rich details of the vessels that are used to complement the vessel cluster I GMM .
I com and I GMM are the complementary image and a marker image, respectively, in the image complementation process. In this section, T seed is the threshold that is used to estimate a broken vessel fragment using the number of pixels in the fragment. T fragment is the threshold for estimating the complementing fragments. First, we search the vessel fragments in the marker image. If the number of pixels in a fragment is less than T seed , the fragment is considered to be a broken vessel fragment; then, we set these pixels to be the seed points and count the number of pixels in the fragments that have the same seed labels in the complementary image. If the number of pixels in a fragment is less than the value of T fragment , then the fragment is regarded as a vessel fragment and complemented in the marker image.
The pseudocode of this vessel complement method is presented below. input: I com (complement image), (marker image) if length(xx)< = T fragment % (Judging whether the number of complement image fragment pixels less than T fragment ) for j = 1:length(xx) I C−GMM (xx(j),yy(j)) = 1; % (Complementing the fragment to the output) end end end end end

Fragments elimination
The result of vessel complementation I C−GMM , is shown in Fig 5(C), and it contains many fragments; however, not all of the fragments belong to the vessel. Because the image complementation process is based on the number of pixels in the fragments, some tiny noise fragments of the binary image I com are considered vessel fragments. Thus, we used the Cemal Kose and Cevat Ikibas's fragments elimination method [5] to solve this problem. First, we search the fragments in the complementary vessel I C−GMM and apply the seed fill algorithm to calculate the number of pixels in the fragments. Then, we calculate the maximum coordinate values along the x-and Y-axes and use the larger of these two values. The number of pixels in a fragment is used to calculate the squareness rate of the fragment, which is described by Eq (16). We determine whether the fragment belongs to a vessel or not using Eq (17).
where S_Rcs is the squareness rate, F_Scs is the number of pixels in a fragment and mx is the larger of the two maximum coordinate values on the x-and y-axes: ( Fig 5(A) shows the complementary image obtained in Section 2.2, which contains many small vessel fragments that could be used to complement the broken vessels. Fig 5(B) shows the outcome of the GMM classifier; the fragments of this image could be regarded as the vessel seed fragments, which are used to locate the vessel position and search the other vessel fragments in the binary image I com . Fig 5(C) shows the outcome of the complementary processing, and Fig 5(D) shows the result after fragments elimination.

Experimental parameters
In the gray-voting process, parameter k has a significant influence on the gray-voting result. Fig 6A-6D show the gray-voting results when k is set to 0, 3, 5 and 10, respectively. As shown in Fig 6, when k is set to 0, both noise and small vessel pixels are detected. As k increases, the noise decreases; however, the pixels of small vessels cannot be detected. We perform the proposed gray-voting algorithm on each pixel of the green channel image I green , as shown in Fig 2(A), and obtain the gray-voting result I vessel , as shown in Fig 2(B); this figure shows that the structure of the vessel, the branch details and some noise are all contained in the gray-voting result I vessel . From multiple experiments, we found that in order to capture more small vessel pixels while suppressing noise, it is a suitable choice that k is set to be 3. The proposed grayvoting algorithm parameters are thus set as follows: In the complementary image I c , the gray-voting algorithm parameters were set as follow to obtain more small vessel fragments: m = 11; k = −5; L ¼ M max ðm 2 À1Þ ; N ¼ À M max ðm 2 À1Þ ; where M max is the maximum gray level of the m × m window. Fig 7. shows the effects of the different values of k (e.g., -10, -5, -3, 0). In the vessel complementation process, T seed is the threshold of the seed fragment's number of pixels. In this study, we set T seed to be a constant value of 30 and T fragment to be a constant value of 100.

Dataset
The two publicly available databases (S1 File), DRIVE   [16] and STARE (Hoover et at., 2000) [17] were used to test the proposed methods. The DRIVE dataset contains 40 images that were obtained from a diabetic retinopathy screening program in The Netherlands. In the database, 33 images do not show any sign of diabetic retinopathy, and 7 show signs of mild early diabetic retinopathy. These 40 images have been randomly selected from the screening population, which consists of 400 diabetic subjects between the ages of 25 and 90 reported by Staal et al. (2004). Each image consists of 584 x 565 pixels. The STARE dataset contains 20 retinal fundus images, which consist of 605 x 700 pixels. Both datasets contain manual segmentation results. For the STARE dataset, we used two sets of manual segmentation results to evaluate the proposed algorithm. In the first manual segmentation dataset, 10.4% of the pixels were marked as vessels; while 14.9% of the pixels were marked as vessels in the second manual segmentation dataset, which contains more small vessel details than the first one.

Algorithm evaluation
We evaluated the proposed algorithm's accuracy, sensitivity and specificity of segmentation results. These evaluation measures are widely used in the vessel segmentation field. The primary concept of the evaluation method is to count the number of pixels that are true positives (TP), which describes the number of pixels that the algorithm segmented as vessel correctly; false positives (FP), which describes the number of pixels that the algorithm segmented as vessels incorrectly; true negatives (TN), which describes the number of pixels that the algorithm segmented as background pixels correctly; and false negatives (FN), which describes the number of pixels that the algorithm segmented as background incorrectly. These values can be obtained by comparing the algorithm's segmentation results with the "gold-standard" manual segmentation results. The evaluation method is defined by the formulae (18)(19)(20): Sensitivity  Fig 9(A) shows the original RGB image, and Fig 9B and 9C show the first and the second manually segmented results. The second manually segmented result shown in Fig 9(C) contained more small vessel details than Fig 9(B). Fig 9D and 9E) are the vessel segmented results by Soares et al. [9] and Zhang et al. [18]. Fig 9(F) shows the result of the proposed algorithm.
Tables 1 and 2 compare the results obtained using the proposed algorithm with those obtained by the other known algorithms with DRIVE and STARE datasets. It is shown that the proposed algorithm has a higher sensitivity than most of the other algorithms when using the  Gray-Voting and GMM to Segment Retinal Vessels DRIVE dataset; the sensitivity of the proposed algorithm is only lower than You's method [11] and Fraz's method [19]), the proposed algorithm also has the highest sensitivity using the STARE dataset.
Compared to the other methods, the results of the proposed algorithm contains more small vessel details, which agree with the fact that the proposed algorithm achieved a higher sensitivity when all methods segment the main vessels relatively well. It is also clear that the accuracy and specificity of the proposed algorithm is below the average levels shown in Tables 1 and 2. There are two primary reasons for this.
First, with the improvement of the retinal vessel segmentation algorithm, the problem about small vessel segmentation which ever to be a further step of vessel segmentation turn out to be realizable. Compared to a main vessel, a small vessel has fewer pixels and a lower vessel/background contrast; thus more sensitive filters are required for small vessel segmentation, which will increase the over-segmentation rate. Over-segmentation will decrease the signal <javascript:void(0);> to <javascript:void(0);> noise <javascript:void(0);> ratio <javascript: void(0);> (SNR <javascript:void(0);>) of the small vessel segmentation results, which lowers the accuracy and specificity of the proposed method.
Second, the evaluation results could be different according to different "gold-standard." With the STARE dataset, there are two manually segmented results. The first manually  segmented results, which are widely used for vessel segment evaluation, ignore many small vessels. Thus, because this approach contains significant amounts of small vessel information, the accuracy and specificity of the proposed method will be lower than the other vessel segment methods when using the first manually segmented result dataset to evaluate this approach. To develop this concept, we mark the TP, FP and FN in green, blue and red, respectively.  (D) and Fig 10(F) show the colored images using the first and the second manually segmented results as "gold standard," respectively. Compared to Fig 10(D), the blue region (i.e., overlapping area) shrinks clearly in Fig 10(F). Fig 10(F) shows that the proposed algorithm can segment some small vessels correctly that do not exist in the first manually segmented results but do exist in the second manually segmented results. In this case, some correctly segmented vessel pixels are incorrectly identified as over-segmentation pixels when the first manually segmentation results are used as the "gold standard." Table 3 shows the comparison of the other two methods with the proposed algorithm using the second manually segmented results in the STARE dataset as the "gold standard." The table shows that the proposed algorithm achieved much higher sensitivity than the other methods, and the highest accuracy; note that this accuracy is lower than those shown in Table 2. Fig 11 shows the accuracy and sensitivity of different algorithms with the two manually segmented results from 20 STARE dataset images. Compared to the other vessel segment methods, our approach achieved a substantially high sensitivity (Fig 11C and 11D)) and a relatively stable accuracy (Fig 11A and 11D).
Small retinal vessel segmentation plays an important role in clinical diagnosis; however, as mentioned above, small retinal vessel segmentation may lead to over-segmentation. Thus new methods to solve this problem must be developed. Moreover, the evaluation algorithm using accuracy, sensitivity and specificity is based on the pixel, yet the overlap rate may not indicate the true topological structure which may be of more important than the pixel overlap rate. For example, if the vessels' structures are all segmented perfectly in topological structures, but the vessels' widths are all thinner than the manually segmented image, the evaluation methods used in this study cannot accurately reflect the real vessel segment affection.
Compared with the other vessel segment methods, the proposed approach showed better segmentation results for both main and small vessels with relatively stable accuracies and high sensitivities.