A new approach for sensitivity improvement of retinal blood vessel segmentation in high-resolution fundus images based on phase stretch transform

,

The eye-fundus photograph is widely used for eye examinations. Accurate identification of retinal blood vessels could reveal information that is helpful for clinical diagnoses of many health disorders. Although many researchers have proposed different methods to segment images of retinal blood vessels, the sensitivity is plausible to be improved. The algorithm's sensitivity refers to the algorithm's ability to identify retinal vessel pixels correctly. Furthermore, the resolution and quality of retinal images are improving rapidly. Consequently, new segmentation methods are in demand to overcome issues from high-resolution images. This study presented improved performance of retinal vessel segmentation using a novel edge detection scheme based on the phase stretch transform (PST) function as its kernel. Before applying the edge detection stage, the input retinal images were preprocessed. During the preprocessing step, non-local means filtering on the green channel image, followed by contrast limited adaptive histogram equalization (CLAHE) and median filtering, were applied to enhance the retinal image. After applying the edge detection stage, the post-processing steps (CLAHE, median filtering, thresholding, and morphological operations) were implemented to obtain the segmented image. The proposed method was evaluated using images from the highresolution fundus (HRF) public database and yielded promising results for sensitivity improvement of retinal blood vessel detection. The proposed approach contributes to a better segmentation performance with an average sensitivity of 0.813, representing a clear improvement compared to several benchmark techniques. processing methods can be applied to overcome problems in manual screening. An automated system that can automatically detect anomalies in the retinal blood vessels, such as diabetic retinopathy markers, will reduce those influential factors while improving screening efficiency [6], [7].
In recent years, retinal blood vessel detection in the retinal fundus image has been extensively studied [8]. Most researchers took dataset samples from two well-known benchmark databases: Structured Analysis of the Retina (STARE) and Digital Retinal Images for Vessel Extraction (DRIVE), to validate their methods [9], [10]. However, these databases provide lower-resolution retinal images than those collected by the latest fundus cameras. Datasets in both databases were dated from the early 2000s when fundus cameras could not produce high-resolution images. Recently, the resolution and quality of retinal fundus images have been rapidly improving. Typical image processing techniques that were developed and evaluated using low-resolution images have demonstrated limitations in clinical applications. Many segmentation algorithms tested using low-resolution retinal images are unsuitable for most applications in clinical cases with retinal images acquired from the latest fundus camera that produces high-resolution images. Thus, new segmentation methods are in demand to overcome issues from high-resolution images. In recent findings, databases containing higher-resolution fundus images, such as the High-Resolution Fundus (HRF) database, have been more prominent [11]. The following highlight several important studies that implement the HRF database, and algorithms in these studies will later be used for performance comparison with the approach proposed in this study.
Budai et al. [11] proposed a Frangi-based vessel enhancement algorithm combined with a hysteresis threshold to increase sensitivity and a multiresolution framework to reduce the computational load. The Frangi method was a mathematical model-based technique that extracted vesselness features based on the eigenvalues of the Hessian matrix. Meanwhile, Odstrcilik et al. [12] employed a technique based on the least error thresholding in matched filtering. Contrast equalization and illumination correction were used in the preprocessing step. To address challenges in pathological images with exudates and other abnormalities, Annunziata et al. [13] employed unsupervised vessel segmentation. During vessel enhancement, an inpainting filter was proposed to paint exudates so that nearby false positives were considerably decreased.
Researchers also employed approaches using conditional random field (CRF) for retinal vascular segmentation. For example, a discriminatively trained extraction model using a fully connected CRF for retinal vascular extraction was presented by Orlando et al. [14]. Zhou et al. [15] employed a dense CRFbased method. The method was improved by implementing thin-vessel enhancement and discriminative feature learning. Likewise, Dharmawan et al. [16] modified the Dolph-Chebyshev function to develop a direction-sensitive matched filter. The algorithm utilized a technique for combining responses of the matched filter bank. In the preprocessing phase, background normalization was applied using an inpainting technique to eliminate the effect of exudates and increase the contrast between the vessel and the background.
Sundaram et al. [17] proposed a hybrid extraction method with a mask generation scheme to segment vessels from retinal images. The proposed method consists of three phases. Phase one is preprocessing to enhance the image using contrast stretching. In phase two, a mask for optic papilla removal was applied. Then in the last phase, the image was binarized using adaptive thresholding, and a multi-scale vessel enhancement was also applied. Meanwhile, Wang and Li [18] employed a level set technique, using the ChanVese region-based term with a distance regularisation and a Gaussian mixture term. Preprocessing steps removed the optical disc's influence and enhanced the blood vessels. At this stage, the morphological closing operation and matched filtering were applied. The post-processing stage removed the background noise using a length filter. Alternatively, the Bar-Combination of Shifted Filter Responses (B-COSFIRE) filter was used by Ali et al. [19] for retinal vessel segmentation. Before applying the segmentation algorithm, Ali's technique added padding to the processed image.
Concisely, many researchers have applied different techniques to segment blood vessels in highresolution retinal fundus images. However, an examination of the sensitivity parameter reveals that these techniques are plausible for improvement. Sensitivity is defined as the probability of identifying vessel pixels, indicating an algorithm's effectivity. Moreover, Chalakkal et al. [20] reviewed studies on screening and diagnosing eye diseases using retinal image analysis and concluded that it was essential to improve sensitivity.
Along with advances in digital image processing, Asghari and Jalali [21] proposed a digital image enhancement and processing method known as phase stretch transform (PST). It was first introduced as an approach to detect edges and sharp transitions in digital images and has been publicly available on Matlab Central File Exchange and GitHub. The algorithm applied a nonlinear dispersive phase operation to the image sample. The algorithm also included image preprocessing and post-processing in its steps. For noise reduction, the Gaussian low pass filter was applied in the preprocessing step. The image postprocessing included thresholding, morphological thinning, and removing isolated pixels for a clean result. The transform's output phase showed changes in image intensity and could be used for feature extraction and edge recognition. PST has been applied as a texture and edge extractor and performs exceptionally well [22], [23]. Hence, this paper used PST as an approach to segment blood vessels from high-resolution retinal fundus images and evaluated its performance by comparing the results with other segmentation techniques.
Retinal fundus images are obtained using a fundus camera, captured at the back of the eye. However, these images could have low illumination and varying contrast, resulting in low-quality images. Improving the quality of retinal images is crucial in acquiring well-segmented images [24]. Considering this issue, we present a novel scheme of preprocessing and post-processing steps in the study to improve the performance of the segmentation. The main contribution of the proposed method is a new improvement scheme of image preprocessing and post-processing of the PST edge detection method in segmenting the blood vessels from high-resolution retinal fundus images. In this study, the proposed preprocessing steps began with non-local means filtering on the green channel for noise reduction. Next, a contrast enhancement was performed, followed by median filtering to reduce the salt and pepper noise. Contrast enhancement and median filtering were also applied in the post-processing steps, followed by Otsu's thresholding to convert the image into a binary image. The last step in image post-processing was the closing morphology operation to omit object holes. In brief, this study proposes a new scheme of preprocessing and post-processing steps in retinal vessel segmentation using the PST method to improve sensitivity.
The remainder of this paper will be written in structure as follows. The second section provides a brief overview of the proposed approach. Implementation results, performance analysis, and discussion are demonstrated in the third section, whereas the conclusion is provided in the fourth section.

Dataset
The study collected 45 images divided into three subsets; 15 fundus images of glaucomatous patients, 15 fundus images of patients with diabetic retinopathy, and 15 images of healthy patients from the HRF Image Public Database as object samples. Fundus image samples were RGB-formatted color images with a resolution of 3504 by 2336 pixels. Ground truth images of retinal vessels marked by experts were available in the database and were used in the research for performance evaluation [11]. A scheme of methodological steps for retinal blood vessel extraction conducted in this study is illustrated in Fig. 1.
A pixel in any RGB color image has three elements: red, green, and blue. However, image segmentation only requires one channel value for each pixel. As a result, in the preprocessing step, a channel with the best representation of vessel images must be selected; this process is called channel extraction.
Next, non-local means filtering was applied to the channel, followed by contrast-limited adaptive histogram equalization (CLAHE) for contrast enhancement and noise reduction using a median filter. After image preprocessing, edge detection using PST was applied for vessel segmentation. Background noise, undesired segments, and incorrectly recognized vessel pixels were removed in the post-processing step. This step also included CLAHE, median filtering, Otsu's thresholding, and morphological operations.  In the morphological operation step, two different procedures were applied. The first procedure was used to omit objects with a size smaller than a given threshold, and the second procedure was used to remove the object holes (morphological closing). Then, the methods used in this research will be explained comprehensively as follows.

Pre-Processing 2.2.1. Non-local means filtering
A non-local (NL)-means filtering was applied to the extracted channel to reduce noise. The main principle of non-local means filtering is to replace a pixel value in an image with an average weight of a selection of other neighborhood pixel values [25], [26]. The term [ ]( ) NL v i has been used to refer to a mathematical expression of the averaged weights of entire pixels, revealed in (1).
where v denotes a discrete noisy image, while w indicates the weight function, which values the resemblance of pixels i to j . These weights can be expressed in a mathematical formulation as in (2).
where ( ) Z i is a constant of the normalizing factor, j N acts as a fixed-sized square neighborhood centered on pixel j . Meanwhile, h is the filtering parameter. The h parameter acts as the decay factor of the exponential function and controls the smoothness of filtering results (degree of smoothing or DoS) [27], [28].

Contrast-Limited Adaptive Histogram Equalization (CLAHE)
In retinal images, non-vessel and vessel segments are challenging to be distinguished due to the small variation gap in image intensity. An intensity transformation widens the gap between non-vessel and vessel pixels, enhancing the image contrast. One way to enhance the contrast is to use the Histogram Equalization (HE) method. CLAHE is a well-known contrast enhancement technique in biomedical image processing due to its high effectiveness in transforming ordinary parts into more visible leapingout parts. This method splits the image into separated regions and applies local histogram equalization in each area, followed by bilinear interpolations to eliminate the boundaries between the regions [29].

Median Filtering
Despite improving the image, the contrast enhancement method also increased the noise in the image. Then, a smoothing technique was introduced in this paper to lessen the noise resulting from the enhancement method. A smoothing method using median filtering was implemented to suppress the blurry edges within the image. The objective of this step was to change the intensity of a pixel in an image with the neighborhood's median intensity. In median filtering, noise can be minimized while preserving the edges. Median filtering has also been used to reduce salt-and-pepper noise [30]

Phase Stretch Transform (PST)
PST highlights edge spatial information in an image by diffracting a 2-dimensional phase kernel function in the frequency domain. In terms of frequency properties, edges in an image have higher frequencies. Meanwhile, the amount of the phase applied to the image by PST is related to the frequency properties; parts of the image with higher frequency will receive a higher amount of phase. In other words, PST detects edges in an image by emphasizing these frequency-dependent properties. This operation is specified in a mathematical expression as in (3) The spatial-frequency operation in PST requires an angle operator ∡〈•〉 and p and q in the operation refer to 2-dimensional frequency variables. Moreover, it performs two different transformations: the two-dimensional Fast Fourier Transform ( 2 FFT ) and the Inverse Fast Fourier Transform ( 2 IFFT ). The two-dimensional FFT is used to transform the image into the frequency domain and perform the phase kernel, while the two-dimensional Inverse FFT retransforms it to the spatial domain. Then, the warped phase kernel, or the nonlinear frequency-dependent phase, which is denoted as [ , ] K p q , is shown in (4).
[ . ] The phase derivative [ , ] PD p q function for edge detection is either linear or sublinear. Such phase derivative profiles can be formed by operating an inverse tangent function, resulting in the following PST kernel phase, as shown in (5).
. . denote the warp and strength of the phase profile applied to the image variables. The phase kernel will be built based on the selection of these variables. Moreover, these variables influence the produced phase kernel and its derivative profile, regulating the magnitude of the phase applied to each frequency [21].

Post-processing 2.4.1. Otsu's thresholding
Otsu's thresholding was used to convert grayscale images to binary images. Otsu presented a nonparametric and unsupervised automatic threshold selection for image segmentation. This method determines the threshold that maximizes the distinctness of the two populations into which the image is divided. This distinctiveness is represented by the interclass variance [31] [32].

Morphological operations
Mathematical morphology in digital image processing is an appropriate tool to analyze shapes in images, which involves two basic processes: erosion and dilation [33]. Dilation is employed to increase the size of the object segment by adding layers around the object. Erosion is the opposite of dilation, in which the object size is reduced by scraping around the object [34]. Morphological operation algorithms can be developed and differentiated based on the sequence of these two processes. An opening algorithm is when erosion is performed first, and dilation is performed subsequently [35]. When the sequence is otherwise, the algorithm is called a closing algorithm. Applying an opening algorithm to an object generally smoothens its contour, breaks narrow isthmuses, and removes thin protrusions. In comparison, applying a closing algorithm tends to smooth out portions of contours in an object. An applied closing algorithm eliminates gaps, fills in tiny holes, combines little breaks, and fuses thin gulfs in the contour

Performance Evaluation
Performance evaluation was carried out by comparing ground truth with the vessel segmentation's results and calculating the sensitivity. It represents the capability ratio of the method to detect the vessel pixels [36]. The sensitivity is defined as (6) [37].
True positive (TP) represents the sum of vessel pixels detected correctly as targetted vessel pixels. Meanwhile, a false negative (FN) denotes the sum of vessel pixels detected wrongly as non-vessel pixels [38].

Results and Discussion
The proposed technique was performed on a PC (Intel Core i5, 2.2 GHz, 12 GB RAM) and was simulated using MATLAB R2018b. The average running time of the proposed technique on an HRF image was 20s. Channel extraction was performed on the original image sample to acquire its channels in each color element, as presented in Fig. 2. Based on the figure, it can be seen that the green channel provided the highest contrast between the background and the vessels [38]. At the same time, the result was not oversaturated or under-illuminated, unlike the other two channels. Hence, the green channel was selected for the following methodological step, the non-local means filtering. In this experiment, the resulting degree of smoothing (DoS) for the non-local means filter was 12, which gave the best retinal blood vessel segmentation result [39]. After image enhancement using the CLAHE method and image smoothing using a median filter, edge detection using PST was applied to the image for vessel extraction. PST, as an operator, can be used to get features in an image. The input of PST was image intensity, while its output was a binary image. Detecting sharp transitions would generate 1 as output, while 0 would be generated elsewhere. The parameters required for the edge detection method are the warp (W) and strength (S) of the phase profile. A combination of these real-numbered parameter values influenced the structure of the phase kernel, along with its derivative profile. Therefore, they also directly affected the phase applied to each frequency. Fig. 3 compares six phase kernels composed of various W and S values. Meanwhile, their generated phase derivative profiles are presented in Fig. 4. Based on Fig. 3 and Fig. 4, the influence of S and W parameters on the structure of the phase kernel and derivative profile can be observed. derivative profile. This can be seen in the figures that different W values with the same S values resulted in almost identical phase profiles. However, the S parameter significantly influenced the structure and amplitude of the phase along with its phase derivative profile. Changing S parameter values would result in a similar structure but a scaled shape; selecting an S value closer to zero would result in a smaller drop-curvature shape. Another qualitative observation that can be noted from the structure of the phase derivative profile is the distinguished pattern of narrowed shape on the p-axis but widened shape on the q-axis. Then, the resulting image was converted from the frequency domain to the spatial domain using 2D IFFT to generate the output phase image. The resulting output phase images would then be compared with ground truth images to calculate the sensitivity of the proposed method. Table 1 highlights the effect of the W and S parameters on the performance of the proposed scheme, which is the sensitivity. Increasing the W parameter did not significantly influence the sensitivity, although there was a slight increase in sensitivity along with the increased value of W.  Moreover, a greater value of W would not only increase the sensitivity but also increase the noise. This precedent follows what was stated by Ang et al. [40] that a greater W produced a sharper edge of the phase derivative but was more noise-prone. However, increasing the S value tended to increase the sensitivity. The highest sensitivity was obtained from W=150 and S=20; these parameter values would later be used for segmentation throughout the study.
From thorough observations of the segmented image, non-vessel pixels were detected in the segmented vessel objects. As shown in Fig. 5(a), white objects had black holes (red circles). A morphological closing procedure, which was a sequence of dilation and erosion, was then employed to improve the image. This procedure was used to close the holes in the object segment. The result showed that black holes could be overcome by converting them into white pixels according to the vessel object's color, as shown in Fig. 5(b). Consequently, the true positive (TP) value increases while the false negative (FN) value decreases, resulting in a higher sensitivity score. The existence of an optic disc can be an obstacle to blood vessel segmentation, so the optic disc needs to be removed [41], [42]. Fig. 6 demonstrates how the PST-based approach could identify blood vessels and remove the optical disc feature without any particular feature removal procedure. Several experiments were conducted empirically on various threshold and disk-shaped structuring elements (SE) values to assess and obtain the most desirable morphological operations. The threshold for the first morphological operation was 4000 to keep objects with a size of ≥4000 pixels. Meanwhile, SE with a radius = 5 was appropriate for the closing algorithm. The visual results of each step of the proposed retinal vessel segmentation method are depicted in Fig. 7.  The sensitivity value was used to measure the performance of the methods. The proposed method was tested using 45 images from the HRF database, and the average sensitivity value acquired was 0.813. In comparison, the initial sensitivity of the technique prior to improvement was 0.432. This indicates a significant improvement in sensitivity. The proposed technique was then compared with other existing vessel segmentation techniques that also employed HRF public databases in their evaluation. Performance comparison was conducted to prove its effectiveness by comparing the sensitivity. The results demonstrated that the proposed new scheme, including filtering, contrast enhancement, and morphological operations, could overcome the problems encountered in the segmentation process related to image contrast and noises. Table 2 summarizes the results obtained from recent study findings.

Conclusion
An important step in developing automated retinal blood vessel analyzers is conducting accurate segmentation of retinal vessels. With accurate blood vessel segmentation, a more detailed observation of vessels can be conducted, such as calculating the vessel width, examining vessel tortuosity, classifying veins and arteries, and clinically diagnosing a disease. The present work's main contribution was to improve the performance of segmenting blood vessels in high-resolution retinal fundus images. The proposed segmentation method outperformed existing techniques for segmenting blood vessels in retinal images from the HRF database, evidenced by a higher sensitivity value. Although this research only suggested a technique for retinal vessel segmentation, this approach provides knowledge support in a more considerable automated retinal image analysis. For the time being, only a single fundus image database, the HRF database, was used to validate the proposed method. However, future works should consider more databases for evaluation to ensure that the proposed method is robust over multiple databases