# Multichannel Retinal Blood Vessel Segmentation Based on the Combination of Matched Filter and U-Net Network

Authors: Yuliang Ma; Zhenbin Zhu; Zhekang Dong; Tao Shen; Mingxu Sun; Wanzeng Kong
Journal: BioMed Research International
 (2021)

---

## Abstract

Aiming at the current problem of insufficient extraction of small retinal blood vessels, we propose a retinal blood vessel segmentation algorithm that combines supervised learning and unsupervised learning algorithms. In this study, we use a multiscale matched filter with vessel enhancement capability and a U-Net model with a coding and decoding network structure. Three channels are used to extract vessel features separately, and finally, the segmentation results of the three channels are merged. The algorithm proposed in this paper has been verified and evaluated on the DRIVE, STARE, and CHASE_DB1 datasets. The experimental results show that the proposed algorithm can segment small blood vessels better than most other methods. We conclude that our algorithm has reached 0.8745, 0.8903, and 0.8916 on the three datasets in the sensitivity metric, respectively, which is nearly 0.1 higher than other existing methods.


---

## Body

## 1. Introduction

The human eyes consist of the following parts: cornea, pupil, iris, vitreous, and retina. Abnormalities in any of these tissue structures may cause vision defects or even blindness. Among them, the study of retinal structure and its blood vessels is significant [1]. The extraction of retinal blood vessels and the characterization of morphological properties, such as diameter, shape, distortion, and bifurcation, can be used to screen, evaluate, and treat different ocular abnormalities [2]. Evaluation of retinal vascular properties, such as changes in width, is used to analyze hypertension, while bifurcation points and tortuosity can help identify cardiovascular disease and diabetic retinopathy [3].The retinal vessel extraction methods, including pattern recognition, are classified into five core classes [4]. The pattern recognition techniques are generally divided into two categories: supervised learning and unsupervised learning. The supervised learning method needs to use manual segmentation images of ophthalmologists for training. This method requires many training images, and the training time is longer than that of other methods, but this method has an excellent generalized effect and can be applied to other images of the same type. Compared with supervised learning, nonsupervised learning methods, such as matched filtering, mathematical morphology operations, blood vessel tracking, and clustering, do not require corresponding image labels but analyze and process based on the existing data. These two types of methods have been applied and innovated by many researchers in recent years.


### 1.1. Unsupervised Learning Methods

Literature [5] proposed a new kernel-based technique, viz, Fréchet PDF-based matched filter. The new method performs a better matching between the vessel profile and Fréchet template. Literature [6] improved the extraction method of blood vessels, using a series of morphological operations to extract small blood vessels, and finally fused with the segmented image to supplement the small blood vessels. Compared with other algorithms, it can segment as many tiny blood vessels as possible. However, the steps of the algorithm are too complicated, and although the final segmentation effect obtains the smallest blood vessels, the small blood vessels are in an intermittent state as a whole, and they are not well connected with thicker blood vessels. Literature [7] proposed a new matched filtering method, which applies contrast-limited adaptive histogram equalization and Gaussian second-derivative-based matched filter in preprocessing and uses an entropy-based optimal threshold method performing binarization. This algorithm effectively improves the sensitivity metric of segmentation, but like literature [6], it does not perform well with accuracy. Literature [8] proposed an automatic segmentation method of retinal blood vessels using a matched filter and fuzzy C-means clustering. The algorithm uses contrast-limited adaptive histogram equalization to enhance the contrast of the image. After using Gabor and Frangi filters for noise removal and background removal, the fuzzy C-means are used to extract the initial vascular network, and the integrated level set method is used to refine segmentation further. The algorithm has good sensitivity and specificity. The problem is that the ability to segment small blood vessels is limited, and many segmentation details are missed. Literature [9] proposed a novel method to extract the retinal blood vessel using local contrast normalization and a second-order detector. The proposed methodology achieves higher accuracy in vessel segmentation than existing techniques. Literature [10] proposed a novel matched filter approach with the Gumbel probability distribution function as its kernel. The reason to achieve the higher accuracy is due to a better matching filter with the Gumbel PDF-based kernel.


### 1.2. Supervised Learning Methods

Literature [11] proposed a method using deep conventional neural networks and a hysteresis threshold method to detect the vessels accurately. The proposed method gives good performance in which more tiny vessels are detected. Literature [12] proposed a multilevel CNN model applied for automatic blood vessel segmentation in retinal fundus images. A novel max-resizing technique is proposed to improve the generalization of the training procedure for predicting blood vessels from retinal fundus images. Literature [13] proposed a new segment-level loss used with the pixel-wise loss to balance the importance between thick vessels and thin vessels in the training process. Literature [14] proposed a cross-connected convolutional neural network (CcNet) to automatically segment retinal vessel trees. The cross connections between a primary path and a secondary path fuse the multilevel features. This method has relatively advanced performances, including competitive strong robustness and segmentation speed. Literature [15] proposed a method for retinal vessel segmentation using patch-based fully convolutional networks. Literature [16] applied dilated convolutions in a deep neural network to improve the segmentation of retinal blood vessels from fundus images. Literature [17] proposed a new improved algorithm based on the U-Net network model. The algorithm integrates the Inception-Res structure module and the Dense-Inception structure module into the U-Net structure. The algorithm dramatically deepens the depth of the network but does not add additional training parameters. It has good segmentation performance in the image segmentation of retinal blood vessels and has strong generalization ability. Literature [18] proposed a new hybrid algorithm for retinal vessel segmentation on fundus images. The proposed algorithm applies a new directionally sensitive blood vessel enhancement before sending fundus images to U-Net. Literature [19] proposed a supervised method based on a pretrained fully convolutional network through transfer learning. This method simplifies the typical retinal vessel segmentation problem into regional semantic vessel element segmentation tasks. Generally, unsupervised methods are less complex and suffer from relatively lower accuracy than supervised methods [13].To solve the problem of insufficient segmentation of small blood vessels in most papers, we have devised a new automatic segmentation framework for retinal vessels based on improving U-Net and a multiscale matched filter. The creative points of this paper are summarized as follows:(1)
We proposed an improved black hat algorithm to enhance the characteristics of blood vessels and reduce the interference of other tissues(2)
An algorithm combining a multiscale matched filter and U-Net neural network is proposed. This paper mainly uses the improved U-Net convolutional neural network combined with a multiscale matched filter to perform multichannel blood vessel segmentation processing on the retinal fundus image(3)
We have devised a new loss function to train the improved U-Net neural network to solve pixel imbalance in the image betterThe rest of this paper is organized as follows. Section2 outlines the proposed method and datasets. The performance of the proposed method and the discussion are described in detail in Section 3. A conclusion is drawn in Section 4.


## 1.1. Unsupervised Learning Methods

Literature [5] proposed a new kernel-based technique, viz, Fréchet PDF-based matched filter. The new method performs a better matching between the vessel profile and Fréchet template. Literature [6] improved the extraction method of blood vessels, using a series of morphological operations to extract small blood vessels, and finally fused with the segmented image to supplement the small blood vessels. Compared with other algorithms, it can segment as many tiny blood vessels as possible. However, the steps of the algorithm are too complicated, and although the final segmentation effect obtains the smallest blood vessels, the small blood vessels are in an intermittent state as a whole, and they are not well connected with thicker blood vessels. Literature [7] proposed a new matched filtering method, which applies contrast-limited adaptive histogram equalization and Gaussian second-derivative-based matched filter in preprocessing and uses an entropy-based optimal threshold method performing binarization. This algorithm effectively improves the sensitivity metric of segmentation, but like literature [6], it does not perform well with accuracy. Literature [8] proposed an automatic segmentation method of retinal blood vessels using a matched filter and fuzzy C-means clustering. The algorithm uses contrast-limited adaptive histogram equalization to enhance the contrast of the image. After using Gabor and Frangi filters for noise removal and background removal, the fuzzy C-means are used to extract the initial vascular network, and the integrated level set method is used to refine segmentation further. The algorithm has good sensitivity and specificity. The problem is that the ability to segment small blood vessels is limited, and many segmentation details are missed. Literature [9] proposed a novel method to extract the retinal blood vessel using local contrast normalization and a second-order detector. The proposed methodology achieves higher accuracy in vessel segmentation than existing techniques. Literature [10] proposed a novel matched filter approach with the Gumbel probability distribution function as its kernel. The reason to achieve the higher accuracy is due to a better matching filter with the Gumbel PDF-based kernel.


## 1.2. Supervised Learning Methods

Literature [11] proposed a method using deep conventional neural networks and a hysteresis threshold method to detect the vessels accurately. The proposed method gives good performance in which more tiny vessels are detected. Literature [12] proposed a multilevel CNN model applied for automatic blood vessel segmentation in retinal fundus images. A novel max-resizing technique is proposed to improve the generalization of the training procedure for predicting blood vessels from retinal fundus images. Literature [13] proposed a new segment-level loss used with the pixel-wise loss to balance the importance between thick vessels and thin vessels in the training process. Literature [14] proposed a cross-connected convolutional neural network (CcNet) to automatically segment retinal vessel trees. The cross connections between a primary path and a secondary path fuse the multilevel features. This method has relatively advanced performances, including competitive strong robustness and segmentation speed. Literature [15] proposed a method for retinal vessel segmentation using patch-based fully convolutional networks. Literature [16] applied dilated convolutions in a deep neural network to improve the segmentation of retinal blood vessels from fundus images. Literature [17] proposed a new improved algorithm based on the U-Net network model. The algorithm integrates the Inception-Res structure module and the Dense-Inception structure module into the U-Net structure. The algorithm dramatically deepens the depth of the network but does not add additional training parameters. It has good segmentation performance in the image segmentation of retinal blood vessels and has strong generalization ability. Literature [18] proposed a new hybrid algorithm for retinal vessel segmentation on fundus images. The proposed algorithm applies a new directionally sensitive blood vessel enhancement before sending fundus images to U-Net. Literature [19] proposed a supervised method based on a pretrained fully convolutional network through transfer learning. This method simplifies the typical retinal vessel segmentation problem into regional semantic vessel element segmentation tasks. Generally, unsupervised methods are less complex and suffer from relatively lower accuracy than supervised methods [13].To solve the problem of insufficient segmentation of small blood vessels in most papers, we have devised a new automatic segmentation framework for retinal vessels based on improving U-Net and a multiscale matched filter. The creative points of this paper are summarized as follows:(1)
We proposed an improved black hat algorithm to enhance the characteristics of blood vessels and reduce the interference of other tissues(2)
An algorithm combining a multiscale matched filter and U-Net neural network is proposed. This paper mainly uses the improved U-Net convolutional neural network combined with a multiscale matched filter to perform multichannel blood vessel segmentation processing on the retinal fundus image(3)
We have devised a new loss function to train the improved U-Net neural network to solve pixel imbalance in the image betterThe rest of this paper is organized as follows. Section2 outlines the proposed method and datasets. The performance of the proposed method and the discussion are described in detail in Section 3. A conclusion is drawn in Section 4.


## 2. Materials and Methods

### 2.1. System Overview

The proposed algorithm consists of three steps: preprocessing datasets, training U-Net in 3 channels, and postprocessing. This algorithm’s main feature extraction framework is based on the improved U-Net model, using three feature extraction channels. It is mainly to perform a whole feature extraction of the image in channel 1 so that some morphological operations are performed in the preprocessing part to reduce image artifacts and noise. On the remaining two channels, matched filters are used to extract retinal vessels of different scales, and then, the improved U-Net model is used to extract features, and the OR-type operator is used to fuse the final output image. Experimental results verify that the image processed by multichannel matched filtering is better than the unprocessed image. The overall flowchart is shown in Figure1.Figure 1
Proposed model.


### 2.2. Datasets

To verify the effectiveness of the algorithm in this paper, this paper chooses three commonly used public datasets for training and testing: DRIVE, STARE, and CHASE_DB1 datasets. These datasets include a wide range of challenging images. The DRIVE contains 40 color retinal fundus images divided into a training set and a testing set. The plane resolution of DRIVE is565×584. The STARE contains 20 color retinal fundus images with a resolution of 605×700 pixels. Unlike the DRIVE, this dataset does not have a training set and a testing set. The CHASE_DB1 contains 28 color retinal fundus images with a resolution of 960×999 pixels, and the training set and testing set are also not divided. Each image in these three datasets has a label of retinal blood vessel image segmented manually by two professional physicians. We randomly selected 5 images in the STARE dataset as test images (im0002, im0077, im0163, im0255, and im0291), and the remaining 15 images were set as the training set. In CHASE_DB1, we select the last 8 images as the test set and the remaining 20 images as the training set. Note that mask images of STARE and CHASE_DB1 are not available, so we extracted the green channel of the images and then used some morphological algorithms and threshold algorithm to obtain the mask images.


### 2.3. Preprocessing

In this paper, the green channel is selected as the input image of the preprocessing part. This is because the retinal blood vessels presented by the green channel have better contrast with the background compared with the red channel and the blue channel [20, 21], as shown in Figure 2.Figure 2
Color fundus image and its different RGB channels: (a) original RGB image; (b) red channel; (c) green channel; (d) blue channel.
(a)(b)(c)(d)It can be seen from Figure2 that the appearance of blood vessels on the green channel of the color image consists of more information compared to that on the red and blue channel images, but the overall image is still dark, and the contrast is not obvious. In order to improve this situation, adaptive histogram threshold processing (CLAHE) [22] and gamma transformation are performed on the extracted green channel grayscale image, as shown in Figure 3. In this part of the process, CLAHE is used to enhance the contrast between the nonvessels and blood vessels, and gamma transformation is used to adjust and reduce the background noise in the image. We can see Tables 1–3 in Supplementary Materials for a comprehensive comparison of blood vessel enhanced algorithms, and these data can prove that the CLAHE method improves the general performance of the proposed method.Figure 3
Typical images after each preprocessing step: (a) original RGB image; (b) red channel; (c) image after CLAHE operation; (d) image after gamma correction.
(a)(b)(c)(d)


### 2.4. Multichannel Feature Extraction

#### 2.4.1. Channel 1

In order to retain all the blood vessel feature information of the image as much as possible, some morphological operations are used in channel 1 to remove background noise, and then, the U-Net network is used for feature extraction. For the artifacts caused by uneven illumination in the image and nonvascular structures, we use the morphological closing operation algorithm to estimate the background and then perform the result using the mathematical operation shown in equation (1).It can be seen intuitively from Figure4 that the brighter video disc structure in the original image is removed, and most of the artifacts are also processed.(1)gx,y=255−Iclosex,y−Ix,y+1m∗n∑x=1m∑y=1nIclosex,y,fx,y=255maxgx,y−mingx,y∗gx,y−mingx,y,where fx,y is the processed image and Iclosex,y is the image after a morphological closing operation. We select disk type structuring elements for the closing operator having a radius of eleven pixels. Ix,y is the original image; m and n are the image pixel size.Figure 4
(a, e, i) The fundus images, (b, f, j) green channel images applying CLAHE and gamma transformation, (c, g, k) background extracted by close operation, and (d, h, l) the final results. The different samples of (a–d) DRIVE, (e–h) STARE, and (i–l) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


#### 2.4.2. Channel 2

By analyzing the gray image of retinal blood vessels, it can be found that the cross-sectional gray intensity of blood vessels is distributed in an inverted Gaussian curve, the gray value of the center line of the blood vessel is low, and the gray value at the edge of the blood vessel is high [5]. Aiming at this remarkable feature of retinal blood vessel images, Chaudhuri et al. [23] designed a Gaussian matched filter and used its distribution to simulate the grayscale intensity distribution of blood vessel cross sections and filter the blood vessels in sections. In this paper, the matched filters are used in channel 2 and channel 3 to separately enhance and extract the large and small blood vessels to realize the comprehensive segmentation of retinal blood vessels.Define the two-dimensional Gaussian kernel function as(2)Kx,y=−e−x2/2s2,y≤l2,where s is the width of the Gaussian kernel and l is the length of the Gaussian kernel. The blood vessel starts from the center of the optic disc and extends in multiple directions. Rotating the Gaussian kernel is used to filter the multidirectional blood vessels.Assuming thatpx,y is a discrete point in the kernel function, the rotation matrix is(3)gi=cosθi−sinθisinθicosθi.θi0≤θi≤p is the angle of the i-th kernel function, and the coordinate value of px,y after rotation is p¯i=u,v; then, the i-th template kernel function is(4)Kix,y=−e−u2/2s2,∀p¯i∈N,where N is the template field, and the value range is(5)N=u,v,u≤3s,v≤l2.In actual algorithm applications, it is often necessary to consider the mean value of the correlation coefficient of the template filter, as shown in(6)mi=∑p¯i∈NKix,yA.Among them,A represents the number of points in the template area. So, the final template kernel function is(7)Ki′x,y=Kix,y−mi,∀p¯i∈N.This paper improves and optimizes the dependence of Gaussian matched filter response on a vessel diameter. The image enhancement result using large-scale Gaussian matched filtering in channel 2 is shown in Figure5, where the parameters are set to l=10.8, s=1.9, and 8 directions which means i=1,2,⋯,8 in equation (3). It can be seen from the image that the algorithm has a better segmentation effect for thicker blood vessels and strong antinoise, but it has a poor segmentation effect on small blood vessels, and there is a problem that the smaller blood vessels cannot be distinguished from the background, and the blood vessels are easily broken. In order to solve this problem, this paper proposes an improved method based on the black hat algorithm, which can effectively reduce the influence of background noise by subtracting the original image before matching filter processing and the obtained image after processing to enhance the characteristics of blood vessels. We performed a series of processing transformations as shown in equations (8) and (9) on the images processed by large-scale matched filtering. We call this algorithm black hat2.(8)Bhatf=fx,y∙bu,v−fx,y,(9)gx,y=255−fx,y−2∗Bhatf,where ∙ is the morphological closing operation and bu,v is disk type structuring element, Bhatf is the black hat transformation, fx,y is the original image, and gx,y is the final processed image.Figure 5
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) large-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


#### 2.4.3. Channel 3

This paper uses a small-scale Gaussian matched filter to enhance the image of small blood vessels, as shown in Figure6. After many experiments, the parameters of the matched filter are set as l=5, s=0.1, and 18 directions which means i=1,2,⋯,18 in equation (3). Using small-scale filters can effectively enhance the small blood vessels in the image, but at the same time, it also enhances much striped noise in the image, and the enhancing effect on the thick blood vessels with central reflection is poor. To reduce the background noise, the black hat2 algorithm used in channel 2 is also used to remove the background in channel 3.Figure 6
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) small-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


### 2.5. U-Net Model

In image semantic segmentation using deep learning, the U-Net network model is the most widely used, which is improved based on the classic full convolutional network (FCN) [24]. U-Net is an image-to-image pixel-level classification network, and its network structure is apparent, as shown in Figure 7. U-Net is different from other standard segmentation networks: U-Net uses an entirely different feature fusion method—splicing. U-Net stitches the features together in the channel dimension. This method fuses the in-depth features extracted from the image with the shallow features to form thicker features, while the fusion operation of FCN only uses corresponding point addition and does not obtain thicker features.Figure 7
U-Net model architecture.Unlike the structure in the original literature [24], this paper sets the padding value of 1 in each layer’s convolution operation, and the convolution kernel size is 3∗3. The purpose is to ensure that the output and input image sizes are consistent and avoid the size increasing operation in the output layer. It is essentially a binary classification operation in the output layer of U-Net. We use an adaptive threshold segmentation algorithm for processing in this paper. The idea of this algorithm is not to calculate the global image threshold but to calculate the local threshold according to different areas of the image, so for different areas of the image, the algorithm can adaptively calculate different thresholds and perform binary segmentation. The specific calculation process is shown in(10)T=−b+12m+1×2n+1∑i=0n∑j=0mgx±i,y±j,where b is the fixed parameter, 2m+1×2n+1 is the area, and T is the area’s threshold.This paper proposes a new loss function that combines the Dice coefficient with the two-class cross-entropy loss function. The Dice coefficient is widely used in the evaluation of image segmentation. In order to facilitate the formation of the minimized loss function, as shown in(11)Ldice=1−2X∩YX+Y,where X∩Y represents the common elements of the prediction graph and the label graph, X and Y represent the number of elements of the prediction graph and the label. In order to facilitate the calculation, approximate ∣X∩Y∣ as the dot product between the predicted probability map and the label, and add the elements in the result. ∣X∣ and ∣Y∣ are quantified by summing the squares of each element. As shown in(12)Ldice=1−2∑iNpk,iqk,i∑iNp2k,i+∑iNq2k,i,where N is the number of pixels, pk,i∈0,1 and qk,i∈0,1 are the predicted probabilities and true labels of the pixel belonging to category k.The cross-entropy loss function used to optimize the network is shown as(13)Lr=−∑iN1−TPNpylogp+1−TNNn1−ylog1−p,where TP and TN are the numbers of true positive and true negative pixels, respectively; Np and Nn are the numbers of segmented pixels and nonsegmented pixels, respectively; y is the label value (y=1, segmentation target; y=0, background); and p is the predicted probability value of the pixel.A coefficientλ is introduced to define the new loss function Loss, as shown in(14)Loss=Ldice+λLr.Notably, the coefficientλ is set to 0.5 in this work, and the flowchart of U-Net is summarized in Algorithm 1.Algorithm 1: Training of U-Net with dynamic learning rate.
Input: Train imagesX, ground truth GInput: Initial epochs E⟵ 30, batch size ⟵ 1, learning rate lr⟵ 0.01Input: Initialize best loss bl⟵ InfOutput: Predicted images P, U-Net parameter1.Xpre⟵ preprocessing (X)2.Xenh⟵ enhancement (Xpre)3.fore⟵ 0 to Edo4.ife=1/3∗Ethen5.lr⟵0.1∗lr6.else ife=2/3∗Ethen7.lr⟵0.1∗lr8.end if9.N⟵ compute the number of train imagesX10.par⟵ initial parameter of U-Net11.whilen⟵1<N+1do12.Yn⟵Unet(Xenhn,par)13.lossn⟵DiceYn,Gn+λCross_entropyYn,Gn14.iflossn<blthen15.bl⟵lossn16.par⟵newpar17.end if18.par⟵SGDpar,lr19.end while20.end for21.P⟵UnetXenh,par22.returnP, par


### 2.6. Postprocessing

In the postprocessing, since the final segmentation image merges the three segmentation images, the noise in the resulting image is also superimposed on all the noises of the three images. Noises will undoubtedly have a significant impact on the actual effect of the segmented image, so this paper addresses this issue in the final postprocessing step. In this paper, a morphological algorithm is used to calculate the size of the connected area of the image. The 8-adjacent connection method is adopted to eliminate the area with the connected area less than 25 pixels, which is to reclassify the area pixels as background. This paper selects a test image in the DRIVE dataset for experimental comparison, and the comparison images are shown in Figure8.Figure 8
Comparison of postprocessing: (a–e) segmentation image without postprocessing; (f–j) segmentation image applying postprocessing.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


### 2.7. Experimental Design

#### 2.7.1. U-Net Implementation Details

The U-Net model used in this paper is slightly different from the structure in literature [24]. In order to keep the input and output image sizes of the model consistent, the convolution structure is adjusted accordingly. The specific model structure parameters are shown in Table 1.Table 1
The parameters of the U-Net architecture.
Block nameLayer nameImage sizeParametersDoubleConvConv (ksize=3, pad=1)3∗3∗C1+1∗C2BN+LReLU2∗C2Conv (ksize=3, pad=1)3∗3∗C2+1∗C2BN+LReLU2∗C29∗C1+C2∗C2+6∗C2Input1×576×5760Encoder block_1DoubleConv_164×576×5769∗1+64∗64+6∗64=37824Maxpooling (ksize=2)64×576×5760Encoder block_2DoubleConv_2128×288×2889∗64+128∗128+6∗128=221952Maxpooling (ksize=2)128×288×2880Encoder block_3DoubleConv_3256×144×1449∗128+256∗256+6∗256=886272Maxpooling (ksize=2)256×144×1440Encoder block_4DoubleConv_4512×72×729∗256+512∗512+6∗512=3542016Maxpooling (ksize=2)512×72×720Encoder block_5DoubleConv_5512×36×369∗512+512∗512+6∗512=4721664Maxpooling (ksize=2)512×36×360Decoder block_1Upsampling (bilinear)512×72×720Concat1024×72×720DoubleConv_6256×72×729∗1024+256∗256+6∗256=2950656Decoder block_2Upsampling (bilinear)256×144×1440Concat512×144×1440DoubleConv_7128×144×1449∗512+128∗128+6∗128=738048Decoder block_3Upsampling (bilinear)128×288×2880Concat256×288×2880DoubleConv_864×288×2889∗256+64∗64+6∗64=184704Decoder block_4Upsampling (bilinear)64×576×5760Concat128×576×5760DoubleConv_964×576×5769∗128+64∗64+6∗64=110976OutputConv (ksize=1)1×576×5761∗1∗64+1=65In training, we set the epoch to 30 and the initial learning ratelr to 0.01, and then, the learning rate is set to update in a three-stage formula, as shown in(15)lr=0.01,epoch>10,0.001,10<epoch≤20,0.0001,20<epoch≤30.Setting a larger learning rate at the beginning is to make the model obtain the vicinity of the optimal global parameters faster, and this operation can reduce the training time of the model. After training for a particular epoch, the learning rate needs to be reduced accordingly in order to make the parameters closer to the optimal value in subsequent updates. The stochastic gradient descent (SGD) algorithm is used in the optimization of the loss function.


#### 2.7.2. Training Image Preparation

We randomly select 15 images from STARE and the first 20 images from CHASE_DB1 as their respective training set. Due to the limited number of images in the existing dataset, to avoid the overfitting phenomenon in the model training, we perform data expansion processing on the training set of each dataset. Thanks to the translation invariance of the convolutional structure, the images in the training set in this paper were flipped horizontally and vertically and rotated 180 degrees to increase the amount of data 4 times.


#### 2.7.3. Measuring Metrics

In order to evaluate the segmentation performance of this algorithm, we use the following metrics to perform a comprehensive evaluation of the segmentation result. These metrics are accuracy (ACC), sensitivity (Se), specificity (Sp), and AUC and calculated as follows:(16)ACC=TP+TNTP+FN+TN+FP,(17)Se=TPTP+FN,(18)Sp=TNTN+FP,(19)AUC=12TPTP+FN+TNTN+FP,where TP is true positive, FP is false positive, TN is true negative, and FN is false negative. Se is the sensitivity, which indicates the degree of classification of blood vessels and nonvascular pixels. In this paper, higher sensitivity indicates that more tiny blood vessels can be detected. Sp is specificity, which is used to express the ability of the algorithm to recognize nonvascular pixels. ACC is the accuracy of algorithm segmentation, reflecting the gap between the algorithm segmentation result and the natural result. AUC is the area under the ROC curve, and we adopt another calculation method to get the AUC, as shown in equation (19) [11].Besides, we also use two other evaluation metrics to measure the effect of segmentation: MCC and CAL.(20)MCC=TP×TN−TP×FNTP+FP×TP+FN×TN+FP×TN+FN.MCC is a correlation coefficient between the segmentation output of the algorithm and ground truth. It comprehensively considers TP, TN, FP, and FN, which is a relatively balanced metric. Finally, it is more suitable for an imbalanced class ratio.CAL can be expressed as the product ofC, A, and L as follows:(21)fC,A,L=C×A×L.SupposeS and SG are the segmentation result and the corresponding ground truth, respectively. These functions are defined as follows:(1)
Connectivity (C): it evaluates the fragmentation degree between S and SG by comparing the number of connected components:(22)C=1−min1,#CSG−#CS#SG,where #C∙ means the number of connected components, while #∙ means the number of vessel pixels in the considered binary image.(2)
Area (A): it evaluates the degree of intersecting area between S and SG and is defined as(23)A=#δαS∩SG∪δαSG∩S#S∪SG,where δα· is a morphological dilation using a disc of α pixels in radius. We set α=2.(3)
Length (L): it evaluates the equivalent degree between S and SG by computing the total length:(24)L=#φS∩δβSG∪δβS∩φSG#φS∪φSG,where φ· is the homotopic skeletonization and δβ∙ is a morphological dilation with a disc of β pixel in radius. We set β=2.According to [26], the CAL metric is essential to quantify thick and thin vessels more equally.


## 2.1. System Overview

The proposed algorithm consists of three steps: preprocessing datasets, training U-Net in 3 channels, and postprocessing. This algorithm’s main feature extraction framework is based on the improved U-Net model, using three feature extraction channels. It is mainly to perform a whole feature extraction of the image in channel 1 so that some morphological operations are performed in the preprocessing part to reduce image artifacts and noise. On the remaining two channels, matched filters are used to extract retinal vessels of different scales, and then, the improved U-Net model is used to extract features, and the OR-type operator is used to fuse the final output image. Experimental results verify that the image processed by multichannel matched filtering is better than the unprocessed image. The overall flowchart is shown in Figure1.Figure 1
Proposed model.


## 2.2. Datasets

To verify the effectiveness of the algorithm in this paper, this paper chooses three commonly used public datasets for training and testing: DRIVE, STARE, and CHASE_DB1 datasets. These datasets include a wide range of challenging images. The DRIVE contains 40 color retinal fundus images divided into a training set and a testing set. The plane resolution of DRIVE is565×584. The STARE contains 20 color retinal fundus images with a resolution of 605×700 pixels. Unlike the DRIVE, this dataset does not have a training set and a testing set. The CHASE_DB1 contains 28 color retinal fundus images with a resolution of 960×999 pixels, and the training set and testing set are also not divided. Each image in these three datasets has a label of retinal blood vessel image segmented manually by two professional physicians. We randomly selected 5 images in the STARE dataset as test images (im0002, im0077, im0163, im0255, and im0291), and the remaining 15 images were set as the training set. In CHASE_DB1, we select the last 8 images as the test set and the remaining 20 images as the training set. Note that mask images of STARE and CHASE_DB1 are not available, so we extracted the green channel of the images and then used some morphological algorithms and threshold algorithm to obtain the mask images.


## 2.3. Preprocessing

In this paper, the green channel is selected as the input image of the preprocessing part. This is because the retinal blood vessels presented by the green channel have better contrast with the background compared with the red channel and the blue channel [20, 21], as shown in Figure 2.Figure 2
Color fundus image and its different RGB channels: (a) original RGB image; (b) red channel; (c) green channel; (d) blue channel.
(a)(b)(c)(d)It can be seen from Figure2 that the appearance of blood vessels on the green channel of the color image consists of more information compared to that on the red and blue channel images, but the overall image is still dark, and the contrast is not obvious. In order to improve this situation, adaptive histogram threshold processing (CLAHE) [22] and gamma transformation are performed on the extracted green channel grayscale image, as shown in Figure 3. In this part of the process, CLAHE is used to enhance the contrast between the nonvessels and blood vessels, and gamma transformation is used to adjust and reduce the background noise in the image. We can see Tables 1–3 in Supplementary Materials for a comprehensive comparison of blood vessel enhanced algorithms, and these data can prove that the CLAHE method improves the general performance of the proposed method.Figure 3
Typical images after each preprocessing step: (a) original RGB image; (b) red channel; (c) image after CLAHE operation; (d) image after gamma correction.
(a)(b)(c)(d)


## 2.4. Multichannel Feature Extraction

### 2.4.1. Channel 1

In order to retain all the blood vessel feature information of the image as much as possible, some morphological operations are used in channel 1 to remove background noise, and then, the U-Net network is used for feature extraction. For the artifacts caused by uneven illumination in the image and nonvascular structures, we use the morphological closing operation algorithm to estimate the background and then perform the result using the mathematical operation shown in equation (1).It can be seen intuitively from Figure4 that the brighter video disc structure in the original image is removed, and most of the artifacts are also processed.(1)gx,y=255−Iclosex,y−Ix,y+1m∗n∑x=1m∑y=1nIclosex,y,fx,y=255maxgx,y−mingx,y∗gx,y−mingx,y,where fx,y is the processed image and Iclosex,y is the image after a morphological closing operation. We select disk type structuring elements for the closing operator having a radius of eleven pixels. Ix,y is the original image; m and n are the image pixel size.Figure 4
(a, e, i) The fundus images, (b, f, j) green channel images applying CLAHE and gamma transformation, (c, g, k) background extracted by close operation, and (d, h, l) the final results. The different samples of (a–d) DRIVE, (e–h) STARE, and (i–l) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


### 2.4.2. Channel 2

By analyzing the gray image of retinal blood vessels, it can be found that the cross-sectional gray intensity of blood vessels is distributed in an inverted Gaussian curve, the gray value of the center line of the blood vessel is low, and the gray value at the edge of the blood vessel is high [5]. Aiming at this remarkable feature of retinal blood vessel images, Chaudhuri et al. [23] designed a Gaussian matched filter and used its distribution to simulate the grayscale intensity distribution of blood vessel cross sections and filter the blood vessels in sections. In this paper, the matched filters are used in channel 2 and channel 3 to separately enhance and extract the large and small blood vessels to realize the comprehensive segmentation of retinal blood vessels.Define the two-dimensional Gaussian kernel function as(2)Kx,y=−e−x2/2s2,y≤l2,where s is the width of the Gaussian kernel and l is the length of the Gaussian kernel. The blood vessel starts from the center of the optic disc and extends in multiple directions. Rotating the Gaussian kernel is used to filter the multidirectional blood vessels.Assuming thatpx,y is a discrete point in the kernel function, the rotation matrix is(3)gi=cosθi−sinθisinθicosθi.θi0≤θi≤p is the angle of the i-th kernel function, and the coordinate value of px,y after rotation is p¯i=u,v; then, the i-th template kernel function is(4)Kix,y=−e−u2/2s2,∀p¯i∈N,where N is the template field, and the value range is(5)N=u,v,u≤3s,v≤l2.In actual algorithm applications, it is often necessary to consider the mean value of the correlation coefficient of the template filter, as shown in(6)mi=∑p¯i∈NKix,yA.Among them,A represents the number of points in the template area. So, the final template kernel function is(7)Ki′x,y=Kix,y−mi,∀p¯i∈N.This paper improves and optimizes the dependence of Gaussian matched filter response on a vessel diameter. The image enhancement result using large-scale Gaussian matched filtering in channel 2 is shown in Figure5, where the parameters are set to l=10.8, s=1.9, and 8 directions which means i=1,2,⋯,8 in equation (3). It can be seen from the image that the algorithm has a better segmentation effect for thicker blood vessels and strong antinoise, but it has a poor segmentation effect on small blood vessels, and there is a problem that the smaller blood vessels cannot be distinguished from the background, and the blood vessels are easily broken. In order to solve this problem, this paper proposes an improved method based on the black hat algorithm, which can effectively reduce the influence of background noise by subtracting the original image before matching filter processing and the obtained image after processing to enhance the characteristics of blood vessels. We performed a series of processing transformations as shown in equations (8) and (9) on the images processed by large-scale matched filtering. We call this algorithm black hat2.(8)Bhatf=fx,y∙bu,v−fx,y,(9)gx,y=255−fx,y−2∗Bhatf,where ∙ is the morphological closing operation and bu,v is disk type structuring element, Bhatf is the black hat transformation, fx,y is the original image, and gx,y is the final processed image.Figure 5
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) large-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


### 2.4.3. Channel 3

This paper uses a small-scale Gaussian matched filter to enhance the image of small blood vessels, as shown in Figure6. After many experiments, the parameters of the matched filter are set as l=5, s=0.1, and 18 directions which means i=1,2,⋯,18 in equation (3). Using small-scale filters can effectively enhance the small blood vessels in the image, but at the same time, it also enhances much striped noise in the image, and the enhancing effect on the thick blood vessels with central reflection is poor. To reduce the background noise, the black hat2 algorithm used in channel 2 is also used to remove the background in channel 3.Figure 6
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) small-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 2.4.1. Channel 1

In order to retain all the blood vessel feature information of the image as much as possible, some morphological operations are used in channel 1 to remove background noise, and then, the U-Net network is used for feature extraction. For the artifacts caused by uneven illumination in the image and nonvascular structures, we use the morphological closing operation algorithm to estimate the background and then perform the result using the mathematical operation shown in equation (1).It can be seen intuitively from Figure4 that the brighter video disc structure in the original image is removed, and most of the artifacts are also processed.(1)gx,y=255−Iclosex,y−Ix,y+1m∗n∑x=1m∑y=1nIclosex,y,fx,y=255maxgx,y−mingx,y∗gx,y−mingx,y,where fx,y is the processed image and Iclosex,y is the image after a morphological closing operation. We select disk type structuring elements for the closing operator having a radius of eleven pixels. Ix,y is the original image; m and n are the image pixel size.Figure 4
(a, e, i) The fundus images, (b, f, j) green channel images applying CLAHE and gamma transformation, (c, g, k) background extracted by close operation, and (d, h, l) the final results. The different samples of (a–d) DRIVE, (e–h) STARE, and (i–l) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 2.4.2. Channel 2

By analyzing the gray image of retinal blood vessels, it can be found that the cross-sectional gray intensity of blood vessels is distributed in an inverted Gaussian curve, the gray value of the center line of the blood vessel is low, and the gray value at the edge of the blood vessel is high [5]. Aiming at this remarkable feature of retinal blood vessel images, Chaudhuri et al. [23] designed a Gaussian matched filter and used its distribution to simulate the grayscale intensity distribution of blood vessel cross sections and filter the blood vessels in sections. In this paper, the matched filters are used in channel 2 and channel 3 to separately enhance and extract the large and small blood vessels to realize the comprehensive segmentation of retinal blood vessels.Define the two-dimensional Gaussian kernel function as(2)Kx,y=−e−x2/2s2,y≤l2,where s is the width of the Gaussian kernel and l is the length of the Gaussian kernel. The blood vessel starts from the center of the optic disc and extends in multiple directions. Rotating the Gaussian kernel is used to filter the multidirectional blood vessels.Assuming thatpx,y is a discrete point in the kernel function, the rotation matrix is(3)gi=cosθi−sinθisinθicosθi.θi0≤θi≤p is the angle of the i-th kernel function, and the coordinate value of px,y after rotation is p¯i=u,v; then, the i-th template kernel function is(4)Kix,y=−e−u2/2s2,∀p¯i∈N,where N is the template field, and the value range is(5)N=u,v,u≤3s,v≤l2.In actual algorithm applications, it is often necessary to consider the mean value of the correlation coefficient of the template filter, as shown in(6)mi=∑p¯i∈NKix,yA.Among them,A represents the number of points in the template area. So, the final template kernel function is(7)Ki′x,y=Kix,y−mi,∀p¯i∈N.This paper improves and optimizes the dependence of Gaussian matched filter response on a vessel diameter. The image enhancement result using large-scale Gaussian matched filtering in channel 2 is shown in Figure5, where the parameters are set to l=10.8, s=1.9, and 8 directions which means i=1,2,⋯,8 in equation (3). It can be seen from the image that the algorithm has a better segmentation effect for thicker blood vessels and strong antinoise, but it has a poor segmentation effect on small blood vessels, and there is a problem that the smaller blood vessels cannot be distinguished from the background, and the blood vessels are easily broken. In order to solve this problem, this paper proposes an improved method based on the black hat algorithm, which can effectively reduce the influence of background noise by subtracting the original image before matching filter processing and the obtained image after processing to enhance the characteristics of blood vessels. We performed a series of processing transformations as shown in equations (8) and (9) on the images processed by large-scale matched filtering. We call this algorithm black hat2.(8)Bhatf=fx,y∙bu,v−fx,y,(9)gx,y=255−fx,y−2∗Bhatf,where ∙ is the morphological closing operation and bu,v is disk type structuring element, Bhatf is the black hat transformation, fx,y is the original image, and gx,y is the final processed image.Figure 5
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) large-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 2.4.3. Channel 3

This paper uses a small-scale Gaussian matched filter to enhance the image of small blood vessels, as shown in Figure6. After many experiments, the parameters of the matched filter are set as l=5, s=0.1, and 18 directions which means i=1,2,⋯,18 in equation (3). Using small-scale filters can effectively enhance the small blood vessels in the image, but at the same time, it also enhances much striped noise in the image, and the enhancing effect on the thick blood vessels with central reflection is poor. To reduce the background noise, the black hat2 algorithm used in channel 2 is also used to remove the background in channel 3.Figure 6
(a, f, k) The grayscale images after preprocessing operation, (b, g, l) small-scale matched filtered images, background extracted by (c, h, m) close operation and (d, i, n) subtraction operation, and (e, j, o) final results. The different samples of (a–e) DRIVE, (f–j) STARE, and (k–o) CHASE_DB1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)


## 2.5. U-Net Model

In image semantic segmentation using deep learning, the U-Net network model is the most widely used, which is improved based on the classic full convolutional network (FCN) [24]. U-Net is an image-to-image pixel-level classification network, and its network structure is apparent, as shown in Figure 7. U-Net is different from other standard segmentation networks: U-Net uses an entirely different feature fusion method—splicing. U-Net stitches the features together in the channel dimension. This method fuses the in-depth features extracted from the image with the shallow features to form thicker features, while the fusion operation of FCN only uses corresponding point addition and does not obtain thicker features.Figure 7
U-Net model architecture.Unlike the structure in the original literature [24], this paper sets the padding value of 1 in each layer’s convolution operation, and the convolution kernel size is 3∗3. The purpose is to ensure that the output and input image sizes are consistent and avoid the size increasing operation in the output layer. It is essentially a binary classification operation in the output layer of U-Net. We use an adaptive threshold segmentation algorithm for processing in this paper. The idea of this algorithm is not to calculate the global image threshold but to calculate the local threshold according to different areas of the image, so for different areas of the image, the algorithm can adaptively calculate different thresholds and perform binary segmentation. The specific calculation process is shown in(10)T=−b+12m+1×2n+1∑i=0n∑j=0mgx±i,y±j,where b is the fixed parameter, 2m+1×2n+1 is the area, and T is the area’s threshold.This paper proposes a new loss function that combines the Dice coefficient with the two-class cross-entropy loss function. The Dice coefficient is widely used in the evaluation of image segmentation. In order to facilitate the formation of the minimized loss function, as shown in(11)Ldice=1−2X∩YX+Y,where X∩Y represents the common elements of the prediction graph and the label graph, X and Y represent the number of elements of the prediction graph and the label. In order to facilitate the calculation, approximate ∣X∩Y∣ as the dot product between the predicted probability map and the label, and add the elements in the result. ∣X∣ and ∣Y∣ are quantified by summing the squares of each element. As shown in(12)Ldice=1−2∑iNpk,iqk,i∑iNp2k,i+∑iNq2k,i,where N is the number of pixels, pk,i∈0,1 and qk,i∈0,1 are the predicted probabilities and true labels of the pixel belonging to category k.The cross-entropy loss function used to optimize the network is shown as(13)Lr=−∑iN1−TPNpylogp+1−TNNn1−ylog1−p,where TP and TN are the numbers of true positive and true negative pixels, respectively; Np and Nn are the numbers of segmented pixels and nonsegmented pixels, respectively; y is the label value (y=1, segmentation target; y=0, background); and p is the predicted probability value of the pixel.A coefficientλ is introduced to define the new loss function Loss, as shown in(14)Loss=Ldice+λLr.Notably, the coefficientλ is set to 0.5 in this work, and the flowchart of U-Net is summarized in Algorithm 1.Algorithm 1: Training of U-Net with dynamic learning rate.
Input: Train imagesX, ground truth GInput: Initial epochs E⟵ 30, batch size ⟵ 1, learning rate lr⟵ 0.01Input: Initialize best loss bl⟵ InfOutput: Predicted images P, U-Net parameter1.Xpre⟵ preprocessing (X)2.Xenh⟵ enhancement (Xpre)3.fore⟵ 0 to Edo4.ife=1/3∗Ethen5.lr⟵0.1∗lr6.else ife=2/3∗Ethen7.lr⟵0.1∗lr8.end if9.N⟵ compute the number of train imagesX10.par⟵ initial parameter of U-Net11.whilen⟵1<N+1do12.Yn⟵Unet(Xenhn,par)13.lossn⟵DiceYn,Gn+λCross_entropyYn,Gn14.iflossn<blthen15.bl⟵lossn16.par⟵newpar17.end if18.par⟵SGDpar,lr19.end while20.end for21.P⟵UnetXenh,par22.returnP, par


## 2.6. Postprocessing

In the postprocessing, since the final segmentation image merges the three segmentation images, the noise in the resulting image is also superimposed on all the noises of the three images. Noises will undoubtedly have a significant impact on the actual effect of the segmented image, so this paper addresses this issue in the final postprocessing step. In this paper, a morphological algorithm is used to calculate the size of the connected area of the image. The 8-adjacent connection method is adopted to eliminate the area with the connected area less than 25 pixels, which is to reclassify the area pixels as background. This paper selects a test image in the DRIVE dataset for experimental comparison, and the comparison images are shown in Figure8.Figure 8
Comparison of postprocessing: (a–e) segmentation image without postprocessing; (f–j) segmentation image applying postprocessing.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


## 2.7. Experimental Design

### 2.7.1. U-Net Implementation Details

The U-Net model used in this paper is slightly different from the structure in literature [24]. In order to keep the input and output image sizes of the model consistent, the convolution structure is adjusted accordingly. The specific model structure parameters are shown in Table 1.Table 1
The parameters of the U-Net architecture.
Block nameLayer nameImage sizeParametersDoubleConvConv (ksize=3, pad=1)3∗3∗C1+1∗C2BN+LReLU2∗C2Conv (ksize=3, pad=1)3∗3∗C2+1∗C2BN+LReLU2∗C29∗C1+C2∗C2+6∗C2Input1×576×5760Encoder block_1DoubleConv_164×576×5769∗1+64∗64+6∗64=37824Maxpooling (ksize=2)64×576×5760Encoder block_2DoubleConv_2128×288×2889∗64+128∗128+6∗128=221952Maxpooling (ksize=2)128×288×2880Encoder block_3DoubleConv_3256×144×1449∗128+256∗256+6∗256=886272Maxpooling (ksize=2)256×144×1440Encoder block_4DoubleConv_4512×72×729∗256+512∗512+6∗512=3542016Maxpooling (ksize=2)512×72×720Encoder block_5DoubleConv_5512×36×369∗512+512∗512+6∗512=4721664Maxpooling (ksize=2)512×36×360Decoder block_1Upsampling (bilinear)512×72×720Concat1024×72×720DoubleConv_6256×72×729∗1024+256∗256+6∗256=2950656Decoder block_2Upsampling (bilinear)256×144×1440Concat512×144×1440DoubleConv_7128×144×1449∗512+128∗128+6∗128=738048Decoder block_3Upsampling (bilinear)128×288×2880Concat256×288×2880DoubleConv_864×288×2889∗256+64∗64+6∗64=184704Decoder block_4Upsampling (bilinear)64×576×5760Concat128×576×5760DoubleConv_964×576×5769∗128+64∗64+6∗64=110976OutputConv (ksize=1)1×576×5761∗1∗64+1=65In training, we set the epoch to 30 and the initial learning ratelr to 0.01, and then, the learning rate is set to update in a three-stage formula, as shown in(15)lr=0.01,epoch>10,0.001,10<epoch≤20,0.0001,20<epoch≤30.Setting a larger learning rate at the beginning is to make the model obtain the vicinity of the optimal global parameters faster, and this operation can reduce the training time of the model. After training for a particular epoch, the learning rate needs to be reduced accordingly in order to make the parameters closer to the optimal value in subsequent updates. The stochastic gradient descent (SGD) algorithm is used in the optimization of the loss function.


### 2.7.2. Training Image Preparation

We randomly select 15 images from STARE and the first 20 images from CHASE_DB1 as their respective training set. Due to the limited number of images in the existing dataset, to avoid the overfitting phenomenon in the model training, we perform data expansion processing on the training set of each dataset. Thanks to the translation invariance of the convolutional structure, the images in the training set in this paper were flipped horizontally and vertically and rotated 180 degrees to increase the amount of data 4 times.


### 2.7.3. Measuring Metrics

In order to evaluate the segmentation performance of this algorithm, we use the following metrics to perform a comprehensive evaluation of the segmentation result. These metrics are accuracy (ACC), sensitivity (Se), specificity (Sp), and AUC and calculated as follows:(16)ACC=TP+TNTP+FN+TN+FP,(17)Se=TPTP+FN,(18)Sp=TNTN+FP,(19)AUC=12TPTP+FN+TNTN+FP,where TP is true positive, FP is false positive, TN is true negative, and FN is false negative. Se is the sensitivity, which indicates the degree of classification of blood vessels and nonvascular pixels. In this paper, higher sensitivity indicates that more tiny blood vessels can be detected. Sp is specificity, which is used to express the ability of the algorithm to recognize nonvascular pixels. ACC is the accuracy of algorithm segmentation, reflecting the gap between the algorithm segmentation result and the natural result. AUC is the area under the ROC curve, and we adopt another calculation method to get the AUC, as shown in equation (19) [11].Besides, we also use two other evaluation metrics to measure the effect of segmentation: MCC and CAL.(20)MCC=TP×TN−TP×FNTP+FP×TP+FN×TN+FP×TN+FN.MCC is a correlation coefficient between the segmentation output of the algorithm and ground truth. It comprehensively considers TP, TN, FP, and FN, which is a relatively balanced metric. Finally, it is more suitable for an imbalanced class ratio.CAL can be expressed as the product ofC, A, and L as follows:(21)fC,A,L=C×A×L.SupposeS and SG are the segmentation result and the corresponding ground truth, respectively. These functions are defined as follows:(1)
Connectivity (C): it evaluates the fragmentation degree between S and SG by comparing the number of connected components:(22)C=1−min1,#CSG−#CS#SG,where #C∙ means the number of connected components, while #∙ means the number of vessel pixels in the considered binary image.(2)
Area (A): it evaluates the degree of intersecting area between S and SG and is defined as(23)A=#δαS∩SG∪δαSG∩S#S∪SG,where δα· is a morphological dilation using a disc of α pixels in radius. We set α=2.(3)
Length (L): it evaluates the equivalent degree between S and SG by computing the total length:(24)L=#φS∩δβSG∪δβS∩φSG#φS∪φSG,where φ· is the homotopic skeletonization and δβ∙ is a morphological dilation with a disc of β pixel in radius. We set β=2.According to [26], the CAL metric is essential to quantify thick and thin vessels more equally.


## 2.7.1. U-Net Implementation Details

The U-Net model used in this paper is slightly different from the structure in literature [24]. In order to keep the input and output image sizes of the model consistent, the convolution structure is adjusted accordingly. The specific model structure parameters are shown in Table 1.Table 1
The parameters of the U-Net architecture.
Block nameLayer nameImage sizeParametersDoubleConvConv (ksize=3, pad=1)3∗3∗C1+1∗C2BN+LReLU2∗C2Conv (ksize=3, pad=1)3∗3∗C2+1∗C2BN+LReLU2∗C29∗C1+C2∗C2+6∗C2Input1×576×5760Encoder block_1DoubleConv_164×576×5769∗1+64∗64+6∗64=37824Maxpooling (ksize=2)64×576×5760Encoder block_2DoubleConv_2128×288×2889∗64+128∗128+6∗128=221952Maxpooling (ksize=2)128×288×2880Encoder block_3DoubleConv_3256×144×1449∗128+256∗256+6∗256=886272Maxpooling (ksize=2)256×144×1440Encoder block_4DoubleConv_4512×72×729∗256+512∗512+6∗512=3542016Maxpooling (ksize=2)512×72×720Encoder block_5DoubleConv_5512×36×369∗512+512∗512+6∗512=4721664Maxpooling (ksize=2)512×36×360Decoder block_1Upsampling (bilinear)512×72×720Concat1024×72×720DoubleConv_6256×72×729∗1024+256∗256+6∗256=2950656Decoder block_2Upsampling (bilinear)256×144×1440Concat512×144×1440DoubleConv_7128×144×1449∗512+128∗128+6∗128=738048Decoder block_3Upsampling (bilinear)128×288×2880Concat256×288×2880DoubleConv_864×288×2889∗256+64∗64+6∗64=184704Decoder block_4Upsampling (bilinear)64×576×5760Concat128×576×5760DoubleConv_964×576×5769∗128+64∗64+6∗64=110976OutputConv (ksize=1)1×576×5761∗1∗64+1=65In training, we set the epoch to 30 and the initial learning ratelr to 0.01, and then, the learning rate is set to update in a three-stage formula, as shown in(15)lr=0.01,epoch>10,0.001,10<epoch≤20,0.0001,20<epoch≤30.Setting a larger learning rate at the beginning is to make the model obtain the vicinity of the optimal global parameters faster, and this operation can reduce the training time of the model. After training for a particular epoch, the learning rate needs to be reduced accordingly in order to make the parameters closer to the optimal value in subsequent updates. The stochastic gradient descent (SGD) algorithm is used in the optimization of the loss function.


## 2.7.2. Training Image Preparation

We randomly select 15 images from STARE and the first 20 images from CHASE_DB1 as their respective training set. Due to the limited number of images in the existing dataset, to avoid the overfitting phenomenon in the model training, we perform data expansion processing on the training set of each dataset. Thanks to the translation invariance of the convolutional structure, the images in the training set in this paper were flipped horizontally and vertically and rotated 180 degrees to increase the amount of data 4 times.


## 2.7.3. Measuring Metrics

In order to evaluate the segmentation performance of this algorithm, we use the following metrics to perform a comprehensive evaluation of the segmentation result. These metrics are accuracy (ACC), sensitivity (Se), specificity (Sp), and AUC and calculated as follows:(16)ACC=TP+TNTP+FN+TN+FP,(17)Se=TPTP+FN,(18)Sp=TNTN+FP,(19)AUC=12TPTP+FN+TNTN+FP,where TP is true positive, FP is false positive, TN is true negative, and FN is false negative. Se is the sensitivity, which indicates the degree of classification of blood vessels and nonvascular pixels. In this paper, higher sensitivity indicates that more tiny blood vessels can be detected. Sp is specificity, which is used to express the ability of the algorithm to recognize nonvascular pixels. ACC is the accuracy of algorithm segmentation, reflecting the gap between the algorithm segmentation result and the natural result. AUC is the area under the ROC curve, and we adopt another calculation method to get the AUC, as shown in equation (19) [11].Besides, we also use two other evaluation metrics to measure the effect of segmentation: MCC and CAL.(20)MCC=TP×TN−TP×FNTP+FP×TP+FN×TN+FP×TN+FN.MCC is a correlation coefficient between the segmentation output of the algorithm and ground truth. It comprehensively considers TP, TN, FP, and FN, which is a relatively balanced metric. Finally, it is more suitable for an imbalanced class ratio.CAL can be expressed as the product ofC, A, and L as follows:(21)fC,A,L=C×A×L.SupposeS and SG are the segmentation result and the corresponding ground truth, respectively. These functions are defined as follows:(1)
Connectivity (C): it evaluates the fragmentation degree between S and SG by comparing the number of connected components:(22)C=1−min1,#CSG−#CS#SG,where #C∙ means the number of connected components, while #∙ means the number of vessel pixels in the considered binary image.(2)
Area (A): it evaluates the degree of intersecting area between S and SG and is defined as(23)A=#δαS∩SG∪δαSG∩S#S∪SG,where δα· is a morphological dilation using a disc of α pixels in radius. We set α=2.(3)
Length (L): it evaluates the equivalent degree between S and SG by computing the total length:(24)L=#φS∩δβSG∪δβS∩φSG#φS∪φSG,where φ· is the homotopic skeletonization and δβ∙ is a morphological dilation with a disc of β pixel in radius. We set β=2.According to [26], the CAL metric is essential to quantify thick and thin vessels more equally.


## 3. Results and Discussion

As shown in Figure9, one test image is selected from each of the three datasets to display the segmentation results of each channel and the fusion results. It can be seen that some of the intermittent blood vessels of each channel are reconnected after fusion, and the number of small blood vessels in the fusion map is significantly higher than that of each channel segmentation map.Figure 9
Performance of each channel’s segmentation result.
(a)
Retinal image(b)
Channel 1(c)
Channel 2(d)
Channel 3(e)
Fusion resultThe DRIVE dataset is selected as the experimental object and compares the three channels’ metric data in this paper. The results show that the overall fusion effect of the three channels is better than the segmentation results of every single channel; in particular, the sensitivity has been dramatically improved, as shown in Table2.Table 2
Segmentation results of improvements on DRIVE.
ChannelSeSpACCAUCChannel 10.81740.97680.96260.8971Channel 20.80080.97410.95870.8875Channel 30.81130.97480.96330.89313 channels of fusion0.87450.96240.95460.9185To illustrate this paper’s segmentation effect, we list various metrics on the DRIVE, STARE, and CHASE_DB1 datasets of different papers in recent years in Tables3–5. It can be seen that the algorithm in this paper is superior to most similar papers in sensitivity and AUC metrics. To have a more comprehensive understanding of the overall segmentation effect of the test set, we show the relevant indicators of the prediction results of all test set images in Table 6. The other essential metrics are MCC and CAL, and they achieved by the proposed method has been contrasted with existing segmentation techniques on the DRIVE, STARE, and CHASE_DB1 datasets shown in Table 7.Table 3
Comparison of the proposed method with other methods on the DRIVE dataset.
MethodSeSpACCAUCKhan et al. (2016) [9]0.73730.96700.95010.8522Khan et al. (2016) [29]0.7800.9720.9520.876Soomro et al. (2017) [11]0.7460.9170.9460.8315Ngo and Han (2017) [12]0.74640.98360.95330.8650Biswal et al. (2017) [30]0.710.970.950.84Yan et al. (2018) [13]0.76530.98180.95420.8736Oliveira et al. (2018) [15]0.80390.98040.95760.8922Wang et al. (2019) [25]0.76480.98170.95410.8733Guo et al. (2019) [31]0.78000.98060.95510.8803Feng et al. (2019) [14]0.76250.98090.95280.8717Ribeiro et al. (2019) [32]0.78800.98190.95690.8850Dharmawan et al. (2019) [18]0.83140.9726—0.902Saroj et al. (2020) [5]0.73070.97610.95440.8534Dash and Senapati (2020) [33]0.74030.99050.96610.8654Biswas et al. (2020) [16]0.78230.98140.95610.8819Budak et al. (2020) [34]0.74390.99000.96850.86702nd human observer0.77600.97240.94720.8742Proposed method0.87450.96240.95460.9185Table 4
Comparison of the proposed method with other methods on the STARE dataset.
MethodSeSpACCAUCKhan et al. (2016) [9]0.73590.97080.95020.8534Khan et al. (2016) [35]0.77280.96490.95180.8689Khan et al. (2017) [36]0.7780.9660.9510.872Soomro et al. (2017) [11]0.7480.9220.9480.835Biswal et al. (2017) [30]0.700.970.950.835BahadarKhan et al. (2017) [37]0.7580.9630.9460.861Yan et al. (2018) [13]0.75810.98460.96120.8714Oliveira et al. (2018) [15]0.83150.98580.96940.9087Wang et al. (2019) [25]0.75230.98850.96400.8704Guo et al. (2019) [31]0.82010.98280.96600.9015Feng et al. (2019) [14]0.77090.98480.96330.8779Dharmawan et al. (2019) [18]0.79240.9827—0.8876Saroj et al. (2020) [5]0.72780.97240.95090.8501Tamim et al. (2020) [38]0.78060.98250.96320.88162nd human observer0.89520.93840.93490.9168Proposed method0.89030.97440.96990.9323Table 5
Comparison of the proposed method with other methods on the CHASE_DB1 dataset.
MethodSeSpACCAUCBiswal et al. (2017) [30]0.760.97—0.865Yan et al. (2018) [13]0.76330.98090.96100.8721Oliveira et al. (2018) [15]0.77790.98640.96530.8822Wang et al. (2019) [25]0.77300.97920.96030.8761Guo et al. (2019) [31]0.78880.98010.96270.8845Soomro et al. (2019) [39]0.80200.9680.8910.885Tamim et al. (2020) [38]0.75850.98460.95770.8716Joshua et al. (2020) [40]0.77960.98640.97220.88302nd human observer0.76860.97790.95600.8733Proposed method0.89160.95960.95610.9256Table 6
Segmentation results of all test images of the three datasets.
ImageACCSeSpAUCDRIVE01_test0.9460.9280.9470.93802_test0.9520.9140.9560.93503_test0.9550.8170.9700.89404_test0.9590.8680.9680.91805_test0.9580.8380.9710.90406_test0.9580.8110.9730.89207_test0.9540.8510.9640.90708_test0.9580.8200.9710.89609_test0.9590.8490.9690.90910_test0.9570.8630.9650.91411_test0.9450.8700.9520.91112_test0.9580.8750.9660.92013_test0.9530.8590.9630.91114_test0.9540.9010.9590.93015_test0.9510.9170.9540.93516_test0.9540.8890.9610.92517_test0.9580.8450.9680.90718_test0.9540.9130.9580.93519_test0.9540.9370.9560.94620_test0.9550.9250.9570.941Avg.0.9550.8750.9620.918STAREim00020.9720.8390.9810.910im00770.9670.9660.9610.964im01630.9610.9760.9600.968im02550.9700.8720.9790.926im02910.9800.7980.9900.894Avg.0.9700.8900.9740.932CHASE_DB111L0.9460.9370.9470.94211R0.9420.9500.9420.94612L0.9530.8780.9590.91912R0.9580.8720.9650.91813L0.9580.8840.9630.92313R0.9560.8500.9630.90714L0.9700.8950.9680.93114R0.9660.8670.9710.919Avg.0.9560.8920.9600.926Table 7
MCC and CAL metrics of existing techniques on the three datasets.
MethodDRIVESTARECHASE_DB1MCCCALMCCCALMCCCALAzzopardi et al. (2015) [41]0.7190.7210.6980.7090.6560.608Orlando et al. (2016) [42]0.7400.6750.7260.6650.6890.571Dharmawan et al. (2017) [18]079910.88340.79590.8181——Yang et al. (2018) [43]0.725—0.662———Strisciuglio et al. (2019) [44]0.7290.7280.6980.7090.6630.620Khan et al. (2020) [45]0.7390.6960.7070.5660.6290.5472nd human observer0.7700.7710.7410.6220.6260.722Proposed method0.7560.7960.7960.8370.5660.733We selected image 19_test from the test set of the DRIVE dataset to display the segmentation results, as shown in Figure10. Literature [5, 27] segmented some small blood vessels, but it is still slightly insufficient compared to this paper’s segmentation diagram. Literature [10] lacks many details, and the small blood vessels are not segmented. The segmentation result of literature [11] contains a lot of edge noise, and there are many intermittent blood vessels. Compared with the existing segmentation methods, the segmentation results in this paper have a good performance in terms of the integrity of the whole blood vessels and the segmentation of small blood vessels.Figure 10
Comparison of different methods on the DRIVE dataset.
(a)
Original image(b)
Ground truth(c)
Literature [5](d)
Literature [11](e)
Literature [14](f)
Literature [10](g)
Literature [27](h)
Proposed resultAs shown in Figure11, we select the test results of the image im0163 in the STARE dataset for comparison. It can be shown that the segmentation results of this paper are similar to those of literature [13, 14], but the background noise in literature [13] is not eliminated. Compared with literature [5, 10, 27], the algorithm in this paper illuminates the optic disc structure in the original image as much as possible in the preprocessing part, so the problem that is incorrectly dividing part of the optic disc structure into blood vessels like these papers did not appear in the final segmentation result.Figure 11
Comparison of different methods on the STARE dataset.
(a)
Original image(b)
Ground truth(c)
Literature [5](d)
Literature [13](e)
Literature [14](f)
Literature [10](g)
Literature [27](h)
Proposed resultThe CHASE_DB1 dataset is not used in most of the papers about retinal blood vessel segmentation. One of the reasons is that the dataset contains half of the abnormal images, which may cause some interference to the trained segmentation model. Meanwhile, this dataset is also a new and challenging dataset compared to the classic DRIVE and STARE datasets. We selected four images image_12R, image_13L, image_13R, and image_14L from the test set of the CHASE_DB1 dataset to compare the segmentation results in order to verify the generalizability of the proposed algorithm, as shown in Figure12. The segmentation result of the algorithm in literature [19] has much noise, and some blood vessels are not effectively separated. Literature [28] does an excellent job in the segmentation of small blood vessels, but there is a problem that some blood vessels are not connected. Due to the postprocessing in this paper, the segmentation result on this dataset contains less noise and guarantees the continuity of most blood vessels. However, compared with the manual label, some tiny blood vessels cannot be completely segmented from the image background.Figure 12
(a, e, i, m, q) Image_12R, (b, f, j, n, r) image_13L, (c, g, k, o, s) image_13R, and (d, h, l, p, t) image_14L from the CHASE_DB1 dataset. (a–d) Original images, (e–h) ground truth, (i–l) literature [19], (m–p) literature [28], and (q–t) proposed segmentation images.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)(q)(r)(s)(t)The source codes of the proposed framework have been running on the PC (Intel Core i5-6300HQ CPU, 2.30 GHz, 12.0 GB RAM, NVIDIA GTX 950M GPU). DRIVE, STARE, and CHASE_DB1 have spent 11.3 h, 7.1 h, and 16.4 h on training separately in each channel. The average testing time of test images was 1.34 s. Table8 shows the parameter comparison of the proposed method with other methods based on U-Net, which can help us compare the framework complexity of different methods. Note that the parameters are not equal to the training time because some methods use slices of a train image as input of the network. For example, literature [19] has 42421 slices as the training set, which means it needs more time to train the network.Table 8
Network comparison of the proposed method with other methods based on U-Net.
MethodInput sizeEpochTraining imagesParametersRonneberger et al. (2015) [24]572∗572∗1N/AN/A28.94 MYan et al. (2018) [13]128∗128∗1Over 3026052 (4)30.96 MJiang et al. (2018) [19]500∗500∗33042421 (4)58.31 MSoomro et al. (2019) [39]Original imageN/A3959 (4)4.71 MJoshua et al. (2020) [40]512∗512∗350113 (4)0.64 MZhang et al. (2020) [17]256∗256∗1120116 (3)3.86 MBudak et al. (2020) [34]48∗48∗3Over 30Over 45840 (2)0.97 MProposed method576∗576∗130220 (3)13.39 Mxy means x training images of y datasets.


## 4. Conclusion

This paper proposes a new retinal blood vessel segmentation method, which combines a multiscale matched filter with a U-Net neural network model of deep learning. First of all, we use an improved morphological image algorithm to effectively reduce the impact of image background in feature extraction. Additionally, in order to avoid ignoring the characteristics of small blood vessels, this paper performs multichannel feature extraction and segmentation on retinal blood vessel images. Finally, the segmented images of the three channels are merged, and various characteristics of retinal blood vessels are obtained as much as possible. In the training of the U-Net model, we used the loss function weighted by the Dice coefficient and the binary cross-entropy to solve the image pixel imbalance problem. The algorithm of this paper is tested on the existing public datasets DRIVE, START, and CHASE_DB1. The experimental results show that there is better performance in four metrics compared with similar papers. The average sensitivity of the algorithm in this paper reached 0.8745, 0.8903, and 0.8916 on the DRIVE, STARE, and CHASE_DB1 datasets, respectively. This result is nearly 0.1 higher than the average sensitivity of other papers. The improvement of the sensitivity metric also reflects that the algorithm in this paper has a good performance in extracting small blood vessels. The focus of this paper is to combine the advantages of unsupervised algorithms and supervised algorithms. We did not make too many improvements to the U-Net network. Therefore, how to prune the deep learning network model structure will be an interesting research direction in the future.



---

*Source: 5561125-2021-05-26.xml*