Retinal Vessel Segmentation via Structure Tensor Coloring and Anisotropy Enhancement

Retinal vessel segmentation is one of the preliminary tasks for developing diagnosis software systems related to various retinal diseases. In this study, a fully automated vessel segmentation system is proposed. Firstly, the vessels are enhanced using a Frangi Filter. Afterwards, Structure Tensor is applied to the response of the Frangi Filter and a 4-D tensor field is obtained. After decomposing the Eigenvalues of the tensor field, the anisotropy between the principal Eigenvalues are enhanced exponentially. Furthermore, this 4-D tensor field is converted to the 3-D space which is composed of energy, anisotropy and orientation and then a Contrast Limited Adaptive Histogram Equalization algorithm is applied to the energy space. Later, the obtained energy space is multiplied by the enhanced mean surface curvature of itself and the modified 3-D space is converted back to the 4-D tensor field. Lastly, the vessel segmentation is performed by using Otsu algorithm and tensor coloring method which is inspired by the ellipsoid tensor visualization technique. Finally, some post-processing techniques are applied to the segmentation result. In this study, the proposed method achieved mean sensitivity of 0.8123, 0.8126, 0.7246 and mean specificity of 0.9342, 0.9442, 0.9453 as well as mean accuracy of 0.9183, 0.9442, 0.9236 for DRIVE, STARE and CHASE_DB1 datasets, respectively. The mean execution time of this study is 6.104, 6.4525 and 18.8370 s for the aforementioned three datasets respectively.


Introduction
Retinal blood vessel segmentation and extraction of different features based on these vessel segments, such as tortuosity, width, length, angle and color are exploited for diagnosis and screening of numerous kinds of diseases such as diabetic retinopathy (DR), arteriosclerosis, hypertension, retinopathy of prematurity and choroidal neovascularization [1,2].The retinal vessel segmentation is also a useful tool for computer-assisted surgery, multimodal image registration and biometric person identification [2][3][4][5].Manual segmentation of retinal blood vessels is a time-consuming, subjective and boring task which must be performed by trained physicians.The medical community generally admits that the retinal vessel segmentation is one of the basic steps for the development of the retinal disease diagnosis systems [6,7].A great number of studies related to retinal vessel segmentation have been published previously and these studies can be categorized into 3 main classes as kernel, classifier and tracking-based studies [7].
Kernel-based studies convolve the image with the kernels having different orientations and sizes.Nonetheless, this type of studies needs long execution times when the kernels get larger and have to be applied with more than one orientation.In addition to these, the response of a specific kernel fits with only the vessels whose standard deviation resemble to Gaussian function of that kernel.Thus, the vessels which have unrelated profiles may not be detected by the kernel-based methods [6].In order to segment the retinal vessels, Chaudhuri et al. came up with a Gaussian function-based two-dimensional linear kernel [8].Hoover et al. performed vessel segmentation by merging local and region-based features of the retinal vessels via applying a threshold probing technique on the response image of a matched filter [9].Gang et al. detected retinal vessels by using a second order Gaussian filter whose amplitude is modified [10].Jiang and Mojon tried to detect the vessels via an adaptive local thresholding scheme using a verification-based multi-threshold probing framework [11].Al-Rawi et al. enhanced Chaudhuri et al.'s matched filter by using optimized parameters derived from an exhaustive search optimization step [8,12].Cinsdikici and Aydin offered a retinal vessel segmentation method which is a combination of the matched filter and the ant colony algorithm [13,14].Zhang et al. performed retinal vessel segmentation by an extended version of the matched filter which exploits the symmetric cross sectional property of the vessels [15].
Classifier-based studies extract a feature vector and classify each pixel as vessel or non-vessel.The supervised and unsupervised methods were applied by the classifier-based studies.Nekovei and Ying offered a backpropagation Artificial Neural Network (ANN) for retinal vessel segmentation in angiography images [16].Sinthanayothin et al. proposed principal component analysis (PCA) and ANN for the detection of the optic disc, fovea, and retinal vessels [17].Niemeijer et al. offered a featured a vector for each pixel by exploiting the information from the green channel of the image and using the first and second order derivatives of Gaussian function-based matched filters having various scales [18].Later, Niemeijer classified the pixels using the k-Nearest Neighbor (k-NN) algorithm [18,19].Soares et al. segmented the vessels via 2-D Gabor wavelet and supervised classification [20].Ricci and Perfetti classified the retinal vessel pixels using Support Vector Machine (SVM) and a feature vector composed by using line operators [21].Marin et al. proposed a retinal vessel segmentation method by classifying a 7-D feature vector derived from gray-level and moment invariant-based features via ANN [22].Tolias and Panas segmented vessels in the retinal angiogram images via a fuzzy C-means (FCM) clustering algorithm [23].Simo and de Ves segmented arteries, veins and the fovea in the retinal angiograms by using the Bayesian image analysis [24].Salem et al. offered a RAdius-based Clustering ALgorithm (RACAL) which uses a distance-based criterion in order to depict the distributions of the image pixels [25].Villalobos-Castaldi et al. segmented the vessels by using the texture information by merging the gray-level co-occurrence matrix (GLCM) and the local entropy data [26].
Tracking-based studies start vessel tracking from at least one seed point and aim to trace whole vascular pattern step by step.This seed point may be set manually by user or automatically via morphological operations [7].This type of technique may be disadvantageous because the seed point requirement may cause the study remain as semi-automatic and missing any bifurcation point may result to lose the entire branch [6].Chutatape et al. performed the segmentation by tracking the vessels using Gaussian and Kalman filters [27].Quek and Kirbas offered a wave propagation and traceback mechanism which labels each pixel's vesselness likelihood by using a dual-sigmoidal filter [28].Delibasis et al. proposed a vessel tracking algorithm which utilizes a parametric model of a vessel composed of a stripe [29].
This study can be classified in the kernel-based category because of the Frangi Vesselness Filter (FVF) and Structure Tensor (ST) steps.However, FVF requires only a few different kernels not for different orientations but only for different scales.More and more, ST algorithm has a lightweight time complexity since it does not execute iteratively for various rotations and scales.The novel parts of this algorithm are utilizing ST as a vessel extraction technique followed by combining anisotropy enhancement function and applying Contrast Limited Adaptive Histogram Equalization (CLAHE) and enhanced mean surface curvature (EMSC) multiplication to the energy space of the tensor field as well as exploiting the coloring part of the tensor visualization method as a new segmentation algorithm.Although a great number of methods were developed, new retinal vessel segmentation techniques are needed still [6].This study proposes a fast method which extracts the main vessel arcs with a high sensitivity and reasonable accuracy.This article is the extended version of a previously presented conference paper [30].The main contribution of this article is reducing the false positive segmentation results especially occurring on the main vessel arc by using EMSC multiplication and ST instead of improved structure tensor (IST).Apart from that, this study is tested not only on DRIVE dataset but also on STARE and CHASE_DB1 datasets and higher performance results are obtained.Some extra post-processing steps such as lesion removal and small hole filling are also added to this study.

Materials
The DRIVE dataset contains 40 images which are captured via a Canon CR5 non-mydriatic 3CCD camera with a 45-degree field of view (FOV) having an approximately diameter of 540 pixels.7 of the images have mild early DR.The manual ground truth segmentation and the mask image corresponding to the FOV are given for each image [31].
The STARE dataset contains 20 images which are captured using a TopCon TRV-50 fundus camera at 35 degree FOV having an approximate diameter of 650 × 500 pixels and 10 of these images include pathology [9].
The CHASE_DB1 dataset contains 28 images which are captured at 30 degree FOV with a resolution of 1280 × 960 pixels from 14 patients in Child Heart And Health Study in England [32].

Methods
The main flowchart of the algorithm is depicted in Figure 1.

Frangi Vesselness Filter
FVF is an algorithm to enhance vessels or tubular structures in medical images which may have different modalities.The Hessian matrix of the target image which is symbolized by H can be obtained via second order Gaussian derivative of image f at point x and scale σ as it is shown Equation ( 1).
The sorted λ 1 and λ 2 Eigenvalues of Hessian matrix of a 2-D image are shown in Equation (2).
Frangi, defined two new metrics for 2-D images as S and R B to measure the background noise in images and the digression from a blob-like structure, respectively.These two metrics are shown Equation (3).
The vesselness function V Frangi which uses the aforementioned metrics is shown in Equation ( 4) where β and c are real positive and user defined parameters [33].
FVF is applied to images in order to emphasize all the retinal vessels with respect to the background.The β and c are set as 0.5 and 15 respectively.The scale range is determined as between 0.01 and 7.01 increasing by 0.5 for each scale.The enhanced filter response of FVF is shown in Figure 2a.

Structure Tensor
ST is basically a covariance matrix which is composed of the partial derivatives of gradients across the edges [34,35].One of the most important benefits of ST is the coherence concept which is calculated via its Eigenvalues.ST has been used for low-level feature analysis such as feature extraction, edge and corner detection after it was introduced by Förstner as well as Harris and Stevens [32,35,36].In the literature, there are basically linear and nonlinear types of ST [37,38].Both types apply smoothing functions or kernels to the basic form of ST which is the covariance matrix of the gradients.Smoothing the ST via Gaussian is a common approach [37].However, the smoothing operation also blurs the edge information as a side effect.In this study, ST is used solely without applying any linear or nonlinear functions since the fact that it is observed the blurring effect highly increases the false positives in the case of retinal vessel segmentation.
ST is an array of (2,2) which is composed of x and y edge derivatives of image f which is an array of (m,n).ST is calculated for each pixel of the image as in Equation ( 5) and a tensor field which is an array of size (m,n,2,2) is obtained.
A frame of the result of the ST is visualized via tensor visualization technique using ellipsoids as shown in Figure 2b.The size of the ellipsoid visualization field is set as (m/3,n/3,2,2) in order to get a basic idea about the structure of the tensor field.The flow of the sticky ellipsoids depicts the tensors which reside on the vessels while the blob-like ellipsoids represent the non-vessel region.

Anisotropy Enhancement Using Principal Eigenvalue
The contrast between vessel and non-vessel tensors on the obtained tensor field via ST is proposed to be improved by enhancing the anisotropy between the principal Eigenvalues.The λ 1 and λ 2 of the tensor field is obtained as arrays of (m,n,1).It is seen that contrast between the Eigenvalues of each tensor is not high enough and an enhancement function Z is applied to the principal Eigenvalue λ 1 as it is shown in Equation (6).The values of ε and α are set experimentally as 0.001 and −50 respectively.The result of the anisotropy enhancement is visualized via ellipsoids again as it is shown in Figure 3a [39].
2.2.4.Applying CLAHE to Energy of Tensor Field At this step, the energy, anisotropy and orientation of the resulting tensor field is calculated for each pixel using the resulting anisotropy enhanced ST field (AEST) as in Equation (7).The, b and c symbols in Equation ( 7) define the elements of the AEST tensor field for each pixel.In order to increase the contrast of the vessels a little bit more, CLAHE algorithm is applied to the obtained energy matrix [40].The energy, anisotropy and orientation values for each pixel are calculated using the values of the matrix elements of the corresponding tensors as in Equation (7).The before and after effect of CLAHE over the energy matrix is shown in Figure 3b,c.

Enhanced Mean Surface Curvature Multiplication
An infinite number of curves can be drawn on a point of a continuous surface M which can be parametrized by u and v axes.The curvature value at a certain point of each curve C is calculated using the second derivative of its arc length parametrization which is symbolized by C (s).The normal curvature, K n , on a certain point of a curve C is calculated by projecting C (s) onto the normal vector N at that point as in Equation (8).The maximum and minimum normal curvature values for each point on a surface M are called as principal curvatures and shown as k 1 and k 2 as in Equation (9).The first fundamental form variables are shown as E, F and G as in Equation (10) whereas the second fundamental form variables are shown as e, f and g as in Equation ( 11) using partial derivatives of surface M. The continuous calculation of Gaussian and mean surface curvatures which are respectively shown as K and H are given as in Equation (12).The discrete calculation of K and H are shown as in Equation (13).The discrete recalculation of k 1 and k 2 is shown in Equation ( 14) [41].
At this step of the study, the CLAHE applied energy space is modeled as a surface in order to enhance the segmentation results by reducing false positive pixels via mean surface curvature information.Firstly, the discrete versions of K, H, k 1 and k 2 are calculated and then EMSC is calculated using γ which is experimentally set as 10 as in Equation (15).Then, it is converted to binary form via decomposing the positive and negative valued pixels which are respectively shown as p and q symbols as in Equation (16).Afterwards, the formerly CLAHE applied Energy space is multiplied by EMSC as in Equation (19).
∀p i,j , ∀q i,j ∈ EMSC, q i,j > 0, p i,j ≤ 0, EMSC q i,j ← 1, At the end, the histogram equalized and EMSC multiplied energy matrix is embedded back to the CLAHE applied AEST tensor field (CAEST) again via Eigen recomposition steps as shown in Equation ( 18) and Figure 3f [42].The a , b and c symbols in Equation ( 18) define the final values of the elements of the CAEST tensor field for each pixel.

Tensor Coloring
The 4-D tensor fields are able to be visualized using ellipsoids.The radius, orientation and color of the ellipsoids are configured with respect to the values of the tensors.In this study, the tensor coloring method of the visualization technique which uses ellipsoids is exploited as a segmentation method since it can be observed that this visualization technique has a reasonable ability to distinguish the vessel and non-vessel tensors.The tensor coloring procedure is described in Equation (19).Firstly, the maximum and minimum of the sum of diagonal tensor values of CAEST are computed.Then, a color map matrix of (256, 3) is generated in MATLAB as it is shown in Figure 4a.Afterwards, the matrix of color map value which is represented as cmv is calculated and its negative values are set as zero.Thereafter, the corresponding values of cmv with respect to each column of the Map matrix constitutes the Colors1, Colors2 and Colors3 matrixes as can be seen in Figure 4b-d, respectively.Finally, these three color matrixes are merged into Colors matrix using cat command of MATLAB as it is shown as a colorful RGB image in Figure 4e [42].
The Colors3 image is selected because of its higher contrast and OTSU method is applied to it in order to determine 3 main gray histogram distributions as it is shown in Figure 5a [43].Then, the vessels are thresholded and the remaining artifacts such as small unconnected pixels are removed using bwareaopen command of MATLAB.The final segmentation result is shown in Figure 5b,c.

Post-Processing
In this study, two different approaches are offered for the lesion removal procedure at the post-processing step.The effects of the post-processing procedures are shown on the segmentation steps of an image having DR in STARE dataset as in Figure 6.The first approach targets especially the spread hard and soft exudate lesions.The highest 10% of the histogram of the product of hue and green color space of the retinal image having these kinds of lesions is removed from the segmentation result as in Figure 6d.
The second approach aims to remove all the other lesions or residual artefacts by taking into account their geometric structures such as solidity and eccentricity.The ratio of the region of a connected component to its convex hull is defined as solidity.The ratio of the distance between the foci of the ellipse which has same the second-moments as a connected component and its major axis length is defined as eccentricity [44].It was empirically inspected with a high probability that the connected components having solidity value above 0.3 and eccentricity value below 0.95 are lesion.After removing the lesion-like connected components as in Figure 6e, all the residual connected components which are smaller than 100 pixels are deleted from the segmentation result.Finally, the small holes occurring on the vessels because of the central vessel reflex are filled using some successive morphological operations as shown in Figure 6f.

Evaluation and Results
Based on the evaluation criteria of the previous retinal vessel segmentation studies, the true positive (TP) is defined as the count of the pixels which are detected as vessel in both the ground truth and segmented images whereas the true negative (TN) is defined as the count of the pixels which are detected as non-vessel in both the ground truth and segmented images.The false positive (FP) is defined as the count of the pixels which are detected as vessel in segmented image but non-vessel in the ground truth image while the false negative (FN) is defined as the count of the pixels which are detected as non-vessel in segmented image but vessel in the ground truth image.The sensitivity (SN) is defined as the ratio of the TP to the sum of TP and FN.The specificity (SP) is defined as the ratio of the TN to the sum of TN and FP.The accuracy (ACC) is defined as the ratio of the sum of the TP and TN to the number of the pixels in the FOV [6].The area under a receiver operating characteristic curve (AUC) is commonly used for the measuring the performance of the classifiers for different thresholds which is not applicable for unsupervised methods.Another definition such as AUC = (Sensitivity + Specificity)/2 is more suitable for unsupervised methods like the proposed study [45,46].
In this study, the sensitivity, specificity, accuracy, AUC and execution time values for the 40, 20 and 28 images of DRIVE, STARE and CHASE_DB1 datasets are respectively shown in Tables 1-3.Some of segmentation results and their ground truth images are shown respectively for DRIVE, STARE and CHASE_DB1 datasets in Figures 7-9.

Discussion
FVF is utilized as a useful preprocessing step in terms of only being able to focus on vessel-like regions.The ST method which is used in this study is successful on the edge extraction by using gradient information.Anisotropy Enhancement step of the 4-D tensor field contributed to vessel enhancement by not deforming the topology of the vessels but only emphasizing the vessel borders and regions.Applying CLAHE to the energy space is useful for normalizing the pixel values after anisotropy enhancement step.The EMSC multiplication step reduces the number of the false positive pixels which reside on the vessel borders even if it also causes a lot of small holes to occur on the vessels which are filled during post-processing.In this study, the proposed lesion removal procedure is also exploited successfully in order to eliminate the miscellaneous artefacts exaggerated by the FVF and ST.
This study does not require parameter calibration except for the parameters of FVF, EMSC and lesion removal criteria as well as can be used automatically for various retinal images having light inhomogeneity and different resolution values.In addition to these, this study offers a new approach for segmentation literature by utilizing tensor visualization technique which uses ellipsoids.As far as we know, this study is the first one in which the tensor coloring method is used for segmentation.
The performance results of the proposed study and the recent state of the art studies in the last two years are compared as in Tables 4-6.This study suggests encouraging sensitivity, specificity, accuracy, AUC and execution time results with respect to the literature for DRIVE, STARE and CHASE_DB1 datasets.The performance results of the proposed study and the results of the other studies which are lower than the corresponding evaluation criterion is written as bold in the below tables.The sensitivity of this study is better than [45,[47][48][49][50][51][52] whereas its specificity is better than [53] and its accuracy and AUC is respectively more successful than [48] and [45,51] for DRIVE dataset as in Table 4.The sensitivity of this study is more successful than [45,49,51,52,54,55] whereas the performance of its specificity is better than [53,55] and its accuracy and AUC is respectively higher than [48,53,55] and [45,51] for STARE dataset as in Table 5.The performance of the sensitivity of this study is higher than [56] whereas its accuracy is more successful than [57] for CHASE_DB1 dataset as in Table 6.Performance comparison of the average execution times for DRIVE and STARE datasets between this study and some of the state of the art studies are given in Table 7 [52].The algorithm of this study is implemented in MATLAB and executed in a PC with 2.2 GHz Intel core i7 CPU and 4 GB RAM.
One of the superiorities of this study with respect to the previous studies is its higher sensitivity values by not trading off the specificity and accuracy values dramatically.The other advantage of this study is that, as far as we know, it has the lowest execution time in the literature as shown in Table 7 [52].The most important reason of the lower AUC values of this study is because of the fact that the previously mentioned calculation method of the AUC for the unsupervised algorithms yields lower values with respect to the supervised methods which can be crosschecked by looking sensitivity, accuracy and AUC values of the studies listed in Tables 4-6 [45,46].The proposed novel ST coloring-based vessel segmentation method can be utilized for different kind of micro tubular segmentation problems.This method can also be easily implemented on parallel software and hardware systems since the fact that the parallelizable structure of the employed tensor operations.The lesion removal algorithm which is proposed in this study can be exploited not only for the retinal image segmentation task but also for the different problems of medical image analysis.
This study can be beneficial for diagnosis of hypertension, retinopathy of prematurity, arteriosclerosis, choroidal neovascularization and diabetic retinopathy by providing the elegant sensitivity, specificity and accuracy metrics as well as execution speed [1,2].In future, we will optimize all the parameters of this proposed retinal vessel segmentation algorithm via hyper parameter optimization techniques which will definitely improve its accuracy and speed.

Conclusions
In this study, the retinal vessel segmentation was performed with mean sensitivity of 0.8123, mean specificity of 0.9342, mean accuracy of 0.9183, mean AUC of 0.8732 and mean execution time of 6.104 s on DRIVE dataset.The vessel segmentation was achieved with mean sensitivity of 0.8126, mean specificity of 0.9442, mean accuracy of 0.9312, mean AUC of 0.8784 and mean execution time of 6.4525 s on STARE dataset.
The vessel detection was performed with mean sensitivity of 0.7246, mean specificity of 0.9453, mean accuracy of 0.9236, mean AUC of 0.8349 and mean execution time of 18.8370 s on CHASE_DB1 datasets.The proposed algorithm provided encouraging results especially by detecting main vessel arcs with reasonable sensitivity, specificity and accuracy values for 3 different publically available datasets and offering a new approach for segmentation literature by exploiting the coloring algorithm of the tensor visualization technique which uses ellipsoids.

Figure 1 .
Figure 1.The flow chart of the proposed retinal segmentation algorithm.

Figure 2 .
Figure 2. Intermediary steps of the algorithm (a) result of FVF; (b) tensor visualization of the result of ST using ellipsoids. λ

Figure 3 .
Figure 3. Intermediary steps of the algorithm (a) the result of the anisotropy enhancement; (b) energy matrix before applying CLAHE; (c) energy matrix after applying CLAHE; (d) EMSC; (e) energy after multiplication of CLAHE result and EMSC; (f) revisualization of the obtained new tensor field.

Figure 5 .
Figure 5. Intermediary steps of the algorithm (a) image after OTSU; (b) image after post processing; (c) final segmentation result.

Figure 6 .
Figure 6.Post-processing procedures on image 2 on STARE dataset (a) retinal image having DR; (b) first segmentation result; (c) product of hue and green colors; (d) effect of histogram-based bright lesion removal step; (e) effect of solidity and eccentricity-based lesion removal step; (f) effect of small hole filling step.

Figure 7 .
Figure 7. Segmentation results of the algorithm on DRIVE dataset (a) segmentation of image 2; (b) segmentation of image 9; (c) segmentation of image 14; (d) ground truth of image 2; (e) ground truth of image 9; (f) ground truth of image 14.

Figure 8 .
Figure 8. Segmentation results of the algorithm on STARE dataset (a) segmentation of image 7; (b) segmentation of image 8; (c) segmentation of image 12; (d) ground truth of image 7; (e) ground truth of image 8; (f) ground truth of image 12.

Figure 9 .
Figure 9. Segmentation results of the algorithm on CHASE_DB1 dataset (a) segmentation of image 5; (b) segmentation of image 9; (c) segmentation of image 27; (d) ground truth of image 5; (e) ground truth of image 9; (f) ground truth of image 27.

Table 1 .
The Statistics of the Obtained Results on DRIVE Dataset.

Table 2 .
The Statistics of the Obtained Results on STARE Dataset.

Table 3 .
The Statistics of the Obtained Results on CHASE_DB1 Dataset.

Table 4 .
Performance Comparison between the Proposed Study and the Recent Studies on DRIVE.

Table 5 .
Performance Comparison between the Proposed Study and the Recent Studies on STARE.

Table 6 .
Performance Comparison between the Proposed Study and the Recent Studies on CHASE_DB1.

Table 7 .
Performance Analysis between the Average Execution Time of the Proposed Study and Some State of the Art Studies on DRIVE and STARE.