Automated Segmentation of Retinal Blood Vessels using Fast Marching Method and local mathematical analysis

In this work, Fast Marching Method (FMM) has been suggested for Retinal blood vessels segmentation. FMM is an optimization technique; and the main advantage of the FMM is its ability to deal with branches and bifurcations without any additional computational cost. This advantage had been used in robotics to find the optimal path for the robot to move from the starting point to the goal with no collisions. Considereing the tree structure of blood vessel, I will use FMM to find the shortest bath between the optic disk and the blood vessels ends to draw the tree of the blood vessels.This method has been implemented using the M language in MATLAB R2016b. In this work local mathematical analysis has been implemented so that we can have an initial estimation of blood vessels distribution in an image in order to minimize the huge amount of noise included in retinal images and to make FMM implementing easier. FMM performance had been compared to other techniques used for retinal blood vessel detection like “Matched Filters”. The results showed that the FMM performance overcame some of those techniques and close to other high resolution methods. The FMM algorithm has been validated using the well-known “DRIVE” database and the resulting resolution ranged between 80% to 93% (depending on the noise amount in image) with iteration number between 500 to 1000 (according to the optic disk position in the image) with an average time of 0.57 seconds for each iteration which mean that the total running time is 5-10 minutes. FMM had also been validated using STARE data set and achieved a TPR of 90% for 700x605 STARE images in 15 minutes, and a TPR of 86% in 2.6 minutes when reducing image size to 350x303.


INTRODUCTION
The retina is a light-sensitive tissue that forms the inner layer of the eye. it contains neural tissues and blood vessels and considered as the only tissue of human body we can use to see blood vessel non-invasively which have great importance in diagnosis and treating many diseases such as diabetes, hypertension, cardiovascular diseases and cancers.segmenting retinal blood vessels is a challenging problem because of the variety of illumination conditions, low contrast, small vessel detection, and the trade between accuracy and computational efficiency. Retina can be visualized as an image using the fundus camera. Those images are often noisy, poorly contrasted and non-uniformly illuminated in addition they suffer from brightness variations along the same image and likewise between different images.
In order to segment retinal blood vessel tree, many methods had been reported including supervised methods [1] and unsupervised methods [2]. Many applications had been presented in the diagnosis of various eye and systemic diseases, such as diabetes, hypertension, and angiogenesis.
Typical blood vessel detection methods apply an initial estimation of the vessel network, after that, tuning methods are used in order to improve the accuracy of the vessel detection [3]. Finally, post processing methods are applied to remove false detected vessels.

Al-Rafidain Engineering Journal (AREJ)
Vol. 24, No 2, December 2019, pp. 1-9 In [3] a new multi-scale line-tracking procedure have been proposed. it starts from a small group of pixels, derived from a brightness selection rule, and terminates when a cross-sectional profile condition becomes invalid. The multi-scale image map is derived after combining the individual image maps along scales. This map contains pixels that certainly belong to a vessel. The initial vessel network is derived after map quantization of the multi-scale confidence matrix. After that median filtering is applied in the initial vessel network, and this restores disconnected vessel lines and eliminates noisy lines. At the end, erroneous areas are removed (during postprocessing) using directional attributes of vessels and morphological reconstruction.Author in [3] declared "The experimental evaluation in the publicly available DRIVE database shows accurate extraction of vessels network. The average accuracy of 0.929 with 0.747 sensitivity and 0.955 specificity is very close to the manual segmentation rates obtained by the second observer". In [4] a spatially adaptive contrastenhancement technique wasintegrated with Tyler Coye algorithm, which hasbeen improved with Hough line transformation based vessel reconstruction method. The proposed approach was evaluated onSTARE and DRIVE data sets and comparedwith fivewidely used contrast enhancement techniques based on wavelettransform, contrast limited histogram equalization, localnormalization, linear un-sharp masking and contourlettransform. The results revealed that the enhancementtechnique is well suited for the application and also outperforms all compared techniques. In [5], a set of automated methods were proposed in order to analyze the retinal vessel network and to quantify its morphologic properties taking into consideration both arteries and veins, in twodimensional color fundus images. The analytical methods included Formation of a well-connected vessel network, Structural mapping, Arteryvenous classification and Blood vessel hemorrhage detection. On the other hand, quantification methods included vessel morphology analysis based on the measurement of tortuosity, width, branching angle, branching coefficient, and fractal dimension. Morphological parameters were measured with respect to arteries and veins separately in a vessel network. Then methods have been validated with the manually annotated retinal fundus images as a ground truth.The accuracy of the method was quantified using two metrics 1-misclassification rate: this metric gave a classification error rate of (17.07%) for single vessel trees (without AV crossing) and 4.96% for paired vessel trees (with AV crozzssing). 2-the histogram of pixel misclassification per image: The average misclassification rate was 8.56% which means that the accuracy of correct classification is 91.44%. [5] declared that "The average running time per image starting at the readily available vessel segmentation to AV classification was 8 minutes including 7 minutes for structural mapping and 1 minute for subsequent AV classification, when processed in MatLab environment on a standard personal computer with Intel core 2 Duo processor, running at 3 GHz". In [6] A supervised method was presented to tackle the problem of retinal blood vessel segmentation, This method combines two superior classifiers: Convolutional Neural Network (CNN) and Random Forest (RF). Here, CNN performs as a trainable hierarchical feature extractor and ensemble RFs work as a trainable classifier. By integrating the merits of feature learning and traditional classifier, the method can automatically learn features from the raw images and predict the patterns. experiments have been conducted on DRIVE and STARE databases, and showed that this method overcomesotherother major studies on the same database. In [7], an automatic unsupervised blood vessel segmentation method for retinal images was proposed.
This method constructed a multidimensional feature vector using the green channel intensity and the vessel enhanced intensity feature by the morphological operation. then, it exploited a self-organizing map (SOM) for pixel clustering, which is an unsupervised neural network. Finally, it classified each neuron in the output layer of SOM as retinal neuron or non-vessel neuron with Otsu's method, and got the final segmentation result. In [8], we propose a new infinite active contour model was proposed, it uses hybrid region information of the image to segment vessels. An infinite perimeter regularizer, provided by using L2 Lebesgue measure of the γ -neighborhood of boundaries, allows for better detection of small oscillatory (branching) structures than the traditional models based on the length of a feature's boundaries (i.e., H1 Hausdorff measure). The local phase based enhancement map is used for its superiority in preserving vessel edges while the given image intensity information will guarantee a correct feature's segmentation. This model wasapplyed to three public retinal image datasets (two datasets of color fundus Al-Rafidain Engineering Journal (AREJ) Vol. 24, No 2, December 2019, pp. 1-9 photography and one fluorescein angiography dataset) and it outperforms its competitors when compared with other widely used unsupervised and supervised methods. It achieved an accuracy of (0.954) on the DRIVE dataset .

THE THEORETICALBACKGROUND
Fast marching method is an optimization process uses iterative exploration and exploitation procedures in order to find the minimum path between the starting points and the goal [9]. Assuming a closed curve (Γ) in the plane propagating normal to its self with speed (F). when (F>0) then the front will always move "outwards". The position of this front can be shown by computing the arrival time T(x,y) of this front whenit crosses each point (x,y) as shown in figure (1) [10]. In multi dimensions, the special derivative of the solution surface T becomes the gradient, which means: As a result, we can describe the front motion as the solution to a boundary value problem. And if the speed (F) depends only on position, then the equation can be reduced to the familiar Eikonal equation.
The Eikonal equation can be written as [11]: Finding the solution of this equation is a difficult task and it,in most cases, finding a global smooth solution of this equation is impossible. Thus, many methods have been explored in order to find the optimal solution of this equation. Among those methods we can find the fast marching method (FMM). FMM uses First order numerical approximation to solve the Eikonal equation based on the following operators (assuming that a function u is given with the value u i,j = u(x i,j ) on a Cartesian grid with grid spacing (h)) [9]. According to Godunov [12], The following upwind scheme is used to estimate the gradient in 2D: Al-Rafidain Engineering Journal (AREJ) Vol. 24, No 2, December 2019, pp. 1-9 … (7) Where: T i,j : is the cost function and it gives the (maximum) cost of moving from one point to its neighboring point.
FMM solves this equation using the following algorithm [9]:  At the initial stage: The starting and the goal points are defined and divided into 3 groups: -Accepted points: this group consists of the starting points (in the initial stage). And for later stages it will contain every point that we had computed its distance function (u) and considered it as fixed. -Considered points: this group contains every point next to the points defined as "accepted". At these points the estimation (v) of the distance function (u) is calculated (this value won't be fixed).
-Far points: includes the rest of points in the image. The estimated distance values (v) of these points are not yet computed. In addition, the distance function (u) would be defined as follows: Where (T) is the cost function from one point to the neighboring points.  The iterative stage: this is done as follows -As long as the goal points don't belong to the accepted group, we keep looking for the considered points that have minimum estimation value of the distance function.
-Points of minimum estimation values are called "trial", and once we get them we assign them as "accepted" and their distance function would be equals to their estimation values. i.e.
-After that we update all the points next to the "Trial" ones so that they become "Considered" and have an estimated distance value of:

Back tracking stage:
In this stage the starting points (x start ) and the goal points (x goal ) are defined in addition to the forward path that we had in the exploration procedure. The goal points are assigned to be part of the back tracking path and then an iteration loop starts. During this loop the neighboring points of the goal points are investigated and the points with the least distance function value are chosen to become part of the back tracking path.

THE PROPOSED ALGORITHM
The goal of the FMM is finding the shortest path between two points without crossing constraints existing in between.
The main advantage of the FMM is its ability to deal with branches and bifurcations without any additional computational cost. This advantage had been used in robotics to find the optimal path for the robot to move from the starting point to the goal with no collisions. Applying FMM on mazes is a good way to show the efficiency of this algorithm and its ability to achieve the goal it was designed for since it can be tested on a complex set of obstacles positioned between the starting and the goal points. After that, the FMM became widely used in image processing especially for medical images. And it is obvious that FMM would be completely efficient for tracking blood vessels in an image because of its ability to deal with bifurcations and thus it would be perfect for retinal blood vessel tree segmentation.
This algorithm was applied completely using the M language of MATLAB to extract the lines of blood vessels. We can say that FMM consists of the following main stages: Initial stage: Illumination properties of retinal images are used to detect the optic disk (OD) position automatically. After that, the pixels of the OD are scanned to get the least intensity pixels, those pixels become the starting points of the FMM. After that local analysis of the retinal pixels is implemented in order to select the local least intensity pixels of the image. Those pixels become the expected goal points of the forward propagation stage. FMM is then implemented so that it matches start points to the maximum available number of the goal points resulting in an approximated map of the retinal blood vessels. Then back tracking procedure starts from the goal points towards start points drawing a single-pixel route from each end point to the start point. This route can be defined as the blood vessel tree of the retinal image. The algorithm was applied on images taken from the familiar DRIVE and STARE data sets. Pre-processing stage: here we try to have an initial approximation of the retinal blood vessel position. The importance of this stage comes from the huge amount of noise in fundus images and the convergence of the intensity valuesamong adjacent pixel at the first place, and from the high sensitivity of the FMM against contrast values of the adjacent pixels. FMM efficiency increases obviously as contrast increase. On the other hand, FMM becomes blind when contrast decrease and tends to cross existing constrains which affects vessel tracing badly.
To overcome these problems, the preprocessing stage included a median filtering procedure, followed by local numerical analysis of the retinal pixels in order to estimate most probable position of the blood vessel. This step is of great importance to magnify the FMM efficiency since it creates a kind of clear borders that can be considered during exploration process. This step can also reduce the number of iterations and time consumption required for exploration process. This initial estimation of the position of blood vessel is shown in figure (7). Applying FMM on pre-processed fundus images: In this stage, image enhancement is first performed to increase contrast between vessel pixels and non-vessel pixels, this enhances FMM performance since it enlarges variations between vessels pixels and background pixels. This stage consists of two main phases: In exploration phase, a forward propagation is implemented from stating points to the goal points. This phase gives every single point that have the probability of being among the vessel pixels.
In exploitation phase, the goal points of the exploration phase are considered as the starting points of exploitation phase, then a back tracking search is implemented in order to find pixel with least distance function values that gives its distancefrom starting points which are the goal points of the exploitation phase. At this phase, we don't take all the pixels that may belong to vessel. Only a single pixel width tracking line is taken into consideration. This line connects the start and goal pixels and provide us the desired tree of the retinal vessels. Figure (8) shows FMM main phases.

RESULTS AND DISCUSSIONS
FMM has been validated on many different images taken from the familiar DRIVE data set. It achieved a TPR of 80% to 93% and in some cases it achieved a TPR of 95%, with an iteration number ranges between 100 to 500 according to the OD position (OD distance from retina's center). The average time of each iteration is about 0.57 sec this results in a running time of 5-10 minutes. FMM hasalso been validated using STARE data set that contain PPM image with 700x605 pixels. For STARE data set, performance had been tested for images with different sizes using different performance measure like iteration number, time and TPR.  The number of iterations had also been studied comparing to the OD displacement from the image center using STARE data set. The results showed an incrimination in iteration number as the OD displacement increases although there is no obvious relationship to describe it, and this is because of the specific illumination properties of each image that have great effect on image quality. In this case, we used single start point with single endpoint to simplify the analysis process. The results showed that this method has a priority comparing to matched filters proposed by Chaudhuri [13] and model based method proposed by Jiang&Mojon [14]. Its performance is also comparable to other high accuracy methods such as rule-based method proposed by Martinez [15] and vessel analysis using morphological and structural properties of the vessels which hadan accuracy of 91.44% [5] in addition to multiscale line tracing method proposed in [3] which had an accuracy value of 92.9%. Speaking of processing time, our method shows a relatively long run time, but it is still close to other researches that declared using MATLAB like the research mentioned in [5] that achieved a processing time of 8 minutes for DRIVE images. While we had a processing time of 6 minutes for DRIVE images.

CONCLUSION
In this article, we introduced a new method for automatic retinal blood vessel tree segmentation based on mathematical analysis of the retinal texture in preprocessing stage thenan FMM stage was usedto trace vessels. We applied this method on images form DRIVE and STARE data sets using MATLAB. The proposed method achieved performance close to the other high accuracy algorithms used in this field. At the end, it is good to refer to the ability of improving this algorithm by searching for more edge-enhancing techniques and taking the advantages of spectral domain to segment arteries and veins and to isolate crossing points of the vessel tree.