Retinal Blood Vessel Segmentation Using Hybrid Features and Multi-Layer Perceptron Neural Networks

Abstract

Segmentation of retinal blood vessels is the first step for several computer aided-diagnosis systems (CAD), not only for ocular disease diagnosis such as diabetic retinopathy (DR) but also of non-ocular disease, such as hypertension, stroke and cardiovascular diseases. In this paper, a supervised learning-based method, using a multi-layer perceptron neural network and carefully selected vector of features, is proposed. In particular, for each pixel of a retinal fundus image, we construct a 24-D feature vector, encoding information on the local intensity, morphology transformation, principal moments of phase congruency, Hessian, and difference of Gaussian values. A post-processing technique depending on mathematical morphological operators is used to optimise the segmentation. Moreover, the selected feature vector succeeded in outfitting the symmetric features that provided the final blood vessel probability as a binary map image. The proposed method is tested on three known datasets: Digital Retinal Image for Extraction (DRIVE), Structure Analysis of the Retina (STARE), and CHASED_DB1 datasets. The experimental results, both visual and quantitative, testify to the robustness of the proposed method. This proposed method achieved 0.9607, 0.7542, and 0.9843 in DRIVE, 0.9632, 0.7806, and 0.9825 on STARE, 0.9577, 0.7585 and 0.9846 in CHASE_DB1, with respectable accuracy, sensitivity, and specificity performance metrics. Furthermore, they testify that the method is superior to seven similar state-of-the-art methods.

1. Introduction

In this paper, an automatic DR diagnosis system is proposed. Pointing to the ‘fact sheet’ of the world health organisation (WHO), over 2.2 billion people all over the world have a vision impairment or blindness. One billion of these cases can be prevented or at least controlled if these cases are detected at an early stage. Diabetic retinopathy (DR) is the fifth leading cause of vision impairment and is estimated at 4.5 million out of the one billion cases [1,2,3,4]. Currently, the need for eye disease screening and check ups has increased rapidly. Despite that, there is a notable shortfall of ophthalmologists, particularly in developing countries. Hence, the automation detection of the DR process by supervised procedures is a benefit and can significantly assist in working out this issue [5,6].

In the last couple of decades, the clinical importance of the segmentation of retinal blood vessels from fundus images has attracted researchers to diagnose many diseases. Consequently, the segmentation of retinal blood vessels is an attractive field of research as it helps in retinal image analysis and computer-assisted diagnosis [7,8].

Retinal vessel investigation is noninvasive, as it can be observed directly without any surgical intervention [9]. Additionally, it is an indispensable element to any automatic diagnosis system for ophthalmic disorders. The variation in the structure of retinal blood vessels can serve as an indicator of serious diseases such as stroke, glaucoma, hypertension, arteriosclerosis, cardiovascular, age-related macular degeneration (AMD), refractive error and DR [7,10]. It is commonly accepted by the retinal medical community that uses $2D$ colour fundus images and 3D Optic Coherent Tomography (OCT) images modality as a clinical indicator for DR diagnosis because of its low cost and the simplicity in computation [11,12]. Fundus image modality is more often used than OCT images in both screening and clinical diagnosis systems. The lesions structure of the retina such as micro-aneurysms (MA), neovascularisation, haemorrhages, hard exudates, and cotton-wool-spots can be more prominent in fundus retinal images. In retinal fundus image diagnosis, there are two objects to be segmented: the anatomical structure (optic disc, blood vessels, macula, and fovea) and the lesions caused by diseases. The identification of the anatomical structure (main region) of the retinal image can assist in the analysis of conditions that affect these regions. For example, glaucoma affects the optic disc by making it more dilated, hypertension affects the vessels by tortuosity, and proliferative DR affects the number of new vessels [13]. The latter is the focal point of the present work. The anatomical attributes of retinal blood vessels, such as width, length, tortuosity, local branching pattern, and angles, play an important role in diagnostic results. They are widely used not only to predict ophthalmological disorders (i.e., DR, glaucoma, and age-related macular degeneration), but also to non-ophthalmological disorders (i.e., hypertension, arteriosclerosis, obesity, and cardiovascular disease) [14,15]. Therefore, the analysis of the vasculature structure in retinal fundus images allows both monitoring and evaluating the progress of these disorders.

Manual segmentation of the vascular tree in retinal fundus images is a tedious task for many reasons [16]. First, the variation in the contrast between blood vessel pixels and background in fundus images can be very tiny. Second, the noise presence and the difference in illumination conditions can cause problems. Third, the blood vessel variations, such as width, shape, branching angles, and brightness, may not be readily observable. Finally, the presence of lesions and haemorrhage in the background of retinal images can introduce confusion. Hence, automatic segmentation is a necessity.

A substantial amount of work is reported in the literature for detecting blood vessels in colour retinal fundus images [9,10,17]. Since retinal image analysis is complicated, it usually fails when a stand-alone image processing technique is used [18]. According to [19], the most classical segmentation techniques on the fundus images encompass six domains: (i) pattern recognition, (ii) mathematical morphology, (iii) matched filter, (iv) vessel tracking, (v) model-based, (vi) parallel-hardware. Pattern recognition is considered the most accurate and robust [18]. Parallel-hardware methods can accelerate the speed of execution, but the accuracy may be relatively lowered [20]. Hybrid methods, which combine one or more different domains to improve system performance, are in common use. Mathematical morphology, matched filter, vessel tracking and model-based techniques can contribute to each other or be applied to pattern recognition techniques and/or parallel-hardware methods to build an efficient, automated, comprehensive blood vessel segmentation system.

Machine learning, supervised and unsupervised, can be used in retinal blood vessel segmentation. In supervised learning, where for each input data, there is a corresponding output label, one can use a feature vector and a gold standard (GS) dataset. The latter includes manual markings for blood vessels, to train a classifier to divide retinal colour image pixels into vessel pixels and non-vessel pixels. On the other hand, where labels cannot exist, unsupervised methods explore interesting patterns of input data. In unsupervised methods, no prior knowledge is required upon the pattern before segmenting; moreover, the run is faster than the supervised methods. Although the supervised approaches are more accurate than the unsupervised method, it requires a craft, features engineering and lots of extracted features to obtain good results in the segmentation of the retinal blood vessel [21,22,23,24].

In the last five years, a new dawn of an era has shone based on deep learning as an automatic features extraction such as Deep convolutional neural networks(CNN). CNN shows the highest robustness and effectiveness in segmenting the blood vessels. Unlike any conventional supervised approach, the deep learning neural network (DNN) methods are intelligent in automatically learning the features by using a predefined filter, which is used to extract features from the input image. Additionally to this, Transfer learning, an advanced approach of Deep learning, is recently used. Rather than building and training the network to learn features from scratch, Transfer learning concepts are uses as a pre-trained model on vast numbers of different images. So the model is already pre-learned by the tremendous number of different features and only need to fine-tune for the more suitable and relevant application. Since the proposed work belongs to the former category, previous work using supervised learning is reviewed next.

A supervised learning method as an ensemble system based on feature-based AdaBoost classifier (FABC) using a collection of a 41-D feature vector is constructed for segmentation, but a high relative number of features is necessary to get the accurate results [25]. A 7-D feature vector comprises a grey level and moment invariant based feature segments. An artificial neural network (ANN) is used as a classifier. The method approaches are unable to achieve the desired results without sophisticated preprocessing steps [26]. An ensemble system of bagged and boosting analysis composed of the 9-D feature vector is constructed from the gradient vector field, morphological transformation, line strength measures, and Gabor filter to handle healthy and pathological colour retinal images. The methods use a number of features, and a decision tree applies where the effects of weak learners are combined using bootstrap aggregation to increase the accuracy of the segmentation [27]. The retinal blood vessel morphology is exploited to identify the successive stage of the severity of the DR disease. A technique based on an ANN is used to automate the screening of DR. An ANN that contains one input layer, five hidden layers, and one output layer is used. The input to the ANN is derived from the Gabor filter and moment invariant based features. The input image is processed using a mathematical morphology operator, mean, and Gaussian filter. This technique is used to diagnose diseases like hypertension and diabetes mellitus, for which a change in the retinal vasculature is a primary symptom. The method is applied only on the one dataset while it fails on another dataset [17]. The retinal vessel is uprooted using Lattice Neural Networks, and dendritic processing introduced relatively better performance in vessels. Nevertheless, the final results are low in DRIVE, STARE datasets [28]. Extensive preprocessing steps, such as zero-phase whitening, global contrast normalisation and gamma correlation on fundus image patches, are used as inputs to CNN. Although the prediction on the DRIVE, STARE, and CHASE_DB1, showed reduced vessels misclassification and central vessels reflex problem, the results are mainly based on preprocessing steps [29]. Five stages of DNN with the autoencoder transform the selected RGB patches of the retinal image into a vessel map image. Even though the prediction is made at every pixel in a patch, final results are based on the boundless cross-modality training on datasets [24]. A supervised hierarchical retinal blood vessel segmentation method is presented to diagnose several ailments, such as cardiovascular diseases, DR, and hypertension. The proposed method uses two superior layers: A convolution neural network, as a deep hierarchical feature extractor, and a random forest, which works as a trainable classifier. This method is tested only on two datasets. Although it outperforms the aforementioned methods, the time cost of this method is relatively high [30]. A supervised method based on simple extreme machine learning uses a 39-D feature vector constructed from local features (29 features), morphological transformation operators (6 features), phase congruency (1 feature), Hessian components (2 features) and divergence of a vector field (1 feature). Regardless of the fact that the method tries only one public dataset and another unknown private dataset, the obtained results are relatively low concerning the other publicly available datasets [12]. Information obtained from a stationary wavelet transform with a multi-scale convolution neural network can catch the tiny vessel pixel. This method uses a rotation operator as a basis of a joint strategy for data augmentation and prediction [31]. A hybrid feature set and hierarchical classification used to segment retinal blood vessels from a retinal colour image. Six groups of discerning features were used to construct the feature vector compose of local intensity variation features, the local binary pattern features, a histogram of gradient features, vector field divergence features, morphological transformation, and high-order-local auto-correlation features [32]. Moreover, a random forest-based technique was used as a feature selection in addition to using two hierarchical trainers, comprised of a support vector machine (SVM) and random forest, used as weak classifiers. At the same time, the final results are from a Bayesian network. The final output results completely depend on the random forest technique, as a feature selection in addition to hierarchical classification. However, it was so low when it tested for each classifier separately. A hybrid segmentation method with a self-mask generated schema was used to extract the retinal blood vessel from the fundus retinal image. However, its performance was very low on the DRIVE dataset [4]. A transfer learning model, VGG16, as a pre-trained model with multilevel/multi-scale deep supervision layers was incorporated to segment the retinal blood vessel. The VGG16 model was improved by convolving a vessel-specific Gaussian kernel with two pre-initialised scales. The final output is an activation map with well-learned vessel-specific feature at multiple scales. The activation map is increasingly obtained by the symmetric feature map. Regardless of using a pre-trained model, VGG16, which trained on ImageNet (15 million of high-resolution labelled images), the method required better post-processing steps [6]. Despite the increasing use of deep learning, there are some drawbacks that face it. First, the need for copious numbers of images to train and test, which does not exist in our case, and second, a high amount of power in processing using a graphical processing unit (GPU). Furthermore, it has a complicated, long structure, which needs a lot of experience, training and effort.

For fulfilling the retinal blood vessel segmentation, this paper uses hand-crafted features and a multi-layer perceptron neural network (MLP), with the following contribution.

• Carefully, 24 hand-crafted features are outlined.
• Special multi-layer perceptron neural network structure for full gauge fundus retinal evaluation.
• Cross entropy loss function.
• Transformed green channel image.
• Morphological transformation are applied twice. First, white top-hat and black bottom-hat are used as filters to extract features and to suppress noise during the extraction. Second, a closing operator is used to connect the isolated vessels and to fill the gap areas in the segmented images during the post-processing step.
• Model for automatic feature extraction and classification of blood vessel in retinal fundus images.
• Maximum and minimum moments of phase congruency are used for the first time.
• Retinal vessel segmentation with improved generalisation on multiple datasets.
• In-depth visual analysis plus comparison with the state-of-the-art methods.

This paper is organised as follows. Section 2 outlines the datasets and the methodology, where feature extracting steps are discussed in detail. In Section 3, the results of the proposed method are presented. The discussion of the results is in Section 4. The conclusions of the paper are described in Section 5.

2. Materials and Methods

Three datasets are commonly in use for vessel segmentation research. These are the Digital Retinal Images for Vessel Extraction (DRIVE) dataset [33], the Structure Analysis of Retina (STARE) dataset [34] and CHASE_DB1 dataset [35]. Because we do some early processing on STARE and CHASE_DB datasets in creating mask images related to each one, the difference between the three datasets must be explained as follows.

1.

DRIVE dataset

The DRIVE dataset consists of 40 images (7 of which are pathology) divided into 20 images for training and 20 for testing (3 of which are pathology). The image size is 565 × 584 pixels captured by a Canon CR5 non-mydriatic 3 charge-coupled-device (CCD) cameras at a 45° field of view (FOV). The image has a FOV approximately 540 pixels in diameter. The DRIVE dataset provides a mask for every image to facilitate the identification of the field of view (FOV) area, as shown in Figure 1.

Every image is accompanied by a manual segmentation label, manually created by three experts (observers) and validated by an experienced ophthalmologist (annotator). The training set is segmented only once and twice for the testing set, resulting in 577,649 pixels ($12.7\%$) marked as a vessel in the first set (first observer), against 556,532 pixels ($12.3\%$) in the second set (second observer). In this paper, the first set (first observer) is used as the ground truth (label) for the performance evaluation, while the second set is serving as a human observer reference for the performance comparison. The DRIVE dataset is depicted in Figure 2 and Figure 3.

2.

STARE dataset

The STARE (Structure Analysis of Retina) dataset comprises 20 images (10 of which present pathology) of size $700\times 605$ pixels. The dataset is captured by a TopCon TRV-50 fundus camera at 35° FOV. The dataset is manually segmented twice by two observers; the first one segments $10.4\%$ of pixels as vessels and $14.9\%$ for the other. Since the blood vessels of the second observer are much thinner than the first, the first is considered the ground truth. Figure 4 illustrates a sample of the STARE dataset. Unlike DRIVE, STARE does not offer masks to determine the FOV, so the method described in [11,34] is leveraged to create them.

3.

The CHASE_DB1 Dataset

The CHASE_DB1 Dataset includes 28 images of size $990\times 960$ pixels, captured from 14 school children for both left and right eyes by a hand-held NidekNM-200-D fundus camera at 30° FOV. The segmentation results of the first observer are deployed as the ground truth. This dataset is characterised by nonuniform background illumination. Thus, the contrast of blood vessels is inferior, with the presence of central vessel reflexes. The CHASE_DB1 is pictured in Figure 5. Like STARE, CHASE_DB1 does not offer mask images, so as in the above, the method described in [11,34] is also leveraged to create them.

The proposed model to segment retinal blood vessels from a colour retinal image contains four stages: (i) green channel extraction, (ii) features extraction, (iii) evaluation of segmentation, (iv) post-processing. In the first stage, the green channel is extracted from the RGB retinal fundus image. The choice of the green channel, in particular, is due to the fact that in this channel, the blood vessel tree is more vivid than in other channels. In the second stage, a 24 feature vector is constructed for pixel representation. The features of this vector encode information on the local intensity value (1 feature), morphological white top-hat value at six different scales (6 features), morphology black-bottom-hat values at six different scales (6 features), maximum-moment of phase congruency (1 feature), minimum-moment of phase congruency (1 feature), Hessian components values (4 features) and difference of Gaussian features at five different scales (5 features). At the end of this stage, a features matrix based on the feature vector and the correspondence manual label is constructed for each pixel in the green channel image. Later, a training sample from this feature matrix is randomly selected and used as input to the classifier to perform the segmentation. In the third stage, a supervised learning approach based on a multi-layer perceptron neural network is adopted. The rule for blood vessel extraction is to learn by the algorithm through the provided training sample. The segmented reference images are usually termed label or gold standard images. The last stage, post-processing, is used to enhance the classifier performance and acts to reduce two types of artefacts. The proposed method is tested on the three publicly available datasets, DRIVE, STARE, and CHASE_DB1, which contain digital colour retinal images with 8 bits per colour channel. A complete block diagram for the proposed model is shown in Figure 6.

The remainder of this section is dedicated to explaining the method in detail.

2.1. Feature Extraction

As mentioned above, in the second stage of our proposed method, a 24 feature vector is constructed for pixel representation. The features composing the vector are as follows.

2.1.1. Local Intensity Feature

In a retinal fundus image, blood vessel pixels appear much darker than the background pixels. So, the intensity value for each pixel in the retinal green channel is considered as a local intensity feature. Figure 7 depicts an image after a local intensity feature is extracted from the colour image.

2.1.2. Morphological White Top-Hat Transformation

In general, mathematical morphology transformation in image processing can combine (merge) two sets of inputs, the first related to the image and the second a predefined linear structuring element. One type of such transformation, white top-hat, is an operator used for lighting objects on a dark background. In the beginning, the opening operator removes objects smaller than the structuring element. Then the difference between the original image and the output of the opening operator is eliminated. Specifically, the opening operator uses the predefined structuring element, which is oriented at a particular angle $\theta$, to eradicate a vessel or part of it when the vessel’s width is shorter than the element and, both are orthogonal [15].

Let I be the green image to be processed, and S be the linear structuring element. Then the white top-hat transformation feature of I is given by

$$\begin{array}{ccc}\hfill T_c^{\theta }& =& I-I\circ S_c^{\theta },\hfill \end{array}$$

(1)

where $\circ$ is the opening operator, $c\in \left\{3,7,11,15,19,23\right\}$ the length of the structuring element and $\theta$ the rotation angle of S, spanning $\left[0:\pi \right]$ in steps of $\pi /15$. The result of the white top-hat transformation in all direction is given by

$$\begin{array}{ccc}\hfill T_c& =& \sum _{\theta \in A}{T}_c^{\theta },\hfill \end{array}$$

(2)

where $A=\left\{x\mid 0\le x\le \pi andxmod(\pi /12)=0\right\}$, and set A is a set that contains the orientation angles.

The length of S is critical to extract the vessel with the largest diameter. Therefore, a set of different lengths is used to formulate six different features. In case the opening operator $\circ$ along a class of S is considered, the summation of white top-hat transformation at each direction will brighten the vessels regardless of its direction. The effectiveness of the summation on the entire retinal image will enhance all vessels regardless of their orientation. That includes the small (tortuous) vessels in the bright zone. Figure 8 shows the image after applying the white top-hat operator.

Following the white top-hat feature extractor computation proposed in [15], six different features are extracted. Not only do vessel-like structures become much lighter but they also acquire intensity values relative to the local background. So, the intensity values of extracted blood vessels are relative to the local intensity of neighbourhood pixels in the original image [36]. The final response of the white top-hat transformation operator is a corrected uniform image with enhanced edge information.

2.1.3. Morphological Black Bottom-Hat Transformation

The counterpart of morphological white top-hat is the morphological black bottom-hat transformation. While the former is used for light objects on a dark background, the latter is used for the opposite. Black bottom-hat transformation is defined as the residue of the closing operator and the original image. In other words, it is a closing • of the image minus the image [12]. The black bottom-hat transformation feature can be defined as follows

Let I be the green image to be processed, S is the structuring element, • is the closing operator and $\theta$ is an orientation angle at a given scale c. Then the black bottom-hat transformation feature on image I is given by

$$\begin{array}{ccc}\hfill B_c^{\theta }& =& I•S_c^{\theta }-I,\hfill \end{array}$$

(3)

Now,

$$\begin{array}{ccc}\hfill B_c& =& \sum _{\theta \in A}{B}_c^{\theta }.\hfill \end{array}$$

(4)

where $B_c$ is the summation of the black bottom-hat transformation in all directions of S at a different $\theta$. Both c and $\theta$ are defined as in the white top-hat transformation.

The image is simplified by both of the above morphological filters. So the cross-curvature computation is easy, and the vessel segments are linearly coherent. Figure 9 shows the image after applying black top-hat operators.

2.1.4. The Principal Moments of Phase Congruency

Instead of extracting features by using intensity gradient-based methods in the spatial domain, the phase congruency method is based on Fourier components. Its importance stems from its ability to detect edges and corner points in an image effectively and invariantly to change in contrast and illumination conditions. Furthermore, phase congruency has a significant advantage over all other feature detectors as it can mark features that discriminate like edges (blood vessels) correctly. It is based on the local energy model, reported by Morrone et al. [37], which suggests that features as line and edges are perceived at points where the Fourier components are maximal in phase. The model is defined as a ratio of local energy to the overall path length used up by Fourier components from the head to tail.

Let x be some location at a signal, the nth Fourier component of which at point x has amplitude $A_n$, and let $\left|E\right(x\left)\right|$ be the magnitude at that point of the vector from the head to tail (considered as the local energy). Then the phase congruency, ($PC)$, of the signal is given by

$$PC=\frac{\left|E\right(x\left)\right|}{{\sum }_n{A}_n\left(x\right)}.$$

(5)

In case all the Fourier components are in phase, then all the complex vectors will be aligned, and the $PC$ will be 1; otherwise, it will be 0. While the model is not based on gradient intensity for detecting features, the construction of a dimensionless measure $PC$ at any point in the image is possible. The most notable drawback of this model is the bad localisation and its sensitivity to noise. A modified model is proposed by Kovesi [38,39], to compensate for the image noise. It provides proper localisation since $PC$ is represented as the phase difference between the amplitudes of the weighted mean of the local phase energy angle $\overline{\varphi }\left(x\right)$ and the local phase of nth Fourier component $\varphi \left(x\right)$, both evaluated at point x. As such, the modified model can be expressed as follows.

$$PC(x,y)=\frac{{\sum }_i^I{\sum }_k^K{w}_i(x,y)⌊A_{i,k}(x,y)\left(cos{(\varphi }_{i,k}(x,y)-\overline{\varphi }(x,y))-\left|sin{(\varphi }_{i,k}{(x,y)-\overline{\varphi }}_{i,k}(x,y))\right|\right)-T⌋}{{\sum }_i^I{\sum }_k^K{A}_{i,k}(x,y)+\eta },$$

(6)

where ($x,y$) is the (green channel) pixel coordinate, I is the total number of orientation. K is the total number of scales, $w_i(x,y)$ is a weighted factor for frequency spread in orientation i (because $PC$ from many frequencies is significantly better than $PC$ over a few frequencies), T is the estimated noise influence (recommended value is 3), and $\eta$ is a constant incorporated to avoid division by zero if the local energy vanishes. In Equation (6), a cosine of the phase difference between the amplitude of $\overline{\varphi }\left(x\right)$ and $\varphi \left(x\right)$ at a particular point is used instead of a sine because when the frequency components are aligned (in phase) with the “local energy”, the phase difference value can be expressed as ${\varphi }_n-{\overline{\varphi }}_n=0$, with the cosine being 1, resulting in a maximum $PC$. The weighting function w used to compute frequencies in ranges is defined as a sigmoid function

$$w(x,y)=\frac{1}{1+{exp}^{\left(\rho \right(c-s(x,y))}},$$

(7)

where c is the filter response cut-off value, below which $PC$ values are ignored, $\rho$ is the value of the gain factor that controls the sharpness of the cut-off (its value set to 10), and the fractional value of $s(x,y)$ of the frequencies is calculated by taking the amplitude of the sum of the filter response and dividing by the maximum amplitude $A_{max}(x,y)$ at each point (pixel) (x,y) of the images as follows.

$$s(x,y)=\frac{1}{K}\left(\frac{{\displaystyle \sum _k^K}{A}_k(x,y)}{A_{max}(x,y)+\eta }\right).$$

(8)

Herein, Equation (6) represents $PC$ for all orientation from $[0:\pi ]$ overall scales K in a 2D image. Unlike the proposed method in [12], the modified model in [39,40] is used by calculating $PC$ independently for each orientation using Equation (6), then the moment of $PC$ covariance is computed. Hence, the variation of the moment within an orientation is taken into account. In such a case, the principal axis is corresponding to the axis about which the magnitude of the moment is minimum (minimum moment), which provides a good indicator for the orientation of the feature. On the other side, the magnitude of the maximum moment corresponding to the moment about the axis perpendicular to the principal axis is considered as an indicator of the significance of the feature. Accordingly, the maximum moment of $PC$ is used directly to establish whether the point in the retinal image has a significant edge feature (blood vessel), where the magnitude has a large value. Meanwhile, if the minimum moment of $PC$ is also large then it is considered as a strong indicator that the feature point has a strong 2D component and can be considered a corner.

The calculation of the maximum and minimum moment of $PC$ is represented in [39,40] based on a classical moment analysis equation. Firstly, the $PC$ information is combined in an orientation matrix ($Om$) given by

$$Om=\left[\begin{array}{cc}{(P{C}_{o}cos\left(\theta \right))}^2& (P{C}_{o}cos\left(\theta \right)(P{C}_{o}sin\left(\theta \right))\\ (P{C}_{o}sin\left(\theta \right))(P{C}_{o}cos\left(\theta \right))& {(P{C}_{o}sin\left(\theta \right))}^2\end{array}\right],$$

(9)

evaluated over all multiple orientations formulated from the $PC$ covariance matrix ($Cm$) given by

$$Cm=\left[\begin{array}{cc}P{C}_x^2& P{C}_{x}P{C}_y\\ P{C}_{y}P{C}_x& P{C}_y^2\end{array}\right],$$

(10)

where $P{C}_o$ is the amplitude of $PC$, which is determined at orientation o, and the sum is performed overall six orientations. Figure 10 vividly depicts the phase congruency filter response at the six different directions. According to the above illustration of the PC model, in case all Fourier components are in phase, all the complex vectors would be aligned, indicating that the feature point has a robust 2D component and can be considered an edge like feature.

A singular value decomposition on a $PC$ covariance matrix is performed, and the maximum moment, (M), and minimum moment, (m), are obtained as follows.

$$\begin{array}{ccc}\hfill M& =& \frac{1}{2}\left(\gamma +\alpha +\sqrt{{\beta }^2+{(\alpha -\gamma )}^2}\right),\hfill \\ \hfill m& =& \frac{1}{2}\left(\gamma +\alpha -\sqrt{{\beta }^2+{(\alpha -\gamma )}^2}\right),\hfill \end{array}$$

where

$$\begin{array}{ccc}\hfill \alpha & =& \sum _o{(P{C}_{o}cos\left(\theta \right))}^2,\hfill \\ \hfill \beta & =& 2\sum _0(P{C}_{o}cos\left(\theta \right))(P{C}_{o}sin\left(\theta \right)),\hfill \\ \hfill \gamma & =& \sum _o{(P{C}_{o}sin\left(\theta \right))}^2.\hfill \end{array}$$

(11)

Figure 11 visualises the extracted features image of the two different moments (M) and (m) features.

2.1.5. Multi-Scale Second-Order of Local Image Structure (Hessian matrix)

The Hessian matrix $\mathbf{H}(x,y)$ for a point $I(x_i,y_j)$ in the 2D retinal image is used for vessel detection and is given by

$$\mathbf{H}(x,y)=\left[\begin{array}{cc}\frac{{\partial }^{2}f}{\partial x^2}& \frac{{\partial }^{2}f}{\partial x\partial y}\\ \frac{{\partial }^{2}f}{\partial y\partial x}& \frac{{\partial }^{2}f}{\partial y^2}\end{array}\right]=\left[\begin{array}{cc}{f}_{xx}& f_{xy}\\ f_{yx}& f_{yy}\end{array}\right],$$

(12)

where the second partial derivatives components are presented as $f_{xx}$, $f_{xy}$, $f_{yx},$ and $f_{yy}$ for the neighbours of the point $(x,y)$. The Hessian matrix is exploited to describe the second-order structure of intensity variation value around a pixel on the retinal image. The information obtained from the $\mathbf{H}$ describes the local curvature of the data in a small neighbourhood surrounding each pixel. The often employed approach to the retinal image is inspired by the theory of linear scale-space [36]. Accordingly, a multi-scale algorithm is designed to identify vascular structure by investigating the signs and magnitudes of the eigenvalues of $\mathbf{H}$, which reflect specific shapes of structures in the retinal image, besides their brightness.

Let $I(x,y)$ be the intensity of the a $2D$ green retinal image at a particular point $(x,y)$ in an image, and $G_{2D}(x,y,\sigma )$ is a $2D$ Gaussian kernel. Then, $2\times 2$ Hessian matrix $\tilde{\mathbf{H}}(x,y,\sigma )$ components of the second order derivatives of $G_{2D}(x,y,\sigma )$ multiplied by $I(x,y)$ at a scale $\sigma$ are defined as:

$$\begin{array}{ccc}\hfill f_{xx}& =& I(x,y)\ast \frac{{\partial }^2}{\partial x}{G}_{2D}(x,y,\sigma ),\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill f_{xy}& =& f_{yx}=I(x,y)\ast \frac{{\partial }^2}{\partial x\partial y}{G}_{2D}(x,y,\sigma ),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill f_{yy}& =& I(x,y)\ast \frac{{\partial }^2}{\partial y}{G}_{2D}(x,y,\sigma ),\hfill \end{array}$$

(15)

where the $2D$ Gaussian kernel is given by

$$G_{2D}(x,y,\sigma )=\frac{1}{2\pi {\sigma }^2}exp\left(\frac{-x^2+y^2}{2{\sigma }^2}\right).$$

(16)

with $\sigma$ the standard deviation of the Gaussian distribution, referring to a scale of the filter, and ∗ the convolution operator. Different values of $\sigma$ are experimentally tested, and the optimal values of $\sigma$ are detected as $\frac{1}{2}$, 1, 2, and $2\sqrt{2}$. The response of second-order derivative of a Gaussian function at a scale $\sigma$ represents a probe kernel that can be a measure of the contrast between regions inside and outside the given range ($-\sigma ,\sigma$) in the derivative direction in the retinal image [41,42]. In this respect, the blood vessel pixels gain the highest response at a suitable scale of $\sigma$, while the background pixels get a low response. It is noteworthy that the nature of the second-order derivative of the Hessian matrix is invariant to linear grey variation of the retinal image pixels, because of the smooth nature of the Gaussian kernel. So only local contrast is assessed along with the actual greyscale values.

As soon as the elements of the Hessian matrix are obtained, the corresponding eigenvalues and eigenvectors are determined. The eigenvalues, ${\lambda }_{{}_{{}_i}}$, have an essential role in the discrimination of the local orientation pattern. Eigenvalues not only represent the magnitude of maximal local contrast change but also the change of the magnitude of the local intensity in the other perpendicular direction. Consequently, the eigenvalue decomposition of Hessian can be exploited to distinguish between the edge-like feature in a fundus retinal image and the background. Concerning blood vessel skeleton, darker features are considered. We are interested in only the case where one eigenvalue has a positive and high magnitude.

Table 1. The relations between the eigenvalues for different features in any 2D image are adopted from Frangi’ Schema where H = High-value; L = low-value; N = Noise usually small; +/ − = positive or negative eigenvalue; Dark = dark feature on a bright background (blood vessel in our case [41]. The representation structure we are interested in is indicated with a *.

2D Image
${\mathbf{\lambda }}_{\mathbf{1}}$	${\mathbf{\lambda }}_{\mathbf{2}}$	The Pattern
N	N	Noisy, no preferred direction
L	H+	Tubular structure (dark)
L	H-	Tubular structure (bright)
H+	H+	Blob-like structure (dark) *
H-	H-	Blob-like structure (bright)

illustrates all possible variants of the structure’s shape and brightness that can be identified through the analysis of Hessian eigenvalues in any 2D medical image. We are only interested in the fundus retinal image modality (in our case, the blob-like structures are dark).

Geometric information for the blob-like structure is obtained through a dissimilarity measure based on the second-ordered derivative that is called blobness ratio, $R_b$, given by

$$R_b=\frac{{\lambda }_1}{{\lambda }_2}.$$

(17)

This ratio is maximal for a blob-like structure and invariant to the grey level. Since the intensity value of the blood vessel pixels is lower (darker) than the background pixels, the magnitude of the second-derivative (eigenvalues) will be high at the blob-like structure and small at the background. So the ratio will be higher at the blood vessel pixels and lower at the background pixels.

The magnitude of the second-order derivatives at the background is small, distinguishing background structure from the blood vessels. The norm of Hessian is used to quantify this property in terms of the structureness ratio, S, and is given by

$$S={\left|\right|\mathbf{H}\left|\right|}_F=\sqrt{(f_{xx}^2+2f_{xy}^2+f_{yy}^2)},$$

(18)

where ${\left|\right|.\left|\right|}_F$ is the Forbenius norm of the components of Hessian matrix. Intuitively, the structureness ratio S is small for the background pixels mostly because of their lack of structure and higher amount of blood vessel pixels.

Finally, a vessel function $V\left(S\right)$ combines the components of the blobness ratio $R_b$, and structureness ratio, $V\left(S\right)$, to map these features into a probability-like estimates $V_o\left(S\right)$ of being a vessel given by

$$V_o\left(S\right)=\left\{\begin{array}{cc}0,\hfill & {\lambda }_2>0\hfill \\ exp(-\frac{R_b^2}{2{\beta }^2})(1-exp(-\frac{S^2}{2c^2})),\hfill & \mathrm{elsewhere}\end{array}\right.,$$

(19)

where $\beta$ is equal to $0.5$ and c is half of the maximum Frobenius norm S. Both V and S are computed at a different scale $\sigma$. Herein, $R_b$ accounts for the eccentricity of the second-order ellipse. Then the maximum value for both V, S are taken as a Hessian feature for each pixel and they are given by:

$$\begin{array}{ccc}\hfill V_{Max}& =& max(V_{0.5},V_1,V_2,V_{2\sqrt{2}}),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill S_{Max}& =& max(S_{0.5},S_1,S_2,S_{2\sqrt{2}}).\hfill \end{array}$$

(21)

Both $V_{Max}$, $S_{Max}$, ${\lambda }_1$ and ${\lambda }_2$ are used as a features for every pixel of a green image.

The visualisation of the four features of the Hessian matrix components are delineated in Figure 12.

2.1.6. Difference of Gaussian

Difference of Gaussian (DoG) kernels are exploited as an edge detector [16]. A DoG kernel is much closer to the second-order of Gaussian function, which is described in detail in Section 2.1.5. To characterise a pixels features regarding its neighbours, a DoG filter is utilised at five different smoothing scales, $\sigma \in \{\frac{\sqrt{2}}{2},1,\sqrt{2},2,2\sqrt{2}\}$, with the base scale set to $\sigma =0.5$. Using DoG as a band path filter not only serves to remove noise in an image but also vividly depicts the homogeneous area in an image. The visualisation of the five DoG features are depicted in Figure 13.

In this paper, a 24-D feature vector is constructed for each pixel in a retinal image and is used to differentiate vessel pixels, lesions, and background. Five groups of features are used. The RGB image is transformed into a green channel image because the contrast between blood vessel pixels and background is more precise and prominent than the other channels. So, the intensity value of each pixel in the green channel image is used as a feature. The moments of phase congruency emphasise the edges of the retinal blood vessel and it is the first time to be used as a feature in the fundus retinal image. The Hessian matrix components emphasise the zero crossing. The white top-hat transformation using multi-scales and multi-orientation structuring elements are efficient in making a corrected uniform image. Black bottom-hat transformation using multi-scales and a multi-orientation structuring element is used as a filter to suppress noise. The difference of Gaussian is used, and the response of the filter is a homogeneous image.

2.2. Segmentation

The procedure can be distinguished in three steps: first, scaling the features by normalisation; second, design of the neural network, in which the configuration is decided and trained (applied what the NN learned); third, testing, to identify each pixel as a vessel or a non-vessel to get the final binary image.

2.2.1. Normalisation

In the end, in the feature extraction process, each pixel in the green channel image is characterised by a vector $V(x,y)$ (in a 24-D feature space) given by

$$V(x,y)=(v_1(x,y),v_2(x,y),v_3(x,y),\dots ,v_{24}(x,y)),$$

(22)

where the features $v_j$ have distinct ranges and values. Normalisation, by narrowing down the range between the different features is necessary and leads to an easier classification procedure of assigning the candidate pixel to the suitable class: class $K_{{}_1}$ of blood vessel pixel or class $K_2$ of non-vessel pixel. One can obtain a normalised version ${\overline{v}}_j$ of features $v_j$ by

$${\overline{v}}_j=\frac{v_j-{\mu }_j}{{\sigma }_j},$$

(23)

where ${\mu }_j$ is the mean value of the jth feature and ${\sigma }_j$ their standard deviation. Finally, the features are normalised with a zero mean and unit variance for each pixel independently. The resulting features are used in the classification procedure.

According to [26], using the results of linear classifiers to separate the two classes of vessel segmentation gives poor results because of the nature of the distribution of the training dataset at the feature space, which is nonlinear. Consequently, segmentation can be made more accurately and efficiently using a nonlinear classifier, such as support vector machine (SVM) [32], neural network [24], a multi-layer feed forward neural network [17], decision tree [43], random forest [32], extreme learning machine [12] and convolution neural network(CNN) [44].

2.2.2. Design of the Neural Network

In the present work, we use a modified MLP neural network of one input layer (with 24 input units), three hidden layers, and one output layer as a binary classifier. In the purpose of getting the optimal topology of hidden layers, several topologies and different numbers of neurons were tested. The test results showed that excellent results are obtained when each of the three hidden layers contains 27 units. A rectified linear unit (Relu) activation function is selected for the hidden layers [45]. Further, in classification, cross-entropy (CE) was found better to use as a loss function than a Minimum Mean Square Error (MSE) loss function, because of the unbalanced distribution of vessel and non-vessel pixels [45], where the CE measures the nearness of two probability distributions.

2.2.3. Training

To train the classifier, a set $S_{Train}$ of M candidates features for which the feature vector [$V-$(22)], and the classification result, $K_1$ (vessel pixel) or $K_2$ (non-vessel) are known. The sample that forms $S_{Train}$ is collected randomly from the manually labelled pixels, vessel and non-vessel, from the training images as proposed in [26,43] for each dataset. Once the classifier passes the training steps, it becomes ready to work and classify images of which the results are not known a priori.

2.2.4. Application of Neural Network (Testing)

The trained neural network is used to test “unseen” dataset images. Specifically, a binary image is generated in which blood vessel pixels are distinguished from background pixels. The mathematical descriptions of pixels are individually input to MLP. In particular, the input layer units receive a vector V of the normalised features, as in Equations (22) and (23). The Relu function has a tremendous out-of-sample prediction ability due to its simplicity. It also alleviates the problem of a vanishing gradient in the learning phase of the ANN [46]. Finally, we obtained a binary image that contains only extracted blood vessel and background pixels.

In this paper, we use a multi-layer neural network (MLP) as a binary classifier, since there are only two classes (blood vessel) and (non-vessel). MLP is an extension of an artificial neural network (ANN), considered to be more nonlinear when mapping a set of the input variable to a set of output variables [47]. Also, it does not make any assumption regarding the probabilistic information about the pattern under consideration.

2.3. Post-Processing

In the post-processing stage, two steps are used to reduce two types of classification artefacts. First, the discontinuous vessels pixels, which are surrounded by vessel neighbours and are misclassified as background pixels ($FN$). Second, the non-vessel pixels, which appear as small isolated regions and are misclassified as vessels ($FP$). For the sake of overcoming both artefacts, a mathematical operator is used. A mathematical closing operator is applied for the first artefact to fill any holes smaller than a predefined structure (a disk of radius 9 in our experiments) in the vessel map. For the other artefact, the acreage of isolated pixels are connected in eight connected areas, then any area below 100 pixels is removed, as proposed by [9]. Once both artefacts are corrected, they will be reclassified correctly as vessel pixels for the first artefact and background pixels for the second.

The experimental results in Section 3 show that the combination of the five groups of features, 24D features, which are encoded to the modified MLP, yields great performance. Our experiments are tested in the environment of 2.4 GHz Intel7-77ooQH CPU and 16 GB memory using Scikit-learn: open source machine learning libraries [48].

3. Results

The segmentation performance of the proposed method on the fundus image can be assessed by comparing the segmentation testing results with the GS image label. To this end, we coded our method and carried out numerous experiments on the three publicly available datasets DRIVE, STARE, and CHASE_DB, as mentioned at the beginning of this section.

3.1. Performance Measures

The classification result of every pixel can be one of the four types shown in the confusion matrix of

Table 2. The contingency table for vessel segmentation.

	Gold Standard (GS) Segmentation
	Vessel Pixels	Background Pixels
Vessel Pixels (Predicted)	True-Positive ($TP$)	False Positive ($FP$)
Background pixels (Predicted)	False Negative ($FN$)	True Negative ($TN$)

, where the entries are as follows. True positive ($TP$) indicates a vessel pixel classified correctly as a vessel pixel. False-positive ($FP$) indicates a non-vessel pixel classified wrongly as a vessel pixel. True negative ($TN$) indicates a non-vessel pixel classified correctly as a non-vessel pixel. False negative ($FN$) indicates a vessel pixel classified wrongly as a non-vessel pixel.

We started by training our classifier using a set of images selected randomly from three datasets: DRIVE, STARE and CHASE_DB1. From each dataset, a sample of around $30\%$ of the $S_{Train}$ are used as a training sample set. Firstly, the RGB image is transformed into a green channel image, since the green channel shows the best contrast between vessel and background [11,18,19,20,25]. Then, for every pixel in the image, a collection of hybrid features including local intensity value, morphology top-hat transformation, morphology bottom-hat transformation, minimum/maximum moment of phase congruency, Hessian matrix components and difference of Gaussian features are extracted from the green band image to generate a 24-D feature vector, which is represented as 24 different features for each pixel.

The design of the feature vector characterises every pixel in terms of some quantifiable measurements that can be easily used to differentiate effectively between a vessel and non-vessel pixel in the classification step, as we described in detail in Section 2. The extracted features are collected in a feature matrix, which represents all features from training images for each dataset. The ($S_{Train}$) are randomly selected from the feature matrix to teach the classifier how to differentiate between vessel pixels and background. After the MLP is trained, it is used to identify and differentiate between the two different patterns in the test images, blood vessel pixel and the background pixel.

To quantify the classification results, we use the following standard metrics: sensitivity ($Se$), specificity ($Sp$), positive prediction value ($Ppv$), negative prediction value ($Npv$), and accuracy ($Acc$), which are given by the equations below [49]. It should be noted, however, that the values of the variables in those equations are cumulative, i.e., denote the total number of pixels classified. For example, $TP$ in those equations denotes the total number of vessel pixels correctly classified as vessel pixels.

$$\begin{array}{cc}Se=\frac{TP}{TP+FN},& \hfill \end{array}$$

(24)

$$\begin{array}{cc}Sp=\frac{TN}{TN+FP},& \hfill \end{array}$$

(25)

$$\begin{array}{cc}Ppv=\frac{TP}{TP+FP},& \hfill \end{array}$$

(26)

$$\begin{array}{cc}Npv=\frac{TN}{TN+FN},& \hfill \end{array}$$

(27)

$$\begin{array}{cc}Acc=\frac{TP+TN}{TP+FP+TN+FN},& \hfill \end{array}$$

(28)

Additionally, the F1 metric is suggested for evaluating classification performance when the dataset is unbalanced, as is the case here, and is given by

$$\begin{array}{cc}F1=\frac{2TP}{FP+FN+2TP}.& \hfill \end{array}$$

(29)

All the above metrics are in [0,1], with 0 being the worst and 1 being the best (though generally not practically attainable.) In our tests, all these metrics are evaluated for each image, then the average for all images tested from each dataset is calculated. In our experiments, we tested images from all three datasets: DRIVE, STARE, and CHASE_DB1. For each image tested, the feature vectors were calculated for all pixels and applied consecutively to the classifier. Based on the classifier’s decision for each pixel, and on our knowledge of the truth, we calculated the total numbers of the $TP$, $FP$, $TN$, $FN$ cases for the image as a whole. These numbers are then plugged into the above equations to obtain the corresponding metrics.

3.2. The Proposed Method Evaluation

The last results are obtained from the probability map image (binary image) and the five performance metrics on the DRIVE, STRAE, CHASE_DB, respectively. The average of the selected measures is tabulated in

Table 3. The performance results on 20 images of the DRIVE dataset.

Image	$\mathit{Se}$	$\mathit{Sp}$	$\mathit{Ppv}$	$\mathit{Npv}$	$\mathit{Acc}$	$\mathit{F}1$
1	0.8198	0.9764	0.8354	0.9723	0.9611	0.8091
2	0.8171	0.9839	0.8975	0.9681	0.9651	0.8117
3	0.7592	0.9773	0.8319	0.9635	0.9535	0.7516
4	0.7940	0.9829	0.8699	0.9694	0.9636	0.7861
5	0.7916	0.9876	0.8943	0.9717	0.9605	0.6572
6	0.6639	0.9915	0.9232	0.9489	0.9548	0.8117
7	0.8189	0.9843	0.7852	0.9782	0.9598	0.7993
8	0.5888	0.9951	0.9391	0.9474	0.9605	0.5829
9	0.7434	0.9785	0.7877	0.9717	0.9565	0.7359
10	0.8063	0.9782	0.8099	0.9767	0.9630	0.7983
11	0.8181	0.9792	0.8425	0.9742	0.9592	0.8100
12	0.7675	0.9824	0.8504	0.9689	0.9623	0.7598
13	0.6566	0.9912	0.9275	0.9429	0.9567	0.6500
14	0.8143	0.9843	0.7852	0.9782	0.9598	0.8062
15	0.8282	0.9699	0.7696	0.9779	0.9625	0.8200
16	0.7588	0.9873	0.8931	0.9659	0.9621	.0.7512
17	0.6754	0.9894	0.8845	0.9611	0.9572	0.6686
18	0.7093	0.9919	0.9302	0.9560	0.9609	0.7022
19	0.7765	0.9862	0.9038	0.9622	0.9711	0.7687
20	0.6782	0.9895	0.9084	0.9511	0.9647	0.6714
Average	0.7542	0.9843	0.8634	0.9653	0.9607	0.7475

,

Table 4. Performance results on 20 images of the STARE dataset.

Image	$\mathit{Se}$	$\mathit{Sp}$	$\mathit{Ppv}$	$\mathit{Npv}$	$\mathit{Acc}$	$\mathit{F}1$
1	0.7904	0.9803	0.7939	0.9556	0.9484	0.7834
2	0.6933	0.9765	0.7446	0.9683	0.9562	0.6869
3	0.8612	0.9787	0.7629	0.9794	0.9683	0.8527
4	0.6868	0.9961	0.9029	0.9582	0.9628	0.6792
5	0.7536	0.9734	0.7772	0.9542	0.9421	0.7467
6	0.8253	0.9821	0.8076	0.9794	0.9716	0.8214
7	0.8369	0.9727	0.8009	0.9864	0.9705	0.8265
8	0.7654	0.9655	0.7472	0.9888	0.9659	0.7577
9	0.8752	0.9842	0.8646	0.9821	0.9752	0.8426
10	0.8512	0.9636	0.7350	0.9776	0.9545	0.8504
11	0.8762	0.9800	0.8168	0.9825	0.9728	0.8741
12	0.8877	0.9816	0.8423	0.9838	0.9757	0.8679
13	0.8382	0.9841	0.8754	0.9721	0.9680	0.8294
14	0.7546	0.9871	0.8974	0.9707	0.9694	0.7498
15	0.6491	0.9924	0.9150	0.9499	0.9538	0.6371
16	0.7546	0.9857	0.8843	0.9487	0.9484	0.7395
17	0.6591	0.9825	0.8665	0.9735	0.9679	0.6487
18	0.6218	0.9975	0.9416	0.9696	0.9726	0.6196
19	0.9678	0.9951	0.8726	0.9710	0.9742	0.9660
20	0.6644	0.9913	0.8346	0.9461	0.9473	0.6548
Average	0.7806	0.9825	0.8341	0.9698	0.9632	0.7717

and

Table 5. Performance results on 14 images of the CHASE_DB1 dataset.

Image	$\mathit{Se}$	$\mathit{Sp}$	$\mathit{Ppv}$	$\mathit{Npv}$	$\mathit{Acc}$	$\mathit{F}1$
1	0.7833	0.9871	0.7762	0.9658	0.9560	0.7754
2	0.8358	0.9778	0.7184	0.9710	0.9526	0.8274
3	0.6551	0.9889	0.8381	0.9673	0.9697	0.6485
4	0.6240	0.9812	0.8720	0.9652	0.9680	0.6177
5	0.8046	0.9702	0.6458	0.9703	0.9449	0.7965
6	0.8015	0.9585	0.6085	0.9786	0.9422	0.8924
7	0.8455	0.9927	0.7891	0.9785	0.9718	0.8370
8	0.7824	0.9963	0.8157	0.9785	0.9708	0.7745
9	0.7826	0.9829	0.7542	0.9741	0.9507	0.7747
10	0.8672	0.9809	0.7642	0.9638	0.9562	0.8585
11	0.6560	0.9951	0.9797	0.9714	09530	0.6495
12	06280	0.9934	0.7784	0.9595	0.9513	0.6220
13	0.7390	0.9958	0.8422	0.9559	0.9602	0.7318
14	0.8145	0.9849	0.7326	0.9626	0.9604	0.8063
Average	0.7585	0.9846	0.7796	0.9687	0.9577	0.7580

. We start with the DRIVE dataset, from which we randomly selected 20 images for testing. The results for the performance metrics are given in Table 3.

Table 4 and Table 5 show, respectively, the same metrics for the STARE and CHASE_DB1 datasets, where 20 images from the former and 14 from the latter were tested. At the bottom of each table, the averages of the metrics for that dataset are placed. Then segmentation results are identified from the obtained metrics on 20 images of the STARE dataset are presented in Table 4.

The best and worst segmentation results are identified from the obtained metrics and are presented in Table 5 on the CHASE_DB1 dataset.

For DRIVE dataset, the best case of segmentation in terms of Se, Sp, Ppv, Npv, Acc, and F1 are 0.8282, 0.9951, 0.9391, 0.9782, 0.9711, and 0.8200 respectively, while the worst case are 0.5888, 0.9699, 0.7696, 0.9429, 0.9535, and 0.5829 respectively. For the STARE dataset, the best case of segmentation in terms of Se, Sp, Ppv, Npv, Acc, and F1 are 0.9678, 0.9975, 0.9416, 0.9888, 0.9757, and 0.9660 respectively, while the worst case are 0.6218, 0.9636, 0.7350, 0.9461, 0.9421, and 0.6196 respectively. For the CHASEDB1 dataset, the best case of segmentation in terms of Se, Sp, Ppv, Npv, Acc, and F1 are 0.8672, 0.9963, 0.9797, 0.9718, and 0.8585 respectively, while the worst case are 0.6240, 0.9585, 0.6085, 0.9559, 0.9422, and 0.6177 respectively.

Table 6. Segmentation results for the best and the worst cases.

Dataset		$\mathit{Se}$	$\mathit{Sp}$	$\mathit{Ppv}$	$\mathit{Npv}$	$\mathit{Acc}$	$\mathit{F}1$
DRIVE	$Best(max)$	$0.8282$	$0.9951$	$0.9391$	$0.9782$	$0.9711$	$0.8200$
	$Worst(min)$	$0.5888$	$0.9699$	$0.7696$	$0.9429$	$0.9535$	$0.5829$
STARE	$Best(max)$	$0.9678$	$0.9975$	$0.9416$	$0.9888$	$0.9757$	$0.9660$
	$Worst(min)$	$0.6218$	$0.9636$	$0.7350$	$0.9461$	$0.9421$	$0.6196.$
CHASE_DB1	$Best(max)$	$0.8672$	$0.9963$	$0.9797$	$0.9786$	$0.9718$	$0.8585$
	$Worst(min)$	$0.6240$	$0.9585$	$0.6085$	$0.9559$	$0.9422$	$0.6177$

summarises the best and the worst results on the whole dataset under study using the proposed method.

3.3. Comparison with The-State-Of-Art Methods

This method is applied with hand-crafted features and a modified MLP with CE, to improve its performance compared with seven other state-of-the-art methods of the same category. The comparison is carried out on the three datasets mentioned above.

Table 7. Comparison with state-of-the-art-methods on the DRIVE dataset. N.A; not available by their authors.

Method	Year	Se	Sp	Acc	F1
Marin et al. [26]	2010	0.7067	0.9801	0.9452	N.A
Fraz et al. [27]	2014	0.7406	0.9807	0.9480	N.A
Zhao et al. [50]	2014	0.7187	0.9789	0.9509	N.A
Fu et al. [28]	2016	0.7444	0.9600	0.9412	0.6884
Somro et al. [51]	2017	0.7382	0.9123	0.9421	N.A
Khowaja et al. [32]	2019	0.7437	0.9596	0.9410	0.6877
Sundram et al. [4]	2019	0.6909	0.9401	0.9301	N.A
Proposed Method		0.7542	0.9843	0.9607	0.7475

,

Table 8. Comparison with state-of-the-art-methods on the STARE dataset, N.A; not available by their authors.

Methods	Year	Se	Sp	Acc	F1
Marin et al. [26]	2010	0.6944	0.9820	0.9526	N.A
Fraz et al. [43]	2014	0.7548	0.9763	0.9534	N.A
Zhoa et al. [50]	2014	0.7187	0.9767	0.9509	N.A
Fu et al. [52]	2016	0.7402	0.9479	0.9336	0.6187
Soomra et al. [51]	2019	0.7382	0.9123	0.9421	N.A
Khowaja et al. [32]	2019	0.7155	0.9746	0.9483	0.7224
Proposed Method		0.7806	0.9825	0.9632	0.7717

and

Table 9. Comparison with state-of-the- art-methods on the CHASE_DB1 dataset, N.A; not available by their authors.

Methods	Year	Se	Sp	Acc	F1
Fraz et al. [27]	2012	0.7224	0.9711	0.9469	N.A
Fraz et al. [43]	2014	0.7259	0.9770	0.9524	N.A
Azzopardi et al. [16]	2015	0.7585	0.9587	0.9587	N.A
Fu et al. [52]	2016	0.7130	N.A.	0.9486	0.7304
Khowaja et al. [32]	2019	0.7559	0.9758	0.9518	0.7646
Sundram et al. [4]	2019	0.7105	0.9601	0.9501	N.A
Proposed Method		0.7585	0.9846	0.9577	0.7580

show the comparison results in terms of the sensitivity ($Se$), specificity ($Sp$), accuracy ($Acc$), and $F1$ metrics for the three data sets. In all three tables, the results for the proposed method are shown in the last row. From the tables, it is evident that the proposed method is better for all metrics than all the other methods except for the $F1$ metric in Table 9. Moreover, visual evidence of the results are pictured by showing some examples for the best and the worst cases are also provided in Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. It is worth mentioning that using different FOVs may give significantly different results. So, to ensure a fair comparison, segmentation metrics are computed on pixels inside the FOV only for all methods, and FOV masks were created for STARE and CHASE_DB1 datasets as mentioned earlier.

Table 7 shows the quantitative results of the proposed method with comparison to the relevant state-of-the-art methods on the DRIVE dataset.

Table 7, Table 8 and Table 9 show an overview of segmentation results for the DRIVE, STARE and CHASE_DB1 datasets. The performance results of the proposed method in terms of average accuracy ($Acc$), sensitivity ($Se$), specificity ($Sp$), and $F1$ score outperform the seven state-of-the-art methods.

A sample of the visual results of the best and worst results for the three datasets in two different images are shown in Figure 14 and Figure 15 for DRIVE, Figure 16 and Figure 17 for the STARE dataset, and Figure 18 and Figure 19 for the CHASE_DB1 dataset. However, this method shows bad performance for CHASE_DB1 compared with DRIVE and STARE where a small decreasing in performance (associated with Table 9) is noticed for all results and for $F1$ in particular. This may be due to the high dimensions of this dataset, the characteristic of the nonuniform background illumination for the images, the bad contrast between vessels and background in addition to the presence of central vessel reflex.

Table 8 shows the quantitative results of the proposed method with comparison to the relevant state-of-the-art-methods on the STARE dataset.

A sample of the visual results of the best and worst results on the STARE dataset in two different images is shown in Figure 16 and Figure 17.

Table 9 shows the quantitative results of the proposed method with comparison to the relevant state-of-the-art-methods on the CHASE_DB1 dataset.

A sample of the visual results of the best and worst results on the CHASE_DB1 dataset in two different images is shown in Figure 18 and Figure 19.

The evaluation performance of the proposed method is computed on pixels inside the FOV only. As with most of the previous work in the DRIVE dataset, the FOV mask is used, while the FOV for both STARE and CHASE_DB1 are created by the same method as that used in [11,15,26,30].

It worth mentioning that using different FOVs may give significantly different results. To compare our results to other relative retinal vessel segmentation algorithms, the results are compared to those of state-of-the-art methods for DRIVE, STARE, and CHASE_DB1 datasets, respectively. The comparison results show the robustness of the proposed method in improving performance measure values of accuracy, sensitivity, and specificity to all datasets used in this study.

4. Discussion

4.1. The Performance Model

While experimenting, we observed that the proposed method performs better than the conventional supervised and faster than the automatic feature extracted models. We proposed an automatic system to extract and differentiate between the blood vessel and background pixels in the fundus image. A modified MLP neural network was employed to act as a semi-deep network. This network is strengthened with the use of Relu, which is used mainly in convolution neural networks (CNN) and CE loss function. CE as a loss function is the most proper to use due to the nature of the distribution for the vessel and background pixels, which is not a Gaussian prior. However, it is independent and identically distributed (I.I.D). The modified MLP learns from the 24 hand-crafted features how to successfully detect and differentiate between the blood vessel and the background pixel. The proposed method profoundly depends on the variety between the type of the selected features and the design of the network itself. In other words, we take our account to make some symmetry for features itself to form a feature vector that can represent each pixel adequately. Despite that preprocessing is always the first stage in most of the state-of-the-art-methods in vessels segmentation model, we preferred to offset it by using five groups of features that can homogenise, uniform, and process the noise in the image intensity. For example, the blood vessel pixel has a blob like-edge, which is detected effectively using the maximum and minimum moment of phase congruency features. The second derivative of hessian matrix features is used to identify the zero-crossing in the final images, the five Gaussian scaled features successfully rendered a more homogeneous image, white top-hat features corrected the uniformity of the image, black bottom-hat features suppressed the noise in the image. Because of the difference in the image dimensions between datasets, each dataset was independently evaluated and tested. We used only one classifier, MLP, and the final results depend on it without any other procedure of features selection or any ensemble or hierarchy classifier. We represented some gaps in the introduction, which were mainly in the complicated structures of the system or the high time cost to segment the vessel. We described some gaps, which were primarily focal to the complicated nature of the system and the high time cost. The datasets have many challenges, such as the scarcity, number of available images, and the variety of image illumination and dimensions. Hence, for more generalised DR diagnosis, more retinal fundus images with high illumination are necessary.

4.2. The Medical Implications of the Model

A modified MLP model acts as a soft and rapid automation system to provide the ophthalmology with primary information and knowledge about the blood vessels (e.g., size, tortuosity, crossing, lesions structure, hard and soft exudates, cotton-wool-spots) and the quick aids. This system can not uproot the role of the ophthalmology doctors, but it can help in understanding the diagnosis results, early diagnosis, and increase the result accuracy as well. Despite getting good results in the proposed model, there is still a chance for further improvement.

5. Conclusions

In this paper, we have proposed an integrated method to segment blood vessels in retinal images. The experimental results prove that the method succeeds both absolutely and in comparison with seven other state-of-the-art similar methods using three well known publicly available datasets. The proposed method encompasses many elements that wholly contribute to its success. The first element is a feature choice, where 24 carefully selected features are employed. The second element is the choice of the classifier, namely an MLP neural network, and the design of this classifier. The third element is the way the classifier is trained, where randomly chosen images from all three datasets participate in the training. The fourth element is post-processing, where classification artefacts can be mitigated. Classification results are provided both visually and quantitatively, wherein the latter case, five standard classification metrics are employed. Both the visual and quantitative results show the good performance of the proposed method vividly, despite the wide variation of vessel sizes. Actually, the method detects both wide and tiny blood vessel with the same efficiency. We avoid using preprocessing to preserve the blood vessel structure and avoid losing any information. Furthermore, our results are only based on the crafty choice of features and the correcting of the outlier values by normalisation. So concentrating on the handling of the outlier values by the feature engineer may help in decreasing the reliability of using feature selection techniques or hierarchical classification as a way to increase system performance. Our experimental results contradict previous thought that the only advantage is in using feature selection and/or hierarchy classification techniques to improve performance. On the contrary, it suggests that using features that are selected craftily or by using feature engineers can give promising results. In the retinal disease diagnosis community, the relevant dataset is scarce. In other words, the availability for investigating and researching using real images from patients who suffer from an eye ailment is a hard task, particularly in a developing country. So, the only available choice is to work on the most publicity available dataset. The average accuracy, sensitivity and specificity values accomplished by the proposed method are 0.9607, 0.7542, 0.9843 on DRIVE, 0.9634, 0.7806, and 0.9825 on STARE, and 0.9577, 0.7585, and 0.9846 on CHASE_DB1, which are comparable with the current blood vessel segmentation techniques. Although 24-D feature vectors are constructed for every pixel in the image, the ensuing space and time complexity is offset by the simplicity of the method, compared to the cutting edge methods like deep learning. The fast implementation and simplicity of the proposed method make it a handy early prediction tool of DR.