Keywords

1 Introduction

Automatic procedures to analyze retinal blood vessels are useful in ophthalmology to allow an early diagnosis of a number of diseases, such as diabetic retinopathy, arteriosclerosis, hypertension, cardiovascular diseases and stroke [2, 3]. The automatic analysis of the structure of retinal vessels is also useful in biometrics [4], since the structure is different for each individual and even for the left and the right eye of the same person. In the literature, segmentation techniques based on matched filters, e.g., [5], wavelet transform, e.g., [6, 7], line detectors, e.g., [8,9,10,11,12], and morphological image processing [8, 9, 13,14,15,16] are available. In this paper we present an unsupervised retinal blood vessels segmentation method based on directional information. Each pixel of the green channel G of a RGB retinal image is assigned a direction out of twelve possible discrete directions. Two different threshold values are then employed to roughly segment the so obtained direction map DM G . The foreground in the under-segmented image G u , obtained in correspondence with the higher threshold value, includes pixels that can be most possibly interpreted as vessel pixels. The foreground in the over-segmented image G o , obtained in correspondence with the lower threshold value, includes several more pixels than G u . Some pixels of G o are detected as missing vessel pixels of G u , and are added to G u to improve the quality of the segmentation. The method has been tested on the DRIVE database [1] producing satisfactory results, and its performance has been compared to that of other unsupervised methods.

2 Building the Direction Map DM G

We work with the green channel G of RGB color retinal images, as done by the majority of authors in the field, since G is the channel characterized by the highest contrast. Gray-levels in G are in the range [0,255], where lighter pixels are those with larger gray-level values. Vessel pixels are thin structures whose pixels are generally darker than their neighboring non-vessel pixels and are aligned along different directions.

The direction of any pixel p of G can be computed by taking into account the gray-levels of the pixels in an n × n window W centered on p, so as to build the Direction Map DM G . Selecting the proper value for n is crucial to obtain a correct segmentation: W should be large enough to include both vessel and non-vessel pixels, even when centered on the most internal pixels of the thickest vessels, so as to have appreciable variations of gray-levels within the window; on the other hand, W should not be too large so as to avoid the inclusion of vessel pixels belonging to close vessels. By taking into account the average width and distance of vessels in the DRIVE database, we set n = 7. Since in an n × n window, 2 × (n−1) discrete straight lines can be built, 12 directions with an angle of 15° between any two successive directions are obtained. See Fig. 1 top, where the twelve directions are shown in different colors.

Fig. 1.
figure 1

The twelve directions d i , top, and the twelve directional templates T_d i , bottom. (Color figure online)

Full size image

To build DM G we use the twelve 7 × 7 directional templates T_d i , i = 1,2,…,12, shown in Fig. 1 bottom. In the i-th template, the pixels aligned along any direction out of d i-1 , d i , and d i+1 are set to 1, and all the remaining pixels are set to 0. Of course, for i = 1 (i = 12) d i-1 , d i , d i+1 are respectively d 12 , d 1 , d 2 (d 11 , d 12 , d 1 ).

Given two arrays DIFF and DM G with the same size as G and initially empty, for each pixel p of G and for any directional template T_d i , i = 1,2,…,12, we compute:

  • the arithmetic mean, Ad i (p), of gray-levels of pixels of G matching 1’s in T_d i ,

  • the arithmetic mean, NAd i (p), of gray-levels of pixels of G matching 0’s in T_d i ,

  • the difference Δi(p) = NAd i (p) – Ad i (p).

The direction d i for which Δi(p) has the maximum value M(p) is taken as direction of p; the score M(p) and the index i are respectively stored in the homologous position of p in DIFF and DM G . For the green channel G of the image 05_test of the DRIVE database, shown in Fig. 2 left, the obtained DM G is shown in Fig. 2 middle, where colors correspond to the directions as in Fig. 1 top.

Fig. 2.
figure 2

The image G, left; the direction map DM G before, middle, and after, right, direction updating. (Color figure online)

Full size image

Actually, some pixels in DM G have been assigned a direction different from the one along which they appear to be aligned. Thus, an updating of DM G is done for each of the twelve directions and by using an auxiliary array AUX with the same size as G and whose pixels are initially set to 255. For the i-th direction d i , the pixels of DM G with assigned direction d i are set to 0 in AUX. Thus, AUX becomes a gray-level image with only two values, where the foreground consists of the pixels with direction d i , and the background includes all the remaining pixels. The idea is to build the direction map of AUX, DM AUX , by following the procedure described to compute DM G and to compare the direction assigned to any pixel p in the two maps. Only pixels for which the same direction is obtained in DM G and DM AUX are confirmed as foreground pixels in DM G . All other pixels of DM G are assigned to the background. Actually, when transferring in AUX pixels with direction d i , also the pixels that in DM G where characterized by the directions d i-1 and d i+1 are set to 0 in AUX, while updating is done only for the pixels with direction d i . The purpose is to avoid that only a few sparse pixels of AUX are foreground pixels, for which almost all directions would be equally possible, so biasing the correct updating of the directions. The updated DM G is shown in Fig. 2 right, where white pixels are those assigned to the background.

3 The Segmentation Method

Segmentation can be achieved by assigning to the background any pixel p whose score M(p) in DIFF is smaller than a threshold θ, which is set based on the directional features of the processed image. To this aim, we compute the arithmetic mean m of the scores M(p) of all pixels that are currently foreground pixels, i.e., pixels that in DM G are different from zero, and set θ equal to a given percentage of m.

Actually, we use two different percentages of m, so as to achieve two different values for θ, θ_u = u × m and θ_o = o × m with u > o, leading to two different segmented images, G u and G o . When the higher threshold value is adopted, only the pixels with higher probability to be true vessel pixels are assigned to the foreground in G u . The resulting image is generally under-segmented, since pixels belonging to slightly lighter vessels or to capillaries may not survive thresholding. With the lower threshold value, also pixels with lower probability to be vessel pixels are selected, so achieving an over-segmented image G o . The number of pixels that are not true vessel pixels is much smaller in G u than in G o . For the running example, the under-segmented image G u and the over-segmented image G o , obtained for u = 1.8 and o = 0.9, are respectively shown in Fig. 3 left and middle. The image difference G o-u includes the pixels in the foreground of G o but not in the foreground of G u . See Fig. 3 right.

Fig. 3.
figure 3

The image G u , left, the image G o , middle, and the image difference G o-u , right.

Full size image

Linking and Cleaning

The image difference G o-u is rather noisy. In fact, the low threshold value θ_o = o × m adopted to obtain G o guarantees the detection of almost all the true vessel pixels, but also causes the wrong assignment to the foreground of a number of false vessel pixels. Thus, we consider as necessary some cleaning of G o-u , based on removal of small size components. Before doing the size-based cleaning, we link with each other small close components that we interpret as parts of vessels resulting in distinct connected foreground components of G o-u since some of their pixels were selected as foreground pixels when adopting both the lower and the higher threshold value.

For the linking process we consider components of G o-u with size smaller than a given maximum value max, since only small size components risk to be removed during the successive cleaning task. On the other hand, we do not consider components with very little size, say with size smaller than a given minimum value min, since these components most possibly consist of false vessel pixels. Thus, we consider for the linking process every component of G o-u whose size s is such that min ≤ s ≤ max. We have experimentally found that the better performance of the method is achieved in the average by setting min = 16 and max = 150. We use an iterated expansion process. At each iteration, the background neighbors of any pixel p in the components of G o-u selected for the linking process are assigned to the components provided that they are foreground pixels in G o , and share with p the same direction. The first requirement guarantees that linking pixels were at least tentatively selected as vessel pixels in the over-segmented image; the second requirement guarantees that linking regards exclusively components that actually are part of vessels aligned along given directions. The number of iterations depends on the maximal distance between two components to be connected. Since G o includes many false vessel pixels, we limit the number of iterations to two, which means that we can link to each other only components with a maximal distance of at most four pixels. The expansion process is followed by an iterated topological removal process, aimed at assigning again to the background all pixels added by the expansion process except those that favored linking of components. When no more pixels are removed, the surviving added pixels are assigned to the foreground in G o-u .

Size-based cleaning is performed, by removing from G o-u all components with area smaller than a threshold σ. We have found that the best performance is obtained in the average for σ = 64. The result obtained after linking and cleaning can be seen in Fig. 4 left. Different colors denote the identity labels of the components of G o-u .

Fig. 4.
figure 4

The image G o-u after linking and cleaning, left, expansion, middle, and shrinking, right. (Color figure online)

Full size image

Selection of Foreground Components Consisting of True Vessel Pixels

We use the main feature characterizing vessels, i.e., the fact that vessels have linear structure, to distinguish in G o-u the components to be maintained in the foreground from those to be assigned to the background. We perform a small number of iterations of an expansion process, set to three in this work, during which thin holes and concavities possibly present in any component of G o-u are filled in, without creating an unwanted merging of foreground components. See Fig. 4 middle showing the resulting G o-u after the three iterations. Then, we perform an equal number of iterations of topological shrinking, during which pixels added to any component by the expansion process are removed provided that they have at least a neighbor in any other component, including the background. The resulting image G o-u is shown in Fig. 4 right. It can be observed that while components originally characterized by linear structure are not remarkably modified by the expansion/shrinking process as regards their size, components with a more complex structure, erroneously assigned to the foreground, have at the end of the process a significantly larger size. Actually, besides the changes in size, we also take into account the changes in maximal width of the components of G o-u . This last feature is easily measured as the maximal value in the distance transform of each component before and after the expansion/shrinking process. Let W in and W exp be the initial maximal width and the maximal width after the expansion/shrinking process, respectively. Moreover, let A in and A exp be the initial area and the area after the expansion/shrinking process. We assign to the background any component for which it results A exp /A in  ≥ 0.40, or it is 0.30 ≤ A exp /A in  < 0.40 and W exp /W in  ≥ 2. All other components remain in the foreground.

All components of G o-u surviving the expansion/shrinking process are recognized as consisting of vessel pixels. The pixels that belonged to these components before the expansion/shrinking process are transferred into G u . The currently obtained G u can be seen in Fig. 5 left.

Fig. 5.
figure 5

The image G u after adding the contribution of G o-u , left; G u after recovery of missed vessel pixels and removal of small components, middle left, G u after foreground concavities and holes filling, middle right, and G u after cleaning of the circular boundaries of retina and optical disc, right.

Full size image

Improving Segmentation

To recover missed vessel pixels, we apply to the components of G u a slightly modified version of the process adopted to link close components of G o-u . As done when processing G o-u , background neighbors of pixels in components of G u are added to the components only if they also exist in G o . Differently from the process applied to G o-u , the number of iterations is not fixed a priori; moreover, we have some tolerance on the direction of the pixels to be added. In detail, the expansion process terminates when no more pixels can be added; at each iteration, if d i is the direction of a pixel p in a given component of G u , any neighbor of p is parallelwise added to G u , provided that its assigned direction is any out of d i-1 , d i , d i+1 . Finally, differently from the process done on G o-u , the expansion process is not followed by any iterated topological removal and we accept as vessel pixels all the recovered pixels. Finally, components of G o-u that notwithstanding the recovery process are still characterized by small size are removed. We use the same threshold σ = 64 already adopted for size-based cleaning of G o-u . The resulting image is shown in Fig. 5 middle left.

By taking into account the average width and distance of the vessels in the DRIVE database, we found adequate to fill foreground concavities up to 2-pixel wide. Since the concavities filling process is the same used in [12], for space limitation we do not describe again it here. Very small size holes (consisting of at most 16 pixels in this work) are also filled in. The resulting image G u is shown in Fig. 5 middle right.

The last process is devoted to cleaning of the circular boundaries of retina and of optical disc. The masks available in the DRIVE database are used to extract only the circular part of G u corresponding to the retina. However, some pixels in proximity of the boundary of the circular mask may have been erroneously interpreted as vessel pixels. Thus, we dilate the circular boundary of the mask by means of a structuring element of radius 20. Foreground pixels reached by the dilation process are marked as removable. Then, an iterated un-marking process is accomplished that removes the marker from any neighbor q of an un-marked vessel pixel p with direction d i if the direction of q is compatible with the direction of p, i.e., the direction of q is any out of d i-1 , d i , d i+1 . The number of iterations of the un-marking process is equal to 20, i.e., to the radius of the structuring element. Pixels that at the end of the un-marking process are still marked are assigned to the background.

To remove pixels erroneously detected as vessel pixels within the optical disc, we first identify in G a region of interest, ROI, that includes the optical disc. It can be noted that the optical disc includes pixels with remarkably lighter gray-levels, and its maximal diameter and position are predictable in the DRIVE database. Thus, we use a sliding window of size 121 × 121 which swipes the image in the rectangular area where the optical disc is expected to be located. The pixels in the sliding window with gray-level larger than 80% of the maximal gray-level in G are counted while the window swipes the image. The position of the window in correspondence of which the number of counted pixels has the highest value defines the ROI. We observe that some vessels exist in the optical disc and these are characterized by rather dark gray-level and almost horizontal direction. Thus, out of all pixels of the ROI that have been assigned to the foreground of G u , we assign to the background those with gray-level larger than a suitable percentage, perc, of the arithmetic mean of the gray-levels within the ROI and with direction i which is not compatible with the horizontal direction. We have obtained the best performance by setting perc = 90% and by considering compatible with the horizontal direction d i , characterized by i = 7, any direction out of d 5 , d 6 , d 7 , d 8 , and d 9 . The final segmentation result is shown in Fig. 5 right.

4 Experimental Results and Concluding Remarks

The suggested segmentation method has been checked on the 40 images of the DRIVE database, by using as ground truth the manually segmented images included in the database. Actually, two manually segmented images are available for each of the 20 retinal images forming the test set, and one manual segmentation for the remaining 20 images forming the training set. Since our segmentation method is unsupervised, we have not done any difference between the test set and the training set.

For a qualitative evaluation of the method, see Fig. 6, where the results for the four images 09_test, 19_test, 27_training, and 32_training are shown.

Fig. 6.
figure 6

Some examples, top, and the segmentation results, bottom.

Full size image

To quantitatively evaluate the method, let TP and TN count the pixels correctly classified as vessel pixels and as non-vessel pixels, and FP and FN the pixels incorrectly classified as vessel pixels and as non-vessel pixels, and compute:

  • Accuracy = (TP + TN)/(TP + FN + TN + FP)

  • Sensitivity = TP/(TP + FN)

  • Specificity = TN/(TN + FP)

The average values of Accuracy, Sensitivity and Specificity have been computed for the 20 test images with respect to the first ground truth to compare the performance of our method with that of other 6 unsupervised segmentation methods in the literature, [8, 11, 12, 17,18,19]. See Table 1. Our method is slightly better as regards Accuracy and Specificity, and inferior for Sensitivity, but generally produces results qualitatively more satisfactory.

Table 1. Performance comparisons.
Full size table