Abstract
We propose an algorithm for vessel extraction in retinal images. The first step consists of applying anisotropic diffusion filtering in the initial vessel network in order to restore disconnected vessel lines and eliminate noisy lines. In the second step, a multiscale line-tracking procedure allows detecting all vessels having similar dimensions at a chosen scale. Computing the individual image maps requires different steps. First, a number of points are preselected using the eigenvalues of the Hessian matrix. These points are expected to be near to a vessel axis. Then, for each preselected point, the response map is computed from gradient information of the image at the current scale. Finally, the multiscale image map is derived after combining the individual image maps at different scales (sizes). Two publicly available datasets have been used to test the performance of the suggested method. The main dataset is the STARE project’s dataset and the second one is the DRIVE dataset. The experimental results, applied on the STARE dataset, show a maximum accuracy average of around 94.02%. Also, when performed on the DRIVE database, the maximum accuracy average reaches 91.55%.
1. Introduction
For decades, retinal images are widely used by ophthalmologists for the detection and follow-up of several pathological states [1–5]. Fundus photographs, also called retinal photography, are captured using special devices called “Charged Coupled Devices” (CCD), which are cameras that show the interior surface of the eye [6–10]. These images directly provide information about the normal and abnormal features in the retina. The normal features include the optic disk, fovea, and vascular network. There are different kinds of abnormal features caused by diabetic retinopathy (DR) such as microaneurysm, hard exudate, soft exudate, hemorrhage, and neovascularization. An example of retinal images obtained by fundus photography is given in Figure 1, where two retinal images are shown. The first one does not show any DR sign (Figure 1(a)) and the second one demonstrates advanced-DR signs indicated by color arrows (Figure 1(b)). However, the manual detection of blood vessels is very difficult since the blood vessels in these images are complex and have low level contrast [11]. Also, not all the images show signs of diabetic retinopathy. Hence, a manual measurement of the information about blood vessels, such as length, width, tortuosity, and branching pattern, becomes tedious. As a result, it increases the time of diagnosis and decreases the efficiency of ophthalmologists. Therefore, automatic methods for extracting and measuring the vessels in retinal images are needed to save the workload of the ophthalmologists and to assist in characterizing the detected lesions and identifying the false positives [12].
(a) Normal retina
(b) Abnormal retina
Several works have been proposed for detecting the 2D complex vessel network, such as single scale matched filter [13–15], multiscale matched filter [16], adaptive local thresholding [17], single-scale Gabor filters [18], and multiscale Gabor filters [19]. Cinsdikici and Aydin [20] put forward a blood vessel segmentation based on a novel hybrid model of the matched filter and the colony algorithm, which extracts vessels perfectly but the pathological areas can affect the result. In [21–23] authors adapted another approach which applied mathematical morphological operators. The suggested method in [21] proved to be a valuable tool for the segmentation of the vascular network in retinal images, where it allowed obtaining a final image with the segmented vessels by iteratively combining the centerline image with the set of images that resulted from the vessel segments’ reconstruction phase using the morphological operator. However, the inconvenience of this method is when a vessel centerline is missing, so the corresponding segmented vessel is normally not included in the final segmentation result. In [22], the authors proved that it was possible to select vessels using shape properties and connectivity, as well as differential properties like curvature. The robustness of the algorithm has been evaluated and tested on eye fundus images and on other images. Gang et al. [24] showed that the Gaussian curve is suitable for modeling the intensity profile of the cross section of the retinal vessels in color fundus images. Based on this elaboration, they proposed the amplitude-modified second-order Gaussian filter for retinal vessel detection, which optimized the matched filter and improved the successfulness of the detection. Staal et al. [25] explained a method for an automated segmentation of vessels in two-dimensional color images. The system was based on extracting image ridges that coincide approximately with vessel centerlines, where the evaluation was done using the accuracy of hard classifications and the values of soft ones. In [26], the authors presented a hybrid method for an efficient segmentation of multiple oriented blood vessels in colour retinal images. The robustness and accuracy of the method demonstrated that it might be useful in a wide range of retinal images even with the presence of lesions in the abnormal images. Dua et al. [27] presented a method for detecting blood vessels, which employs a hierarchical decomposition based on a quad tree decomposition. The algorithm was faster than the existing approaches. In the recent years, alternative approaches for an automated vessel segmentation have used the Hessian-based multiscale detection of curvilinear structures, which has been effective in discerning both large and small vessels [28–31].
In this paper, we propose a multiscale response to detect linear structures in 2D images. We will use the formulation, which was suggested in [36, 37]. The presented detection algorithm is divided into two steps. First, we present a flux-based anisotropic diffusion method and apply it to denoise images corrupted by an additive Gaussian noise. In order to extract only the pixels belonging to a vessel region, we use a Gaussian model of the vessels for interpreting the eigenvalues and the eigenvectors of the Hessian matrix. Then, we compute the multiscale response from responses computed at a discrete set of scales. The method has been evaluated using the images of two publicly available databases, the DRIVE database [34] and the STARE database [33]. Prior to analysing fundus images, we have used the green channel alone, since it gives the highest contrast between the vessel and the background.
2. Methodology
2.1. Preprocessing Technique
In the ocular fundus image, edges and local details between heterogeneous regions are the most interesting part for clinicians. Therefore, it is very important to preserve and enhance edges and local fine structures and simultaneously reduce the noise. To reduce the image noise, several approaches have been proposed using techniques such as linear and nonlinear filtering. In linear spatial filtering, such as Gaussian filtering, the content of a pixel is given by the value of the weighted average of its immediate neighbors. This filtering not only reduces the amplitude of noise fluctuations but also degrades sharp details such as lines or edges, so the resulting images appear blurred and diffused [24, 38]. This undesirable effect can be reduced or avoided by designing nonlinear filters. The most common technique is median filtering. With it the value of an output pixel is determined by the median of the neighborhood pixels. This filtering retains edges but results in a loss of resolution by suppressing fine details [39]. In order to perform this task, Perona and Malik (PM) [18] developed an anisotropic diffusion method, a multiscale smoothing, and the edge detection scheme, which were a powerful concept in image processing. The anisotropic diffusion was inspired from the heat diffusion equation by introducing a diffusion function, , which depended upon the norm of the gradient of the image:where and denote gradient operation and image intensity, respectively, is the divergence operator, and denotes the magnitude. The variable represents the spatial coordinate, while the variable is used to enumerate iteration steps in the discrete implementation. Perona and Malik suggested the following diffusion functions:where is a parameter of the norm gradient. In this method of anisotropic diffusion, the norm gradient is used to detect edges or frontiers in the image as a step of intensity discontinuity. To understand the relation between the parameter and the discontinuity value , can be defined as the following product , called the flow diffusion.(i)If , then and we have a filter pass-all.(ii)If , then and we obtain an isotropic diffusion filter (like a Gaussian filter), which is a low-pass filter that attenuates high frequencies.
The one-dimensional discrete implementation of (1) is given bywhere and .
The above result is generalized in -dimensional:if , and .
Up to now, the anisotropic diffusion has been defined as the case where the diffusivity is a scalar function varying with the location in the image. As described earlier, the PM diffusion (Figure 2) limits the smoothing of an image near the pixels with a high gradient magnitude (edge pixels). As the diffusion near an edge is very weak, the noise smoothing near the edge is also small. To address this, diffusions using matrices instead of scalars have been put forward [36, 40, 41], where the anisotropic diffusion allows the diffusion to be different along various directions defined by the local geometry of the structures in the image (Figure 3). Thus, the diffusion on both sides of an edge can be prevented while allowing the diffusion along the edge. This prevents the edge from being smoothed and then being removed during denoising.
The flux of the matrix diffusion (MD) form can be written as where is a positive definite symmetrie matrix that may be adapted to the local image structure, which can be written in terms of its eigenvectors and and eigenvalues and , as follows:Subsequently, the gradient vector field can be written asFollowing the eigenvalues and eigenvectors that we have chosen, different matrix diffusions can be obtained [36, 41]. The diffusion matrix proposed by Weickert et al. [41] had the same eigenvectors as the structure tensor, with eigenvalues that are a function of the norm of the gradient [41, 42]. In our work, we have used a 2D basis which corresponds, respectively, to unit vectors in the directions of the gradient and to the minimal curvature of the regularized (or smoothed) version of the image, which is the image convolved with a Gaussian filter with a standard deviation . This basis is of particular interest in the context of small, elongated structures such as blood vessels, where the minimal curvature holds for the axis direction orthogonal to the gradient. These directions are obtained as two of the eigenvectors of the Hessian matrix of the smoothed image: (further details are described in Section 2.3). Therefore, the eigenvectors are defined as follows: where is the gradient of the image convolved with a Gaussian filter with a standard deviation , gives an estimation of the vessel direction, and is its orthogonal. Also, we have used the eigenvalues in (6) as a diffusion function associated to each vector of the basis depending on the first order derivative of the intensity in this direction, instead of the traditional norm of the smoothed gradient. Furthermore, the diffusion can be decomposed as a sum of diffusions in each direction of the orthogonal basis and the divergence term can be written as [36]where and indicate the first order derivative of the intensity in the direction and the th diffusion function, respectively. Also, can be chosen to be any of the diffusivity functions from the traditional nonhomogeneous isotropic diffusion equation, which depends on the first order derivative of the intensity in this direction, as and , with , being only a diffusing function to allow smoothing in a direction. For further details, the reader could refer to [36, 43].
As in [36], we use a data attachment term with a coefficient which allows a better control of the extent to which the restored image differs from the original image (at ) and of the result of the diffusion process at convergence. The anisotropic diffusion equation becomes
In order to evaluate the denoising effects of the directional anisotropic diffusion , we have added a Gaussian white noise to each of the images in Figure 4. Once the diffusion method is applied to these noisy images, its effectiveness in reducing the noise is got by calculating the peak signal to noise ratio (PSNR) relative to the original image as follows: where and MSE is the mean-squared error which is written as where refers to the original image without noise and is the image after the denoising process.
(a)
(b)
(c)
(d)
(e)
(f)
The higher the PSNR is, the better the effect of the denoising is. Note that this measure does not necessarily imply that an image with a higher PSNR is also more visually gratifying. However, based on our experiments using the three test images with an additive white Gaussian noise, we can draw some observations. First, all the techniques we have tried have several parameters that must be selected carefully to obtain the best results. Since we have a “clean” original image, as well as one with noise, we can use the increment in the PSNR value to guide our choice of the parameters. These parameters and the obtained results are indicated in Tables 1, 2, and 3, where we can observe that for the images corrupted with an additive Gaussian noise, the DAD method performs better than the PM method. It gains a higher PSNR (, , and ) and a smaller MSE (, , and ) than the aforementioned three methods.
Figure 4 represents some of the best results for the different methods (GF, MF, PM, and DAD) on the presented three test images (Vessels, phantom, and Lena). For instance, the results recorded after applying the DAD method show that this latter improves much more the visual rendering of the image compared to other methods. As shown in the images of the first row, a DAD filter can effectively improves the quality of a noisy image and also well enhances edges and preserves more details than other filters. Indeed, the Gaussian filter smooths very strongly the planar areas which causes loss of information regarding the fine structures of the image, and it blurs the image. The Median filter, compared to the Gaussian filter, preserves edges but losses details. Comparing the results of the DAD method to those obtained by the PM diffusion in Figures 5 and 6, we can derive several observations. The denoising of PM diffusion model is sensitive to the value of the conductance parameter , and, therefore, smoothing is performed along ridges but not across a ridge line which causes enhancing the desired ridges as well as the noise. To be compared to the DAD diffusion filter, the diffusivity is a tensor-valued function varying with the location and orientation of edges in an image. So, when this filter is applied to a ridge line smoothing is performed along ridges as across a ridge line while preserving the details.
2.2. Multiscale Medialness
The general approach of multiscale methods is to choose a range of scales between and (corresponding to and ), which are discretized using a logarithmic scale in order to have more accuracy for low scales and to compute a response for each scale from the initial image [36, 43, 47]. The user specifies the minimal and maximal radius of the vessels to extract. Thus, the computation of the single scale response requires different steps. First, a number of points are preselected using the eigenvalues of the Hessian matrix. These points are expected to be near a vessel axis. Then, for each preselected point, the response is computed at the current scale . The response function uses eigenvectors of the Hessian matrix of the image to define at each point an orientation orthogonal to the axis of a potential vessel that goes through . From this direction, the two points located at an equal distance of in the direction , noted and (Figure 7). The response at is taken as the maximum absolute value, among these two points, of the first derivative of the intensity in the direction:where is the unitary vector of the direction , that is, , and is the gradient of the image at the scale . is obtained by the convolution with the first derivative of a Gaussian function of the standard deviation , where multiplying the derivatives by ensures the scale invariance property and allows comparing the responses obtained from different scales. The gradient vector can be computed by a bilinear interpolation for better accuracy, which is especially needed when looking at small vessels [37, 39].
A vessel of a radius is detected at a scale , so we use the scales corresponding to each radius for the multiscale processing. For a fixed scale , we calculate a response image where is the initial image. Then we calculate the multiscale response for the image which is the maximum of the responses over scales: for each point and a range of scale:This response can be interpreted as an indicator that the point belongs to the center line of a vessel, and can be interpreted as an indicator that the point belongs to the center line of a vessel with a radius . Finally, this response is normalized to give a multiscale response that combines interesting features of each single scale response.
One difficulty with the multiscale approach is that we want to compare the result of a response function at different scales, whereas the intensity and its derivatives are decreasing scale functions. So far, all considerations have been made at a single scale defined by the scale parameter . In his work, about scale space theory, Lindeberg and Fagerström [48] showed the need for a multiscale analysis to take the varying size of objects into account. He also showed the necessity of normalizing the spatial derivatives between different scales. Thus, the normalized vesselness response is obtained by the product of the normalization term and the final vesselness: The parameter can be used to indicate the preference for a particular scale (Figure 8). If it is set to one, no scale is preferred. Besides, the multiscale response is got by selecting the maximum response over a set of different scales between and .
(a)
(b)
(c)
2.3. Extraction of Local Orientations
The proposed model assumes that the intensity profile of the vessels in the cross section is Gaussian (Figure 9). This is a common assumption that it is employed in numerous algorithms [28, 35, 49]. It is also commonly assumed that the intensity does not change much along vessels [49–51]. Recently, the Hessian matrix could be used to describe the local shape characteristics and orientation for elongated structures [35, 52]. The eigenvalues of this matrix, when the gradient is weak, express the local variation of the intensity in the direction of the associated eigenvectors. Subsequently, we assume that we want to characterize the dark vessels (low intensity) on a white background (high intensity).
(a)
(b)
Let us denote and as the eigenvalues of the Hessian matrix with and , being their associated eigenvectors (Figure 10). For a linear model with a Gaussian cross section, the vessel direction is defined by the eigenvector with the smallest eigenvalue at the center of the vessel, but it is less determined at the contours because both eigenvalues of the Hessian matrix are zero.
(a)
(b)
(c)
To summarize, for an ideal linear structure in a 2D image,
In retinal images, some large vessels may have a white line in their center and some elongated and disjoint spots (Figures 11(a), 11(b), and 11(c)); accordingly, the vessels do not invalidate the Gaussian profile assumption. So, such lines are usually lost after the preselection of vessel pixels using the Hessian eigenvalue analysis and classified as background pixels. Therefore, the responses of the gradient magnitude are a task which is of particular importance in improving the detection vessels (Figure 11). The experimental results are demonstrated in Figure 11, which shows hand labeled “truth” images, and segmented images obtained, respectively, by the Hessian eigenvalue analysis and the gradient magnitude. From these results we can deduce that responses based on the gradient magnitude can availably detect white lines as vessel pixels an removes some noise spots.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
3. Results
In this section, the proposed method has been evaluated on two publicly available retinal image databases, the STARE database [33] and the DRIVE database [25]. The STARE dataset contains twenty fundus colour retinal images, ten of which are from healthy ocular fundi and the other ten are from unhealthy ones. These images are captured by a Topcon TRV-50 fundus camera at a 35 Field Of View (FOV), which have digitized with a 24-bit gray-scale resolution and a size of pixels. The dataset provides two sets of standard hand-labeled segmentations, which are manually segmented by two eye specialists. We create for this dataset a binary mask of the gray channel of the image using a simple threshold technique (Figure 12). We adapt the first eye specialist hand labelled as the ground truth to evaluate our vessel detection technique. The DRIVE dataset consists of 40 fundus ocular images, which have been divided into a training set and a test set by the authors of the database. These images are captured by the Canon CR5 camera at 45 FOV, which have been digitized at 24 bits with a resolution of pixels. The dataset also gives two sets of standard hand-labeled segmentations by two human experts as a 9-ground truth.
(a)
(b)
The first expert hand labelled segmentation has been adapted as a ground truth to evaluate segmentation techniques on both STARE and DRIVE datasets. It is a common practice to evaluate the performance of retinal vessel segmentation algorithms using receiver operating characteristic (ROC) curves [25, 35]. An ROC curve plots the fraction of pixels correctly classified as vessels, namely, the true positive ra te (TPR), versus the fraction of pixels wrongly classified as vessels, namely, the false positive rate (FPR), by varying the rounding threshold from to (Figure 13). The closer the curve approaches the top left corner, the better the performance of the system. In order to facilitate the comparison with other retinal vessel detection algorithms, we have selected the value of the area under the curve (AUC), which is for an ideal system.
(a)
(b)
(c)
To measure the performance of the proposed enhancement filter, we ran our multiscale analysis filter with the following set of parameters:(i), , , and the minimal and maximal radii used in this application are and , discretized using scales;(ii)the parameter set to one to indicate no scale is preferred;(iii)the value is a constant threshold on the norm of gradient on the image;(iv) is the number of iterations for the anisotropic diffusion filter.
The computing time of our algorithm for an image of the STARE database is about seconds, including anisotropic diffusion filtering, and about the same time for the DRIVE database. The implementation of the filter has been done in MATLAB, on a personal computer with a Intel Core Duo processor and 4 GB of memory. In the first experiment, we apply a preprocessing task such as filtering data with an anisotropic diffusion version, cited above, in order to remove or at least reduce noise. The DAD filter denoises the original image by preserving edges and details. To show that the segmentation works better with anisotropic diffusion, Figure 14 presents a segmentation result before and after the application of the anisotropic diffusion scheme. In this figure, we show the improvements provided by the DAD model, which tends to remove noise effects and, unfortunately, smaller objects. So, it preserves efficiently the vessels while making the background more homogeneous.
(a)
(b)
(c)
(d)
On the other hand, for computing the response, it is possible to retain the mean of the two calculated values (the gradient of the two points located at an equal distance from the current point), like in the 3D case proposed by [36], or the minimal calculated value in the 2D case [37]. We prefer retaining the maximum of these two values. Figure 15 shows a synthetic image which consists of 100 × 100 pixels with an 8-bit resolution. We have chosen this image because it contains an object close to the vessel form. The latter figure shows the segmentation results by maximum, average, and minimum response functions. We note that for the case of minimum or average responses, the ring is not completely detected like in the original image, since we can see it has been missing pixels belonging to the edges, in contrast to maximum case where the extraction of the ring is complete. Table 4 presents the calculated with our method for the test set of the STARE database, using the green channel images. As given in the table, the experimental results show that the maximum model performs much better than the average or minimum model .
Figure 16 presents the obtained response image of a real retinal image, where four scales have been used for radii of vessels ranging from to 7: . This figure shows that small and large vessels can be better distinguished in the maximum case than the minimum or average ones.
Although the contrast is not very high in the original figure (Figure 14(a)), the method detects most vessels, over a large size range. For example, in Figure 17, an image of the retinal tree vasculature is presented, where different responses recorded at increasing scales are represented. The last image shows a quite good performance of the vessel subtraction. Yet Figure 18 proves that it is possible to design a system that approaches the performance of human observers.
In order to evaluate the suggested method, experiment results of the 20-image sets of the STARE database are shown in Table 5. In Table 6, our method is compared to the most recent methods in terms of , , and maximum accuracy average where the maximal accuracy indicates how to extract a binary image that matches the vessel images to a high degree. The accuracy is estimated by the ratio of the sum of the number of correctly classified foreground and background pixels, divided by the total number of pixels in the image. In this latest table, the performance measures of Staal et al. [25], Zhang et al. [14], Mendonça and Campilho [21], Chaudhuri et al. [13], Martinez-Perez et al. [45], and Hoover et al. [35] have been reported by their original papers. In addition, these performance results are the average values for the whole set of images, except the method of Staal [25] which used out of images of the STARE images, among which ten were healthy and nine were unhealthy. Table 5 presents our results on all 20 images in the STARE database, estimated using the hand-labeled images set of the first human expert designated as a ground truth. The estimated experimental results are the average corresponding to an of around and a maximum average accuracy . The results show that our method has a competitive maximum average accuracy value where it performs better than the matched filter [13] and remains close to the others.
The results of the proposed method are also compared with those on twenty images from the DRIVE database, and the result is depicted in Table 7. The hand-labeled images by the first human expert have been used as ground truth. The experimental results show an around of . Also, we have compared the performance of the suggested technique with the sensitivities and specificities of the methods cited in Table 7. It has been found that for the DRIVE database the method has provided a sensitivity of and a specificity of . We have shown that the proposed method performs well with a lower specificity even in the presence of lesions in the abnormal images.
4. Conclusion
The purpose of this work is to detect linear structures in real retinal images in order to help the interpretation of the vascular network. We put forward to combining an anisotropic diffusion filter to reduce the image noise with a multiscale response based on the eigenvectors of the Hessian matrix and on the gradient information to extract vessels from retinal images. The main advantage of this technique is its ability to extract large and fine vessels at various image resolutions. Furthermore, the directional anisotropic diffusion plays a vital role in denoising images and in decreasing the difficulty of vessel extraction especially for thin vessels. Our first results show the robustness of the method against noise as well as its applicability to detect blood vessels. The is used as a performance measure, and the values achieved with our algorithm are competitive compared to the existing methods. Therefore, from the experimental results, it can be seen that the number of classified pixels has been slightly lower compared to the other methods using the same database mainly due to the weakness of blood vessels, causing missing vessels, and also because of lesions, resulting in a detection error. Also, the retinal images suffer from nonuniform illumination and have a poor contrast. Thus, to avoid wrong classified pixels or miss classified ones, caused by an occasional false measurement, this system can very well be improved in the future with adding, for instance, some postprocessing tasks to reach more accurate measurement for blood vessels.
Conflict of Interests
The authors declare that there is no conflict of interests regarding to the publication of this paper.