AUTOMATIC RETINAL VESSEL DETECTION AND TORTUOSITY MEASUREMENT

As retinopathies continue to be major causes of visual loss and blindness worldwide, early detection and management of these diseases will help achieve significant reduction of blindness cases. However, an efficient automatic retinal vessel segmentation approach remains a challenge. Since efficient vessel network detection is a very important step needed in ophthalmology for reliable retinal vessel characterization, this paper presents the study on the combination of difference image and k-means clustering for the segmentation of retinal vessels. Stationary points in the vessel center-lines are used to model the detection of twists in the vessel segments. The combination of arc-chord ratio with stationary points is used to compute tortuosity index. Experimental results show that the proposed k-means combined with difference image achieved a robust segmentation of retinal vessels. A maximum average accuracy of 0.9556 and a maximum average sensitivity of 0.7581 were achieved on DRIVE database while a maximum average accuracy of 0.9509 and a maximum average sensitivity of 0.7666 were achieved on STARE database. When compared with the previously proposed techniques on DRIVE and STARE databases, the proposed technique yields higher mean sensitivity and mean accuracy rates in the same range of very good specificity. In a related development, a non-normalized tortuosity index that combined distance metric and the vessel twist frequency proposed in this paper also achieved a strong correlation of 0.80 with the expert ground truth.


INTRODUCTION
Diabetic retinopathy (DR), an eye disease secondary to diabetes, is a major cause of visual loss and blindness worldwide (WHO; 2010).The focus in the evaluation on DR is usually on hard exudates, microaneurysms and vessel morphology.Retinopathy of prematurity (ROP) has also become a major cause of blindness to children in many middle income countries (O'Sullivan et al.;1997;Varughese et al.;2008;Chen et al.;2008;Gergely and Gerinec;2009;Hakeem et al.;2012;Adio et al.;2014).
In the face of the global prevalence of DR and ROP, the cases of vision loss and blindness tend to increase in the absence of efficient diagnosis and management approaches of the diseases (Rechtman et al.;2007).The families, communities and countries affected by these epidemics are also likely to suffer serious economic setbacks caused by the financial burden, reduced-earnings and reduced-productivity due to visual impairment and blindness (Atlas;2013).Ophthalmologists, with the help of detected vessel network, focus on retinal vessel feature analysis during the diagnosis of these diseases.Manual detection and analysis of the retinal vessels has however been a very tedious and time consuming task that requires trained and skilled personnel who are often scarce (Sussman et al.;1982;Verma et al.;2002).
Cost-effective interventions to manage diabetes can prevent the economic setbacks arising from complication such as blindness (Klonoff and Schwartz;2000).Regular retinal examinations for diabetic patients can help in an early detection of DR and significant reduction of the cases of blindness.Automatic segmentation and analysis of the retinal blood vessels are very important for the automatic detection of DR.Given the fact that the screening of large population of patients can be tedious and time consuming since there is a shortage of specialists, automated approach is cost-effective as it has the potential to increase the productivity of the ophthalmologists in such situations (Abràmoff et al.;2010a,b).Mobile-based diagnostic systems can also be used by field workers attending to larger number of people with diabetes in the rural areas where there are no ophthalmologists (Prasanna et al.;2013).This will yield a cost-effective and efficient early management of DR as patients who need further medical care can be referred to ophthalmologists at the early stage rather than advanced stage of the DR.
Automated vessel segmentation has however been a challenging problem due to complexities such as noise due to nonhomogeneous illumination caused by different reasons such as; curved surface of the retina, the degree of dilation is highly variable across patients, involuntary movements of the patients eye due to the bright flash-light during the acquisition of the retinal fundus images (Joshi and Sivaswamy;2008;Li et al.;2012).Although existing methods have made great progress in this field, it remains the subject of on-going research as there is need for further improvement in the detection of both large and thin vessels.Supervised vessel segmentation techniques are computationally expensive since training time is required (Condurache and Aach;2006;Lupascu et al.;2010).Another major drawback of the supervised vessel segmentation techniques is that their performance is highly dependent on the labelled training sample by experts which could sometimes be extremely difficult, expensive or unavailable (Li and Li;2013).Furthermore, a new set of labelled training sample is often required for retraining to obtain good segmentation performance when detecting vessels on a new set of retinal images (Staal et al.;2004;Marín et al.;2011).
This paper presents an investigatory study on the combination of difference image and k-means clustering for retinal vessel network segmentation.A study on the characterization of the detected vessels is also presented.
The rest of this paper is organized as follows.Section two describes several works previously done on the detection of retinal vessels and vessel tortuosity.Section three describes the methods and techniques used for the detection and tortuosity measurement of vessels in this paper.Section four explains the experimental setup, results and discussion, while the conclusion is drawn in section five.
Several retinal vessel segmentation and tortuosity measurement techniques have been proposed in the literature.This section gives a detailed review of segmentation and tortuosity measurement of vessels in retinal images.
Different segmentation methods have been previously proposed for the segmentation of retinal vessels (see (Fraz et al.;2012b) for extensive review up until 2012).These techniques can be divided into two major approaches namely unsupervised and supervised methods.
In supervised vessel segmentation methods, different algorithms (Nekovei and Sun;1995;Sinthanayothin et al.;1999;Niemeijer et al.;2004;Staal et al.;2004;Soares et al.;2006;Ricci and Perfetti;2007;Osareh and Shadgar;2009;Lupascu et al.;2010;Marín et al.;2011;You et al.;2011;Fraz et al.;2012c, 2013) are used for learning the set of rules required for the retinal vessel extraction.A set of manually segmented retinal vessels by experts are considered as the reference images.These reference images are used for the training phase of the supervised segmentation techniques.The availability of reliable labelled training samples by experts for the supervised segmentation methods could sometimes be expensive or unavailable.Fraz et al. (2012cFraz et al. ( , 2013) ) implemented a supervised segmentation technique based on ensemble classifier of bootstrapped decision trees for the segmentation retinal vessel network.Lupascu et al. (2010) implemented a supervised segmentation technique for detecting vessels using Ada-Boost classifier.A feature vector comprising of local and spatial properties of the vessels was generated from the responses of various filters (matched filters, Gabor wavelet transform, Gaussian filter and its derivatives).Ada-Boost classifier was further trained and used to classify each pixel as either vessel or non-vessel.The drawback of this method is the high computational time of 125 seconds required to segment vessels in each retinal image.Another drawback of this method is its inability to detect the thin vessels as retinal vessels width vary from 15 pixels in very large vessels to 3 pixels in small vessels across the fundus images (Li et al.;2012).Marín et al. (2011) generated a 7-D vector composed of gray-level and moment invariantsbased features for pixel representation, while a multilayer feed forward artificial neural network (ANN) classifier was used for the vessel segmentation.Nekovei and Sun (1995) used a back-propagation ANN for the segmentation of blood vessels in angiography.ANN was trained using ground truth images of manually labelled angiograms.The ANN was finally applied to segment the blood vessels in the angiograms.Sinthanayothin et al. (1999) combined principal component analysis (PCA) with ANN for the localization of anatomical structures in retinal images.Niemeijer et al. (2004) implemented vessel segmentation method based on pixel classification.Each pixel of the green plane of the retinal image and responses of Gaussian matched filter were used to construct feature vectors.Consequently, these feature vectors were classified using a kNN classifier.Osareh and Shadgar (2009) applied multi-scale Gabor filters, and PCA for the extraction of features in the retinal image.The features generated were consequently used to classify the pixels of the retinal images as either vessels or non-vessels using Gaussian mixture model (GMM) and support vector machines (SVM).Ricci and Perfetti (2007) proposed automated vessel segmentation based on line operators.Two segmentation methods were considered.One of the segmentation methods used two orthogonal line detectors with the gray level of the target pixel to construct a feature vector for supervised classification using a support vector machine.Although the technique had a good performance, there are however high false detection around the border of the optic disc.The need for a new set of training samples and the retraining of classifier before applying it on a new dataset also remains the drawback of the method.Soares et al. (2006) implemented the combination of two-dimensional (2-D) Gabor wavelet transform and Bayesian classifier for retinal vessel segmentation.A feature vector comprising a multi-scale 2-D Gabor wavelet transform responses and pixel intensity was generated from the retinal images.Each pixel in the retinal image was further classified as vessel or background tissue using a Bayesian classifier.Although the technique had a good performance, segmentation of thinner vessels as well as false detection around the border of the optic disc remains a challenge.The need for a retraining of classifier before applying it on a new dataset also remains the limitation of this method.Staal et al. (2004) utilized ridge information and kNN classifier with sequential forward feature selection for retinal vessel segmentation.One of the drawbacks of the technique in (Staal et al.;2004) is its inability to segment thinner vessels.This method is also expensive as it requires a time consuming manual labelling (about 2 hours per image) and a retraining of the classifier before applying it on a new dataset.The need for a retraining of classifier before applying it on a new dataset also remains the limitation of this method.You et al. (2011) combined radial projection with SVM using a semisupervised self-training approach for the segmentation of vessels.Radial projections were used to locate the vessel centre-lines and the low contrast blood vessels.Having enhanced the vessels using modified steerable wavelet, feature vector was generated using line strength measures.SVM classifier was further used in a semi-supervised self training manner for the classification of each pixel as either vessel or nonvessel.A supervised classification approach utilizing multi-scales vesselness and texture features was proposed by Poletti and Ruggeri (2012) to segment the retinal vessel network in wide-field images acquired with RetCam fundus camera.Poletti and Grisan (2014) proposed a supervised technique based on ADABoost classifier using an array of optimal discriminative convolution kernels for detection of vessels in retinal images.
The approaches based on unsupervised classification, on the other hand, attempt to find inherent patterns of blood vessels in retinal images that can be further used to determine that a particular pixel belongs to vessel or not.The training data or hand labelled ground truths are not needed for the design of algorithm in these approaches.These approaches are matched filter-based (Chaudhuri et al.;1989;Hoover et al.;2000;Chanwimaluang and Fan;2003;Cinsdikici and Aydın;2009;Zhang et al.;2010;Chakraborti et al.;2014), scale space-based (Martínez-Pérez et al.;1999;Grisan et al.;2004;Vlachos and Dermatas;2010;Poletti et al.;2011;Wang et al.;2013), tracking-based (Yin et al.;2013), model-based (Szpak and Tapamo;2008;Xiao et al.;2013), adaptive thresholding-based (Jiang and Mojon;2003;Cornforth et al.;2005;Akram and Khan;2013;Mapayi et al.;2014), mathematical morphology-based (Zana and Klein;2001;Mendonca and Campilho;2006;Jiménez et al.;2010;Fraz et al.;2011;Miri and Mahloojifar;2011;Fraz et al.;2012a) and clustering-based (Tolias and Panas;1998;Yang et al.;2008;Kande et al.;2010;Lupas ¸cu and Tegolo;2011a,b;Sun et al.;2011;Saffarzadeh et al.;2014).Chaudhuri et al. (1989) implemented matched filter response (MFR) 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.The proposed technique however had a poor detection of the vessels.Hoover et al. (2000) segmented retinal vessels by applying a threshold probing technique combining local vessel attributes with region-based attributes on MFR image.While compared to the technique used in (Chaudhuri et al.;1989) where a basic thresholding of an MFR was used, the method proposed in (Hoover et al.;2000) reduced the false positive rate by as much as 15 times.Chanwimaluang and Fan (2003) proposed the combination of matched filter and entropy for the segmentation of retinal vessels.The performance measure of the proposed technique was only visual.Cinsdikici and Aydın (2009) implemented the combination of matched filter and ANT colony algorithm for the detection of retinal vessels.Zhang et al. (2010), having identified that the general matched filter responds to both vessels edges and the nonvessel edges, extended the general matched filter with the first-order derivative of the Gaussian properties of the retinal vessels.The improved technique was, however, faced with the inability to segment the thinner vessels.Chakraborti et al. (2014) implemented an unsupervised segmentation technique that combines vesselness filter and matched filter using orientation histogram for the segmentation of retinal vessels.
Martínez-Pérez et al. (1999) used a combination of scale space analysis and region growing to segment the vessel network.Vlachos and Dermatas (2010) proposed a retinal vessel segmentation method based on multi-scale line-tracking procedure and morphological post-processing.Wang et al. (2013) proposed multi-wavelet kernels and multi-scale hierarchical decomposition.Vessels were enhanced using matched filtering with multi-wavelet kernels.Yin et al. (2013) implemented a probabilistic trackingbased method for vessel segmentation.A Bayesian method with maximum a posteriori (MAP) was used for detecting the retinal vessel edge points.Li et al. (2006) combined multi-scale analysis based on Gabor filters, scale multiplication, and region based thresholding to achieve adaptive thresholding for vessel segmentation.Szpak and Tapamo (2008) used gradient based approach and level set technique for the segmentation of retinal vessels.The proposed technique was however unable to detect the thinner vessels.Xiao et al. (2013) proposed a Bayesian method with spatial constraint for the segmentation of retinal vessels.The spatial dependence of the posterior probability of each pixel in relation to their neighboring pixels was utilized.An energy function was further defined and a modified level set approach was used for the vessel segmentation.Jiang and Mojon (2003) implemented an adaptive local thresholding based on a verification-based multi-threshold probing scheme for the detection of vessel network.Vessels' properties such as contrast, curvilinear angle, width, and size were modelled into the verification phase, and a number of thresholds were used to probe and detect vessel network.The combination of the resulting detected vessel networks obtained from probed thresholds followed by post processing techniques was used to generate the final segmented vessel network.The technique was however faced with the limitations of some unconnected vascular structures and the inability to detect the thinner vessels.Akram and Khan (2013) enhanced the vascular pattern using 2-D Gabor wavelet and followed by a multilayered thresholding technique that applied different threshold values iteratively to generate gray level segmented image.Cornforth et al. (2005) applied wavelet analysis, supervised classifier probabilities and adaptive threshold procedures, as well as morphology-based techniques.A comparative study of global thresholding techniques previously implemented by Mapayi et al. (2015) showed that global thresholding approaches are limited at efficiently segmenting thin retinal vessels.In another study, Mapayi et al. (2014) implemented an adaptive thresholding technique utilizing different types of local homogeneity information for retinal vessel segmentation.
Zana and Klein ( 2001) implemented a vessel segmentation method based on the use of mathematical morphology and cross-curvature evaluation.Mendonca and Campilho (2006) combined differential filters for center-line extraction with morphological operators for filling vessel segments considering intensity and morphological properties.Miri and Mahloojifar (2011) found ridges by applying multi-structure elements to enhanced retinal image.Morphological operators by reconstruction was further applied to get rid of the ridges outside the vessel tree.Fraz et al. (2011Fraz et al. ( , 2012a) ) (2004).Poletti et al. (2011) proposed an automatic extraction of vessel centerline in wide-field ROP retinal images using a sparse tracking scheme.An hybrid approach comprising of a fuzzy clustering and mathematical morphology was proposed by Yang et al. (2008).A morphological top-hat operation was used for retinal image smoothing and background information removal.The vessels were then extracted through fuzzy clustering.
Retinal vessel tortuosity is one of the early symptoms of some retinopathies (Wasan et al.;1995) and are of great importance in the diagnosis of ROP (Davitt and Wallace;2009).Retinal vessel network tortuosity is the measure of twists and curvature of a vessel.Retinal vessel tortuosity measures are used to determine the state of retinal images as either healthy or diseased even when there are no visible pathologies as shown in Fig. 1.One of the first changes in vessels morphology to occur in DR patients is the increase in vessel tortuosity (Scheie;1953;Kristinsson et al.;1997;Owen et al.;2008;Taarnhøj et al.;2008;Dougherty et al.;2010;Onkaew et al.;2011;Chakravarty and Sivaswamy;2013).Tortuosity has been shown to be a more reliable vascular parameter in differentiating ROP severity than vessel width (Capowski et al.;1995).High blood pressure has been linked to diabetes due to its relationship with hypertension (Bedell;1945;Kylstra et al.;1986;Patasius et al.;2007) and tortuosity has also been identified as an indicator of high blood pressure.Degeneration of retinal vessel walls and changes in their elastic properties have been identified to be major causes of vessel tortuosity (Azegrouz et al.;2006).The duration of diabetic disease has a strong effect on changes in vessel tortuosity (Lim et al.;2010;Sasongko et al.;2010).
Although vessel tortuosity are visually analyzed and determined by ophthalmologists based on the curvature and twisting rate of vessels, qualitative grading of vessel tortuosity suffers from inter-observer and intra-observer variations (Wolffsohn et al.;2001).Several methods have been applied to quantitatively determine the measure of tortuosity of a vessel (Bullitt et al.;2003).The term often called the tortuosity index (TI), is used for the estimation of the vessel tortuosity.TI provides a reproducible measure of vessel characteristics.
Vessel tortuosity was computed as the ratio of the actual vessel length to the length of the underlying chord (Lotmar et al.;1979;Hart et al.;1999;Heneghan et al.;2002;Martinez-Perez et al.;2002;Wilson et al.;2008).This technique assumes a vessel to be non-tortuous if it is straight line and tortuous while the radius of curvature is longer than the chord length of the vessel.The techniques however failed to differentiate varying vessels with the same length but with different tortuosity.Hart et al. (1997) implemented automatic tortuosity measurement using integral curvature methods.Dougherty and Varro (2000) used second derivatives along central axis of the blood vessels for automatic tortuosity measurement.Bullitt et al. (2003) applied the number of inflection points for the automatic vessel tortuosity measure.Grisan et al. (2003Grisan et al. ( , 2008) ) implemented automatic tortuosity measures using arc to chord ratio and the points of changing curvature sign.The vessel extraction and inflection point placement was however done manually (Sinthanayothin et al.;2010;Onkaew et al.;2011).Capowski et al. (1995) addressed tortuosity measurement using the number of direction changes in vessel course.Azegrouz et al. (2006) proposed a vessel-thickness tortuosity measurement method utilizing the average of the two curvature values computed at two boundary points of the vessel.
Generally, automatic tortuosity measurement methods based on local tortuosity measurement performs better while compared to the techniques based on global tortuosities (Wilson et al.;2008).However, the contribution of each local tortuosity measure is considered for computing the global tortuosity measure.
Although, much has been achieved in the previous works, the performances obtained from literature suggest the need for further work to address the robust segmentation of large and thinner retinal vessels.There is also need to investigate an efficient approach for vessel tortuosity measurement.

MATERIALS AND METHODS
Since efficient vessel network detection is a very important step needed in ophthalmology for reliable retinal vessel characterization, a robust segmentation technique that performs the detection of large and thin vessels in a timely manner is highly needed.This section presents an unsupervised vessel segmentation approach to address the problems of inability to detect the thinner vessels (Martínez-Pérez et al.;1999;Vlachos and Dermatas;2010;Saffarzadeh et al.;2014;Mapayi et al.;2015) connectivity loss in vessel network (Jiang and Mojon;2003).Since unsupervised vessel segmentation methods are not dependent on labelled training set, the method investigated in this paper also overcomes the major drawback of supervised segmentation methods which is their high dependence on the labelled training (Li and Li;2013) with the retraining of the classifiers (Soares et al.;2006;Marín et al.;2011).The detailed description of the materials with the vessel segmentation methods and tortuosity measurement techniques investigated in this paper are further presented in sections Materials, Vessel Segmentation and Vessel Tortousity Measurement respectively.

MATERIALS
Experiments were carried out using MATLAB 2010a on an Intel Core i5 2410M CPU, 2.30 GHz, 4GB of RAM.The proposed method was evaluated using the retinal images on the publicly available databases (DRIVE;2004;STARE;2000). DRIVE ( 2004) is a publicly available database consisting of 40 colour retinal fundus images.These images were obtained from a diabetic retinopathy screening program in the Netherlands.The images were acquired using a Canon CR5 non-mydriatic 3CCD camera with a 45 degree field of view (FOV).Each image was captured using a 24-bit RGB colour at 768 × 584 pixels.The FOV of each image is circular with a diameter of approximately 540 pixels.This set of forty colour retinal fundus images are divided into two groups.The first group of the images is made up of twenty images which are a training set for supervised segmentation techniques.The training set are, however, not used in this paper because we applied an unsupervised segmentation method.The second group is a testing set made up of twenty images.All the human observers that manually segmented the vessels in the retinal images were instructed and trained by an experienced ophthalmologist.A set of manual segmentations of the vessel network is available for the training images.Two manual segmentations sets X and Y are available for the test cases.In set X, 577,649 pixels were marked as vessel and 3,960,494 pixels as background (12.7% vessel).In set Y , 556,532 pixels were marked as vessel and 3,981,611 as background (12.3% vessel).It was, however, observed that it took a human observer an average time of 7200 seconds (2 hours) to segment the vessels in each of the retinal images.Set X of the manual segmentations of the test cases is used as gold standard (ground truth), while set Y is often compared to the performance of the automatic segmentation techniques on the database.STARE (STructured Analysis of the Retina) STARE ( 2000) is a publicly available database that consists of 20 retinal images captured with the use of TopCon TRV-50 fundus camera with 24-bit grey-scale resolution and spatial resolution of 700×605 pixels.The approximate diameter of the FOV is 650 pixels.Two observers manually segmented all the images.The first observer segmented 10.4% of pixels as vessel, against 14.9% vessels for the second observer.The segmentations of the two observers are fairly different in that the second observer segmented many more of the thinner vessels than the first one.The 20 manually segmented images provided by the first observer is used as the ground truth for the comparative performance evaluation of different vessel segmentation algorithms in the literature.
It is however important to note that one of the major reason why the first and second observer disagreed in the manual segmentations of the retinal vessels was because of the difficulty to distinguish vessel pixels due to JPEG-artefacts.The decisions of the human observers at distinguishing pixels at vessel borders as either vessel or non-vessel are also subjective.Although the second observer labelled more pixels as vessels when compared to the first observer on the STARE database, the first observer took a more conservative view of the boundaries of vessels and in the identification of small vessels when compared to the second observer.Hence, the manual segmentation of the first observers of DRIVE and STARE databases are used as the ground truth for the comparative performance evaluation of different vessel segmentation algorithms in the literature.
In a related development, 50 digital fundus images containing 50 different vessel segments used for the vessel tortuosity measurement.Expert's ground truth that determined the state of each vessel segment as either tortuous or normal was also obtained from a consultant ophthalmologist.

VESSEL SEGMENTATION
The green component of the coloured retinal image is used for segmentation since it provides the best vessel-background contrast (Staal et al.;2004).A detailed description of the vessel segmentation approach investigated in this paper is given as: 1. Extraction of the green channel of the coloured fundus image.
2. Process the green channel of the retinal image using different filtering techniques.
3. Generate the difference image.
4. Segment the retinal vessels from the generated difference image using k-means clustering.
5. Implement a post-processing phase that combines median filter and morphological opening for the removal of misclassification.

Filtering Techniques
The green channel of the retinal image is enhanced using different filtering techniques.Linear filters such as mean filter and Gaussian filter are used for smoothing images.Although these filters reduce image noise, they are however weak at preserving edges in an image.Non-linear filter, particularly median filter, is efficient at removing image noise as well as preserving edge information in images.It is however worth noting that the selected filter window sizes should not be too large in order to efficiently manage the noise due to illumination variation conditions of the retinal image.It is also important to carefully select window sizes that have sufficient data points for good enhancement.The window sizes of the filters applied in this work are empirically selected.For the purpose of investigation in this paper, mean and Gaussian filters are considered for the convolution while median filter is used for the filtering of the retinal image.The image padding does not affect the retinal image as there are black background pixels after the circular field of view.We replicated the border pixel values (which are zeros) for the image padding at the border of the retinal image (i.e including the background), and this has no effect as the background pixels around the border are also black.

Difference Image
A difference image is generated by subtracting the green channel of the coloured retinal image from the convolved of filtered retinal image to balance the illumination of the retinal image.The difference image D(x,y) is given below as: where U is the convolved or filtered retinal image, V is the green channel of the retinal image.D(x,y)= {D ρ (x, y), D υ (x, y), D σ (x, y)}, where D ρ (x, y) is the difference image based on median filter (DIMDF), D υ (x, y) is the difference image based on mean filter (DIMNF) and D σ (x, y) is the difference image based on Gaussian filter (DIGF).A model that combines two possible difference images was also investigated.The combinations obtained are:

K-Means Clustering
K-means clustering is an unsupervised segmentation technique used in defining the natural group of pixels in an image.This is achieved by classifying input image data points into different classes through a set of distances computed using the image data points and centroids.K-means clustering technique is used for dividing n sample input data points of X = {x 1 , x 2 , ...., x n } into a group of k clusters.This is achieved by considering the similarities among the input points within the same cluster as well as the differences among the different clusters.Sum of squared errors is a very useful criterion measure for clustering.Given k clusters, the sum of squared errors is computed as: such that µ i is the centroid of the i th cluster S i , {i = 1, 2 . . .k.}.

POST-PROCESSING
Median filter and morphological opening are used for the post-processing phase to remove noisy pixels appearing to be broken vessels and restore the connectivity of the background tissues (black background) while the several vessel networks are still preserved.Median filtering is a non-linear smoothing operation often used in image processing to reduce noise and preserve edges at the same time.In this research, the median filtering is performed by moving a 2 × 2 sliding-window through all the pixels of the binarized retinal image containing detected vessels.When the chosen mask is greater than 2 × 2, the median filter categorizes some of the smaller retinal vessels as noise.This is followed by morphological opening.In the application of morphological opening to the filtered image, erosion is applied to remove the remaining noisy pixels and dilation is subsequently applied.The morphological opening with line structuring elements orientated in five different directions namely 0 • , 30 • , 60 • , 120 • , 150 • was applied.Another mophological operator applied is the area opening.This technique removes the remaining noisy pixels appearing to be broken vessels using the connectivity value 8.After this, the FOV mask is subtracted from the segmented image to obtain the detected vessel in the circular field of view.Given that the filtered segmented image is S filtered im , the final output image after the application of morphological opening is computed as: where γ(S filtered im ) is the final segmented vessel network after post-processing and ε(S im ) describes the application of erosion and δ (S filtered im ), the application of dilation.

VESSEL TORTOUSITY MEASUREMENT
Retinal vessel network tortuosity is the measure of twists and curvature of a vessel (Grisan et al.;2003, 2008;Onkaew et al.;2011;Chakravarty and Sivaswamy;2013).Several clinical studies have described the relation between vessel tortuosity and retinopathies such as DR, ROP (Scheie;1953;Capowski et al.;1995).Fig. 5 shows the variation in the directional changes, frequency of twists and the differences in the vessel lengths of the different vessels a, b, c, d, e, f and g.

Description of Tortuosity
With the help of ophthalmologic rules and clinical definition, the metric that provides numeric index for vessel tortuosity measurement is utilized (Bullitt et al.;2003).Rather than inflection points alone (Bullitt et al.;2003), the stationary points of the vessel are considered to be very important for computing the tortuosity measure in this paper.A combination of arcchord ratio with stationary points are used to compute tortuosity index using the local tortuosity measurement approach.Fig. 6 shows a sample vessel of interest is marked, detected vessel of interest and skeletonized form of detected vessel of interest.
The tortuosity measure applied in this paper utilizes the chord length, arc length and the frequency of vessel curve using the stationary points as described in Algorithm 1.The blood vessel is skeletonised to extract the center line of the vessel.
The chord length L chord of the skeletonised blood vessel is defined as where L chord is the length of the straightline connecting the two end-points of the skeletonised blood vessel.
The arc length L arc of the skeletonised blood vessel is defined as where L arc is the actual length of the skeletonised blood vessel.
The pixel locations of skeletonised blood vessels are extracted to plot a vessel description graph as shown in Fig. 7.A set of frequent smaller swings occurred on the skeletonised vessel segment as a result of noise.In order to efficiently detect the number the stationary points, a moving average filter with span n = 5 was applied to smoothen the smaller swings that occurred on the skeletonised vessels to remove the noise.The stationary points of the vessel description graph are computed using the gradients of the vessel description graph.We applied a 1-point difference scheme.With such scheme, each gradient is computed as either 0, 1 or -1.With the help of the 1-point difference, a change of direction in the vessels can also be detected.The stationary points are points on the vessel description graph where the gradient is zero.There are however three major types of stationary points: minimums, maximums and inflection points as illustrated in Fig. 8 and Fig. 9 (Thomas;1996).Given the skeletonised blood vessel T with vessel coordinates X and Y , the gradients of the vessel description graph are computed as: such that X = {x 1 , x 2 , . . ., x n }, Y = {y 1 , y 2 , . . ., y n }, G = {g 1 , g 2 , . . ., g n } and −1 ≤ G ≤ 1 where GεI.
The nature of each type of the stationary points is shown in Fig. 9     The stationary points in the vessel description graph are computed using the vector of gradient obtained from Eq. ( 9) and the nature table in Table 1.In Algorithm 1, the gradient patterns D i of the stationary points detectors S1, S2, S3, S4 are computed.The vessel twist counter is initialize to zero and updated using UpdateRealVesselPattern as the presence of twists are detected in the vessel by matching the detector patterns D i with vessel pattern E j obtained from the thinned segmented vessel.The total number of twists, SPC count , is then computed.We investigate a normalized TI metric TI freq1 that combines distance metric and the vessel twist frequency obtained using the stationary points.This is computed as: where SPC count is the twist frequency of the vessel curves.
Several literature have applied a distance metric by dividing the arc length by the chord length (see Kalitzeos et al.;2013).While the previously presented distance metric divided the arc length by the chord length (Kalitzeos et al.;2013), a different distance metric that divides the chord length and the number of stationary points by the arc length is used for the non-normalized TI in this study.The non-normalized TI is computed as

RESULTS
Empirically, we established that using window sizes 11 × 11, for mean and Gaussian filters, 15 × 15 for median filter were effective for the detection of the vessel network.The average time taken to segment the vessel networks in each of the retinal image on DRIVE and STARE databases ranges from 3.4 to 4.0 seconds.
The performance measures used are sensitivity, specificity and accuracy.The measures are described in the Eqs.12-14 below as: where TP is True Positive, TN is True Negative, FP is False Positive and FN is False Negative.An event is said to be TP when a pixel is rightly segmented as a vessel and TN when rightly segmented as background.Conversely, 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.Sensitivity measure indicates the ability of a segmentation technique to detect the pixels belonging to vessel while specificity measure indicates the ability of a segmentation technique to detect background pixels.Therefore, the higher the sensitivity rate, the higher the rate of vessel detection.The accuracy measure indicates the degree of conformity of the segmented retinal image to the ground truth.(1999) 0.9181 0.6389 0.9496 Zana and Klein (2001) 0.9377 0.6971 0.9769 Jiang and Mojon (2003) 0.9212 0.6399 0.9625 Niemeijer et al. (2004) 0.9416 0.7145 0.9801 Staal et al. (2004) 0.9442 0.7345 0.9773 Mendonca and Campilho (2006)  A receiver operating characteristic (ROC) curve performance measure is a plot of the rightly classified pixels, referred to as true positive rate (TPR) versus the fraction of the wrongly classified pixels as vessels, referred to as false positive rate (FPR).Area under the curve (AUC) is a performance measure computed from the ROC curve.The ROC are determined by setting the threshold on the posterior probability.The image with k = 10 in Fig. 3 serve as the posterior probabilities.
The differences between the various results achieved by the best methods investigated in this paper are assessed by means of z-test statistics (Foody;2009).The sensitivity is the proportion of the rightly detected vessel pixel while the accuracy rate is the proportion of the rightly detected pixel in the retinal image.To evaluate the variability of the performance of the different methods on the same dataset, we used the same set of data point samples.The significant difference in average sensitivity rate indicates the significant difference in the detection rate of the vessels and the significant difference in accuracy rate indicates the significant difference in the accuracy rate of the retinal vessel segmentation.Both measures are very important as good vessel segmentation method for efficient vessel analysis in ophthalmology requires that sensitivity and accuracy rates of the segmentation be good (Mapayi et al.;2015).The statistical significance of a difference between two proportions is evaluated as where p = (x 1 + x 0 )/(n 1 + n 0 ) with x 0 and x 1 representing the number of cases of the rightly detected pixel in the classifications of sample datapoints of size n 0 and n 1 respectively.The statistically difference is computed at an asymptotic confidence level of 95% (α = 0.05).K-means with DIMDF has the best average accuracy and average sensitivity rates on DRIVE and STARE databases.DIMDF combined with k-means clustering detected majority of large and thin vessels, while very few thinner vessels remain undetected.In a related development, DIMDMNF and DIMDGF while combined with k-means yielded higher average sensitivity rates but lesser average accuracy rates while compared to DIMDF combined with k-means clustering technique.The increase in sensitivity is as a result of the integration of the other difference image with DIMDF.DIMDF however has the highest average accuracy rate while compared to DIMDMNF and DIMDGF while combined with k-means clustering technique.Fig. 11 gives a visual description that compares the results obtained from vessels segmented using DIMDF DIMDMNF and DIMDGF while combined with k-means clustering technique on DRIVE database.Fig. 10 and Fig. 12 show the visual results obtained from vessels segmented using DIMDF DIMDMNF, DIMNGF and DIMDGF while combined with k-means clustering technique in comparison with the ground truth and the result obtained by Hoover et al. (2000) on STARE database.As shown in Figs. 10, 11 and 12, majority of large and thin vessels are detected, a very few thinner vessels still remain undetected.The false detection around the border of the optic disc are higher on the results obtained by using DIMDMNF and DIMDGF while combined with k-means clustering but lesser on the results obtained while using k-means combined with DIMDF as shown Fig. 11.This is also responsible for the higher average accuracy rate achieved by DIMDF while combined with k-means clustering technique.The result obtained while DIMNGF is combined with k-means clustering technique in Fig. 12 however shows that the thinner vessels are not detected.
The performance measures of the different previously proposed techniques on DRIVE and STARE databases are also compared to the proposed techniques in Tables 2 and 3.The average accuracy rate is written by Avg.Acc., average sensitivity rate as Avg.Sens. and average specificity rate as Avg.Spec.

DISCUSSION
This section presents a detailed discussion of the results obtained in this study.

COMPARISON WITH EXISTING SEGMENTATION METHODS ON DRIVE DATABASE
The k-means clustering techniques investigated achieved average sensitivity rates of 0.7079, 0.7399, 0.7518 and 0.7581 with corresponding average accuracy rates of 0.9523, 0.9556, 0.9531 and 0.9516 respectively on DRIVE database.All the techniques that combined difference image with k-means clustering presented significantly higher average sensitivity rates of 0.7079, 0.7399, 0.7518 and 0.7581 with significantly higher average accuracy rates of 0.9523, 0.9556, 0.9531 and 0.9516 when compared with the average sensitivity rates of 0.3357, 0.6389, 0.6399, 0.6971, 0.6562, 0.6565 and 0.6522 with their corresponding average accuracy rates of 0.8773, 0.9181, 0.9212, 0.9377, 0.9459, 0.9482 and 0.9267 achieved by Chaudhuri et al. (1989), Martínez-Pérez et al. (1999), Jiang and Mojon (2003), Zana and Klein (2001), Lupas ¸cu and Tegolo (2011a), Lupas ¸cu average accuracy rates of 0.9480 and 0.9483 when compared with the average accuracy rates of 0.9509, 0.9500 and 0.9492 obtained from three of the k-means clustering segmentation techniques with significant differences.The human observer STARE (2000) presented a higher average sensitivity of 0.88949 when compared with all the average sensitivity rates of k-means clustering segmentation techniques.The higher average sensitivity rate achieved by the human observer STARE (2000) showed the higher detection rate of both large and thin vessels by the human observer in the retinal images.The human observer STARE (2000), however, presented a significantly lower average accuracy rate of 0.9354 when compared with the average accuracy rates of 0.9509, 0.9500 and 0.9492 obtained from three of the k-means clustering segmentation techniques.

COMPARISON OF TORTUOSITY MEASURES WITH EXPERT'S GROUND TRUTH
The tortuosity measures computed in the previous section utilized chord length, arc length and the frequency of twists using the stationary points of the vessel curves.
The non-normalized tortuosity index (TI) combined a different distance metric and the vessel twist frequency as described in Eq. 11.The normalized TI as described in Eq. 10, utilized the chord length for the normalization previously presented distance metric combined with the twist frequency measure.A threshold value 1.0 is applied to differentiate tortuous from non-tortuous retinal vessels.Hence, TI values less than 1.0 are considered non-tortuous while TI values greater than or equal to 1.0 are considered to tortuous.Table 5 shows the correlation of the tortuosity measures obtained with the expert ground truth using Spearman's rank correlation coefficient.The result obtained shows that the non-normalised metric (TI freq2 ) has a stronger correlation with the expert ground truth with significant differences when compared with the other techniques investigated.
In conclusion, this paper presents unsupervised vessel segmentation approach that combines difference image and k-means clustering for the segmentation of retinal vessels.K-means clustering combined with difference images addressed the major drawback of supervised segmentation methods which is their high dependence on the labelled training with the retraining of the classifiers.K-means clustering combined with difference image based on median filter also addressed the segmentation of large and thinner retinal vessels as well as the reduction of false detection around the border of the optic disc.It is also shown that an integration of difference images based on median filter with difference image based on mean filter and difference image based on Gaussian filter while combined with k-means clustering technique segmented both large and thin vessels in retinal images.The reason for the good performance of these segmentation techniques are due to the very good edge-preserving property of median filter.The vessel segmentation techniques that combined median filtered based difference image with k-means clustering technique present higher average sensitivity and average accuracy rates when compared with the previous techniques on DRIVE and STARE databases.The vessel segmentation techniques investigated in this paper are also computationally fast ranging from 3.4 to 4.0 seconds on DRIVE and STARE databases.
Having designed and implemented the detection of twists in vessel segments through the use of stationary points in the vessel center-lines, we also showed that the proposed non-normalized TI that combined the distance metric and the vessel twist frequency has a stronger correlation with the expert's ground truth when compared to the previously presented distance metric and the normalized TI investigated in this paper.
combined vessel centre-lines detection and morphological bit plane slicing for the detection of vessel network in retinal images.Jiménez et al. (2010) implemented vessel centre-lines detection and a combination of morphological operations for retinal vessel segmentation.Kande et al. (2010) combined matched filtering and a spatially weighted fuzzy C-means for vessel segmentation in retinal images.Lupas ¸cu and Tegolo (2011a,b) having trained a self-organizing map (SOM) on retinal images, proposed two clustering techniques for vessel segmentation.The map was divided into two classes using k-means clustering technique (Lupas ¸cu and Tegolo; 2011a) and modified Fuzzy C-Means clustering algorithm (Lupas ¸cu and Tegolo; 2011b).A post-processing technique based on hill climbing strategy on connected components was further applied to detect the vessel network (Lupas ¸cu and Tegolo; 2011a,b).Saffarzadeh et al. (2014) implemented a preprocessing phase based on k-means followed by the use of multi-scale line operators for the detection of retinal vessel network.With the use of k-means, the visibility of the vessels was enhanced and the impact of bright lesions reduced.The retinal vessels were finally detected using the line detection operator in three scales.Sun et al. (2011) implemented morphological multi-scale enhancement in combination with fuzzy filter and watershed transformation for vascular segmentation.Tolias and Panas (1998) implemented a fuzzy C-means algorithm to segment vessels in retinal angiogram images.The technique however did not segment thinner vessels due to their low contrast against background.An automatic extraction of the vascular structure in retinal images using a sparse tracking technique was also proposed in Grisan et al.
combination of median filter and mean filter based difference images (DIMDMNF), D υ σ is the combination of mean filter and Gaussian filter based difference images (DIMNGF) and D ρ σ is the combination of median filter and Gaussian filter based difference images (DIMDGF).The results obtained in Eqs.2-4 are normalized to the interval [0, 255].Fig. 2 shows the coloured retinal image, the green channel of the coloured retinal image and the difference image obtained using the DIMDF investigated in this paper.

Fig. 2 :
Fig. 2: Coloured retinal image on DRIVE database (a), green channel of the coloured retinal image (b), difference image based on DIMDF (c).

Fig. 4 :
Fig. 4: Detected vessels based on k-means and difference image before post-processing(a) k-means and DIMDF (b) k-means and DIMDGF (c) k-means and DIMDMNF (d) k-means and DIMNGF.

Fig. 6 :
Fig. 6: Coloured retinal with vessel of interest marked (a), detected vessel of interest (b), and skeletonization of detected vessel of interest (c).

Fig. 7 :
Fig. 7: Vessel graph describing X and Y coordinates of skeletonised blood vessels.

Fig. 8 :
Fig. 8: Graph describing different stationary points between point A & point B.

Fig. 11 :
Fig. 11: The red rectangular box indicates the noise rate around the border of the optic disc (a) DRIVE database coloured fundus image (b) DRIVE database gold standard (c) Segmented vessels using k-means with DIMDF (d) Segmented vessels using k-means with DIMNGF (e) Segmented vessels using k-means with DIMDGF (f) Segmented vessels using k-means with DIMDMNF.

Fig. 13 :
Fig. 13: ROC curves showing the performance of the various difference images combined with k-means clustering on DRIVE.

Fig. 14 :
Fig. 14: ROC curves showing the performance of the various difference images combined with k-means clustering on STARE.

Table 1 :
and Algorithm 1: Vessel Twist Detection Based on Stationary Points Algorithm Input: S1, S2, S3, S4 are the sets of stationary point detectors.T is the thinned segmented vessel.Output: SPC count is total number of vessel twists, D i is the set of gradient pattern obtained from the detectors and E j is the set of gradient pattern obtained from the thinned segmented vessel / * Computes and Updates the computed gradient patterns of the detectors based on stationary points.

Table 1 :
Stationary points nature table.

Table 2 :
Performance of different segmentation methods on DRIVE database.

Table 3 :
Performance of different segmentation methods on STARE database.

Table 4 :
Comparison of AUC of the proposed techniques with previous works on DRIVE and STARE databases.

Table 5 :
Correlation of tortuosity measures with expert's ground truth.