Comparative Study of Retinal Vessel Segmentation Based on Global Thresholding Techniques

Due to noise from uneven contrast and illumination during acquisition process of retinal fundus images, the use of efficient preprocessing techniques is highly desirable to produce good retinal vessel segmentation results. This paper develops and compares the performance of different vessel segmentation techniques based on global thresholding using phase congruency and contrast limited adaptive histogram equalization (CLAHE) for the preprocessing of the retinal images. The results obtained show that the combination of preprocessing technique, global thresholding, and postprocessing techniques must be carefully chosen to achieve a good segmentation performance.


Introduction
Diabetic retinopathy (DR) accounts for about five percent of the causes of blindness globally, representing almost five million blind as stated by World Health Organization [1]. An early detection of DR is ensured through the regular examination of retinal images in diabetic patients, thus reducing the incidence of blindness cases. Automatic vessel segmentation has a great potential to assist ophthalmologists in the early detection of DR [2].
There have been various works done on the segmentation of vessels in retinal images. These works can be classified into two major categories. The first category is the unsupervised methods. This comprises vessel tracking [3][4][5], matched filter responses [6][7][8], morphology-based techniques [9,10], and locally adaptive thresholding [11]. The second category is the supervised methods. This category requires manually labeled images for training. This includes the use of neural networks [12], Bayesian classifier [13], k-nearest neighbor classifier [14], and SVM classifier [10,15], for the classification of the image pixels as either blood vessel or background tissue pixels. The method proposed in this paper belongs to the unsupervised method.
Chaudhuri et al. [6] implemented a matched filter by initially approximating the intensity of gray level profiles of the cross sections of retinal vessels using a Gaussian shaped curve. An Otsu thresholding technique was further applied to the matched filter response image to segment the retinal vessels. Hoover et al. [8] segmented retinal vessels by applying a threshold probing technique combining local vessel attributes with region-based attributes on matched filter response (MFR) image. Compared to [6] where a basic thresholding of an MFR was used, the method proposed by [8] reduced the false positive rate by as much as 15 times. Zhang et al. [16], having identified that the general matched filter responds to both vessels edges and nonvessel edges, extended the general matched filter with the first-order derivative of the Gaussian properties of the retinal vessels. Martinez-Perez et al. [17] applied the combination of scale space analysis and region growing to segment the vasculature. The technique proposed in [17] was, however, unable to segment the thin vessels. Zana and Klein [18] implemented a vessel segmentation method based on the use of mathematical morphology. Although the result achieved in [18] was good, the vascular structures were not always connected to one another. Jiang and Mojon [11] implemented an adaptive local thresholding based on a verification-based multithreshold probing scheme. The proposed technique in [11] was, however, faced with the limitations of some unconnected vascular structures and the inability to detect the thinner vessels.  [19] for vessel segmentation. A 7D vector composed of gray level and moment invariants-based features for pixel representation was computed, while a neural network classifier is used for the vessel segmentation. Soares et al. [13] generated a feature vector computed from the measurements at different scales of two-dimensional (2D) Gabor wavelet transform on each pixel. Bayesian classifier with Gaussian mixtures was further used to classify the resulting feature space as either a vessel or nonvessel pixel. Staal et al. [14] implemented a ridge-based vessel segmentation method. The retinal image ridges which cooccur approximately with vessel centre-lines were extracted. Primitives in the form of line elements were further composed of the ridges. The feature vectors computed for every pixel were classified using a k-nearest neighbour classifier and sequential forward feature selection. Niemeijer et al. [20] implemented vessel segmentation method based on pixel classification. Each pixel of the green plane of the retinal image was used to construct a feature vector. Consequently, these feature vectors were trained using a kNN-classifier. A filtered output and the pixel values within a neighborhood were compared. The best results were obtained from the filter output. Niemeijer et al. [20] further did the comparative study of the proposed vessel segmentation technique with the techniques proposed in [11,14]. Fraz et al. [21,22] implemented a supervised segmentation technique based on ensemble classifier of bootstrapped decision trees for the segmentation retinal vessel network. Lupaşcu et al. [23] implemented a supervised segmentation technique for detecting vessels using Ada-Boost classifier. A feature vector comprising local and spatial properties of the vessels were generated from the responses of various filters (matched filters, Gabor wavelet transform, and Gaussian filter and its derivatives). Ada-Boost classier was further trained and used to classify each pixel as either vessel or nonvessel. Ricci and Perfetti [24] proposed two different automated vessel segmentations based on line operators. The best of the two segmentation methods constructed feature vector for supervised classification using a support vector machine.
Szpak and Tapamo [25] used gradient based approach and level set technique. The proposed technique in [25] was, however, unable to detect the thinner vessels. Vlachos and Dermatas [26] implemented a multiscale line-tracking combined with a morphological postprocessing technique. Wang et al. [27] proposed multiwavelet kernels and multiscale hierarchical decomposition for the segmentation of retinal vessels. Mendonça and Campilho [28] combined differential filters for center-line extraction with morphological operators for the detection of retinal vessel network. Xiao et al. [29] proposed a Bayesian method with spatial constraint with level set for the segmentation of retinal vessels. Yin et al. [30] implemented a probabilistic tracking-based method for vessel segmentation.
Lupascu and Tegolo [31,32] trained a self-organizing map (SOM) on retinal images. The map was further divided into two classes using -means clustering [31] and modified fuzzy -means [32] techniques. The entire image is fed into SOM again and the class of the best matching unit on SOM is assigned to each pixel. A postprocessed technique based on hill climbing strategy on connected components was used to detect the vessel network. Saffarzadeh et al. [33] implemented a preprocessing phase based on -means followed by the use of multiscale line operators for the detection of retinal vessel network. With the help of k-means, the visibility of the vessels was enhanced and the impact of bright lesions was reduced. The retinal vessels were finally detected using the line detection operator in three scales.
Setiawan et al. [34] used contrast limited adaptive histogram equalization (CLAHE) to enhance the green channel color retinal image in order to enhance color retinal fundus image. The enhancement was achieved using histogram manipulation to get the uniform distribution of the intensity of the green channel. Contrast limited adaptive histogram equalization spreads the intensity distribution and adjusts the intensity of the original image. The red, green, and blue channels were finally combined as an enhanced color retinal image. Phase congruency on the other hand is a technique that is not affected by uneven illumination and contrast of the retinal image. A bank of log-Gabor filters was used by Kovesi to compute the phase congruency of an image and a binary segmentation was obtained by universal thresholding [35,36].
Amin and Hong [37] implemented the detection of retinal blood vessels using phase congruency at an high speed. Although the technique performed well in terms of speed, there is a need for a higher accuracy rate and a dynamically computed thresholding approach. Tagore et al. [38] used phase congruency to improve the contrast of vessel segments against the retinal background. A hierarchical clustering based histogram thresholding was then used to segment the contrast enhanced vessels. In related development, vessels cross-sectional profiles in the Fourier domain were represented and characterize using phase congruence by Zhu [39]. A bank of Gabor filter was used to transform the input image. The performance of the proposed technique in [39] was only described using visual results.
Although global thresholding technique has been used in [6], it has, however, been said to be inefficient for the retinal vascular segmentation [8,11]. This might have also resulted from certain limitations of the preprocessing phase [16]. In order to effectively produce good vessel segmentation, there is a need for an efficient preprocessing phase to enhance the vessels, good global thresholding technique, and an efficient postprocessing technique. This paper presents a study on the use different global thresholding techniques combined with different preprocessing and postprocessing techniques. The rest of this paper is organized as follows. Section 2 describes the methods and techniques used in this study. Section 3 explains the experimental setup, results, and discussion, while the conclusion is summarized in Section 4.

Methods and Techniques
Retinal fundus images are often characterized by noise due to illumination and contrast variation. Due to this, the use of global thresholding techniques for the detection of vessels in these noisy retinal images becomes challenging. In order to solve this problem, the need for an efficient preprocessing Computational and Mathematical Methods in Medicine 3 (1) Input the colored retinal image (2) Compute the gray scale of the colored retinal image (3) Image preprocessing through CLAHE using (1) and (2) (4) Enhance image using any of the following smoothing filters: (i) Gaussian filter (ii) Average filter (iii) Adaptive filter (iv) Average and Gaussian filters (5) Segment image using any of the following thresholding techniques: (i) Otsu threshold using (9) (ii) ISODATA threshold using (10) (6) Perform background subtraction (7) Perform morphological operation for noise removal (8) Obtain segmented vascular network Algorithm 1: Algorithm for CLAHE global-based thresholding technique.
technique is highly desirable. This section describes the two different preprocessing techniques and the different filtering techniques which are used to enhance the vessels. The different thresholding techniques and postprocessing techniques used in this paper are also described in this section. For the purpose of simplification, we group these techniques into two major approaches, namely, CLAHE global-based thresholding approach and phase congruence global-based thresholding approach as described in Algorithms 1 and 3.
(1) Preprocessing Phase. The different techniques used in the preprocessing phase are described below.
(a) CLAHE: CLAHE algorithm is used for partitioning the image into contextual regions and it applies the histogram equalization to each one. Figure 1 shows the colored, the gray scale, and the green channel of the retinal fundus image. CLAHE computes the local histogram at each pixel of the retinal image and performs histogram clipping, histogram renormalization, and output pixel mapping to an intensity proportional to its rank within the histogram. Given that ℎ is the histogram bin and ( × ) is the contextual region, the rank for a pixel with intensity is computed as follows: where the clip limit determines the contrast enhancement limit and ∑ =0 min( , ℎ ) describes the rank in a clipped histogram. Since each region will have a different number of clipped pixels, it is, however, beneficial to redistribute the part of the histogram that exceeds the clip limit evenly among all histogram bins to normalize the ranks computed in different regions. This normalization is provided by ∑ =0 ((∑ =0 max(0, ℎ − ))/( × )), where ℎ is the histogram bin in the different region. The rank of intensity in at ( , ) is computed and scaled to produce a fractional rank , such that 0.0 ≤ ≤ 1.0.
The output intensity level out is then computed in some grey scale ranging between 1 and 2 as follows: (b) The phase congruence model proposed in [35,36] has been very promising in the detection of object boundary in the presence of noise. The green channel is enhanced using phase congruence to minimize retinal image noise due to nonuniform illumination and contrast. Phase congruency is computed as follows: where ( ) is the amplitude, ( ) is the phase, and | ( )| is the local energy, given that In order to apply phase congruence to images, (3) is modified to be as follows: where ( , ) is the position of the pixel in the green channel of the retinal image, while and are the given scale and orientation, respectively. is the weighing factor for the distributed frequency, while estimates the image noise. The energy is computed using , Δ , ( , ), while is added to the denominator such that the divisor will be nonzero. The visual results of both CLAHE and phase congruence preprocessing techniques can be seen in Figure 2.
(c) Filters: the resulting images from CLAHE preprocessing technique are still affected to some extent by noise. In order to further enhance the retinal images, different filters are considered. The different filters considered are adaptive filter, average filter, and Gaussian filter. The combination of (1) Input the colored retinal image (2) Compute the gray scale of the colored retinal image (3) Compute the GLCM using (11) (4) Construct the GLCM-IDM feature matrix using (13) and (14) (5) Compute the IDM-based threshold using (16) Algorithm 2: Algorithm for computing IDM-based threshold.
(1) Input the colored retinal image (2) Extract the green channel of the colored retinal image (3) Perform preprocessing on image using phase congruence as indicated in (5) (4) Further enhance the preprocessed image using average filter (5) Segment image using any of the following thresholding techniques: (i) Otsu threshold using (9) (ii) ISODATA threshold using (10) (iii) IDM-based threshold computed in Algorithm 2 (6) Perform morphological operation for noise removal (7) Subtract image mask to obtain segmented vascular network Algorithm 3: Algorithm for phase congruence global-based thresholding technique. average filter and Gaussian filter was also used to further enhance the output of CLAHE preprocessing technique. Each of these different filtering approaches was considered in order to investigate their suitability for further enhancement of the retinal image. In related development, the resulting images from phase congruence were also enhance using average filter. The performance of each of the filtering approaches was, however, measured after the final vessel segmentation. The visual results from DRIVE database can be seen in (2) Global Thresholding. Automatic thresholding is potentially useful to dynamically select an optimal gray level threshold value for the segmentation of retinal vessels in the image from the background tissue based on their intensity distribution. The different global thresholding techniques studied in this paper are as follows.
(a) Otsu thresholding: global thresholding technique based on Otsu [40] is used on the results computed from phase congruence and CLAHE with filters for the initial estimation of the vessel network. The threshold that minimizes the intraclass variance as a weighted sum of variances of the two classes is explored in Otsu's method. The weighted sum of variances of two classes is expressed as follows: such that weights describe the probabilities of the two classes separated by a threshold and 2 variances of the classes. The class probability 1 ( ) is then computed from the histogram as follows: while the class mean 1 ( ) is computed as follows: such that ( ) is the value at the center of the th histogram bin. 2 ( ) and ( ) can likewise be computed on the histogram for bins greater than . Otsu further showed that minimizing the intraclass variance is the same as maximizing interclass variance; thus the desired threshold 2 ( ) is given as follows: where 1 ( ) and 2 ( ) are the means of the first and second group, respectively. The visual result of Otsu threshold on the image obtained from CLAHE preprocessing technique can be seen in Figure 3.
(b) ISODATA threshold selection: ISODATA threshold technique divides the histogram of the image output from phase congruence method into two using an initial threshold value 0 . The threshold is computed as follows: where 1 and 2 are the mean values of the two different parts of the histogram. This process continues until ℎ ≈ ℎ−1 .
(c) Inverse difference moment (IDM)-based binary thresholding: image signal statistics, particularly first-and second-order statistics, are good texture feature descriptors used for supervised segmentation techniques. Moments, first-order statistics, are concerned with individual image pixel properties while second-order statistics such as gray level cooccurrence matrix (GLCM) are concerned with individual pixel properties as well as the spatial interdependency  Computational and Mathematical Methods in Medicine of the two pixels at particular relative positions. The IDM texture information is computed using the GLCM of the gray scale of the retinal fundus image. The GLCM for the retinal fundus image is computed in the relative distance " " between the pixel pair and their relative orientation "Φ" across four directions (horizontal: 0 ∘ , diagonal: 45 ∘ , vertical: 90 ∘ , and antidiagonal: 135 ∘ ) as where ( , ) = means is the gray level of the pixel ( , ), and is defined as The IDM feature across the different distances " " and varying relative orientation "Φ" is defined as follows: where ( , ) is the ( , )th entry in a normalized gray scale spatial dependence matrix ( , ) / .
A multiscale IDM feature measurement across the varying distance " " and relative orientation "Φ" is used in the computation of an IDM feature matrix as follows: where 4 . The range measure of is given below as where 1 ≤ ≤ 4 and Φ is a row vector containing the range of each column of matrix ( ). The threshold value that will be used for the binarization of the output image from the phase congruence and average filter is computed as follows: (3) Postprocessing Phase. The different techniques used in the postprocessing phase are described below.  (a) Median filtering and morphological opening: median filter is used to restore the connectivity of several vessel lines by revealing some hidden pixels that belong to vessel lines. It is also used to get rid of the remaining noisy pixels. The choice of applying a 2 × 2 median filter has a good performance. This is referred to as (MO) in Tables 4 and 5. This is followed by the use of morphological opening in removing part of the remaining noisy pixels. The use of morphological opening 8 Computational and Mathematical Methods in Medicine Number of standard deviations of the noise energy beyond the mean at which we set the noise threshold point.

3
The fractional measure of frequency spread below which phase congruency values get penalized. cutOff 0.5 0.5 Controls the sharpness of the transition in the sigmoid function used to weight phase congruency for frequency spread.   (MO) alone and the combination of morphological opening and median filter was used. This is referred to as (MOMF) in Tables 4 and 5.
(b) Morphological directional filtering and reconstruction: the morphological directional filtering described in [12] is used to handle the several misclassifications that still remained. Morphological openings with line structuring elements orientation in five various directions, namely, 0, 30, 60, 120, and 150 degrees, are used. This paper adopts length of 1 pixel to keep vessel like structures with length of greater or less than 1. A logical OR for the responses of the five different directions and morphological reconstruction were performed on the image to remove a few erroneous regions before producing the final vessels network. The (MOMF) described in (a) above is combined with morphological directional filtering and morphological reconstruction for the purpose of performance investigation. This is referred to as (ATC) in Tables 4 and 5.   proposed method was evaluated using the retinal images on the publicly available DRIVE [41] and STARE databases [8]. DRIVE database is made up of 40 images captured with the use of Canon CR5 camera with 24-bit gray scale resolution and a spatial resolution of 565 × 584 pixels. The 40 images were divided into two groups. The first group of the DRIVE images is a training set made up of twenty images. The second group is a testing set made up of twenty images. DRIVE database also provides gold standard images as the ground truth for vessel segmentation for the comparative performance evaluation of different vessel segmentation algorithms. STARE database on the other hand consists of retinal images captured with the use of TopCon TRV-50 fundus camera with 24-bit gray scale resolution and spatial resolution of 700 × 605 pixels. The database provides 20 coloured retinal images and 20 hand-labeled images as the ground truth for the comparative performance evaluation of different vessel segmentation algorithms. The outcome of retinal vessel segmentation is a pixelbased classification result. Each pixel is either classified as vessel or background. Different events such as true positive (TP), true negative (TN), false positive (FP), and false negative (FN) take place during the pixel classification. An event is said to be TP if a pixel is correctly segmented as a vessel and TN when a pixel is correctly segmented as background. In related development, an event is said to be FN if a vessel pixel is segmented to be a background and a FP when a background pixel is segmented as a pixel in the vessel. The statistical performance measures commonly used for the evaluation of segmentation techniques are sensitivity, specificity, and accuracy. Sensitivity measure indicates the ability of a segmentation technique to detect the vessel pixels while specificity measure indicates the ability of a segmentation technique to detect background pixels. The accuracy measure, however, indicates the degree of conformity of the segmented retinal image to the ground truth. The measures are described in the equation below as Sensitivity = TP (TP + FN) ,

Experimental Results and Discussions
, where TP is true positive, TN is true negative, FP is false positive, and FN is false negative. Table 1 gives an overview of the parameter description and optimum parameter values of the phase congruence technique. In related development, different optimal values were empirically selected for the parameters used in CLAHE global thresholding approaches as described in Table 2. Figure 4 shows retinal images and their segmentation results obtained through phase congruence using different global-based thresholding techniques on DRIVE database. Figure 6 also shows the segmentation result obtained from a diseased retinal from DRIVE database using phase congruence combined with IDM thresholding technique. Figure 12 shows the result of phase congruence-based global thresholding approach on STARE database. Figure 5 shows different segmentation results obtained through CLAHE combined with different filters using Otsu thresholding technique while Figure 7 shows different segmentation results obtained through CLAHE combined with different filters using ISODATA thresholding technique on DRIVE database. Figure 11 also shows the result of CLAHEbased global thresholding approaches on STARE Database. Figure 8 describes the average sensitivities, specificities, and accuracies of the segmentation results obtained from phase congruence-based global thresholding approaches    [6] 0.8773 0.3357 0.9794 Staal et al. [14] 0.9442 0.7345 0.9773 Niemeijer et al. [20] 0.9416 0.7145 0.9801 Zana and Klein [18] 0.9377 0.6971 0.9769 Jiang and Mojon [11] 0.9212 0.6399 0.9625 Marín et al. [19] 0.9452 N/A N/A Ricci and Perfetti [24] 0.9646 N/A N/A Martinez-Perez et al. [17] 0.9181 0.6389 0.9496 Soares et al. [13] 0.9466 N/A N/A Vlachos and Dermatas [26] 0.9285 0.7468 0.9551 Akram and Khan [42] 0.9469 N/A N/A Amin and Hong [37] 0.9191 0.6608 N/A Mendonça and Campilho [28] 0.9463 0.7315 N/A Tagore et al. [38] 0 The best results achieved using phase congruence-based global thresholding approaches are obtained using IDMbased thresholding compared to ISODATA and Otsu thresholding. Phase congruence combined with IDM-based thresholding has very good accuracies due to the accurate segmented vessels and possesses good sensitivities due to the ability to segment some thin vessels. It is, however, still unable to segment the thinnest vessels. The best average accuracy of 0.94302 and average sensitivity of 0.71520 are achieved using morphological opening postprocessing technique. The next in rank of phase congruence-based global thresholding approaches is IDM-based thresholding combined with all postprocessing techniques combined giving average accuracy and sensitivity results of 0.93772 and 0.73910. IDMbased thresholding combined with morphological opening combined with median filter gives average accuracy and sensitivity results of 0.93596 and 0.74247. Phase congruence combined with IDM-based thresholding generally performed better than all the CLAHEbased global thresholding approaches. Phase congruence combined with IDM-based thresholding also gives better performance compared to Otsu and ISODATA thresholding combined with phase congruence. The performances of Otsu and ISODATA thresholding coupled with phase congruence are, however, at their best when morphological opening and median filter are combine for the post processing phase.
Tables 4 and 5 describe the performances of the best techniques from the different segmentation methods investigated in this paper and other previously published works using DRIVE and STARE databases.
Phase congruence combined with IDM-based thresholding using morphological opening postprocessing technique presents a higher average accuracy rates on DRIVE and STARE databases compared to the previously proposed phase congruence based technique by Amin and Hong [37]. Tagore et al. [38] achieves a lower average accuracy rate on DRIVE database but a higher average accuracy rate compared to the best phase congruence-based global thresholding approach presented in this paper. The technique proposed in [38], however, did not present the sensitivity and specificity measures. In related development, the phase congruence technique proposed by Zhu [39] discussed only the visual performance. Tables 4 and 5 also compare the results obtained in this paper with other results achieved in other literatures.

Conclusion
The performance of different vessel segmentation approaches based on combination of different preprocessing techniques, global thresholding, and postprocessing techniques has been investigated. It has also been shown that the combination of preprocessing technique, global thresholding, and postprocessing techniques must be carefully chosen to achieve a good segmentation performance. It is, however, important to state that the paper shows that sensitivity, specificity, and accuracy measures must all be high to ascertain a good segmentation performance. It was also shown that phase congruence combined with IDM-based thresholding generally performs better compared to phase congruence combined with ISODATA and Otsu threshold. Phase congruence combined with IDM-based thresholding is at its best on DRIVE database but did not have a better performance compared to the best of the CLAHE-based global thresholding approaches on STARE database. CLAHE-based global thresholding approaches were, however, shown to have maintained high accuracy rates across DRIVE and STARE databases. Although good accuracy and specificity rates were achieved, the sensitivity rate shows that global thresholding approach is still limited at efficiently segmenting the thin vessels. Our future work shall investigate the use of more robust segmentation techniques for the detection of both large and thin vessels in retinal images.