Segmentation of blood vessel structures in retinal fundus images with Logarithmic Gabor filters Sebastian Gross, Monika Klein, Dorian Schneider Institute of Imaging & Computer Vision, RWTH Aachen University, 52056 Aachen, Germany Internal Medicine III, University Hospital Aachen, 52057 Aachen, Germany [email protected] Abstract. The analysis of blood vessel structures in the retinal fundus imagesisimportantforthediagnosisofmanydiseases.Vesselsegmenta- tioncanassistinthedetectionofpathologicalchangeswhicharepossible indicatorsforarteriosclerosis,retinopathy,microaneurysmsandmacular degeneration. In this article, two approaches to blood vessel segmentation are pre- sented. Both of them are based on the evaluation of phase symmetry information using complex logarithmic Gabor wavelets. In the first ap- proach, a phase symmetry filter is combined with the front propagation algorithmfastmarching,thesecondmethodusesahysteresisthreshold- ing step. Theapproacheshaveshownexcellentresultsforthevesselsegmentation oncolonpolyps.Althoughtheywereadaptedtostructuresinretinalfun- dusimaging,neithereyespecificknowledgenorsupervisedclassification methods are used. For high comparability with previous publications in the field, the al- gorithms are evaluated on the two publicly available image databases DRIVE and STARE. The hysteresis thresholding approach which per- formsslightlysuperiorachievesanaverageaccuracyof94.92%(sensitiv- ity:71.22%,specificity:98.41%)fortheDRIVEand95,65%(sensitivity: 71.87%, specificity: 98.34%) for the STARE database. 1 Introduction In modern medicine, the analysis of retinal fundus images is of particular im- portance to various medical diagnoses [1]. Pathological changes in retinal blood vessel systems can be an indication for several eye diseases like retinopathy, mi- croaneurysms,andmaculardegeneration.Inaddition,theriskforcardiovascular diseases like hypertension, arteriosclerosis, and stroke [2] can be estimated with minimal effort, since the observation of the ocular fundus is the easiest way to examine finest blood vessels, where the first signs of vascular damage become visible. 2 S. Gross, M. Klein, D. Schneider Thesegmentationofretinalbloodvesselscanbehelpfuleithertoremovethe vesselsfromtheimagetodetecthemorrhagesorfortheexaminationofthevessel tree itself. Features like cross sections of arteries and veins, contrast, curvature, andconcentrationofthevesselsaswellasthesizeofareaswithinsufficientblood vessel support can be used to evaluate the condition of the blood circulation. Furthermore, segmentation can assist in localizing the optic nerve or ease the registration. Inthispaper,twoapproachestoautomaticvesselsegmentationarepresented and compared to previously published methods. They have, in similar versions, shown excellent results for the segmentation of vessels [3–5] Although they have been adapted to retinal fundus images, neither eye specific knowledge nor a supervised classification methods were used. The text is structured into multiple sections. Section 2 gives an overview of previousworkinthisarea.Insections3and4,ourmethodsforbloodvesselseg- mentation are described and the used image databases and evaluation methods are highlighted. In section 5, results are presented and section 6 gives a short conclusion and discussion on the topic. 2 Overview of related work A common approach to blood vessel segmentation is the analysis of ridge infor- mation[6–10],wheretheeigenvectorsoftheHessianmatrixareusedtodetermine thedirectionofhighestorientation.Matchedfiltersare frequentlyused[8–12]as well. To reduce computational efforts, a lot of publications propose the Hessian matrixtodiscovertheassumedvesseldirectionandsubsequentlocalizationwith matched filters [8–10]. Further approaches are feature extraction and classifica- tion [13] (for example based on local binary patterns or morphological opera- tions), line detection [14], line segment detection [6] with subsequent segment classification based on forward feature selection [15], Gabor wavelets [16,17], or thresholding[7,18].Numerouspublicationsapplysupervisedclassificationmeth- ods, in which, based on defined features, each pixel is classified into one of the two possible classes (blood vessel or background). 3 Algorithms for blood vessel segmentation We propose to apply an approach developed for blood vessel segmentation on colon polyps [3–5] to fundus images. A graphical overview of the approach is given in Fig. 1. 3.1 Pre-processing As the green image channel provides the highest contrast between blood vessels and background, it is chosen for the segmentation [14,16]. To minimize background variations, background equalization and contrast enhancement on the green color channel are performed with an 11×11 tophat Retinal blood vessel segmentation with logarithmic Gabor Filters 3 Fig.1. Schematic process chain of the algorithms 4 S. Gross, M. Klein, D. Schneider filter[19]aspresentedin[3].Astructureelementisselectedandrankorderfilter operations [20] are performed with this structure element. First, the element is aligned horizontally and the maximum of the underlying pixels is stored at the center pixel’s position in a new image. Each of the following filters uses the preceeding filter results as input image. The next step is filtering with a vertically aligned structure element. Subse- quently, a second horizontal filter operation is performed. However, this time, the minimum value is stored at the center pixel’s position. Finally, a vertical minimum filter operation completes the tophat filter operation. The tophat filter results are used to scale the original image. To this is end, each pixel is divided by the value of the corresponding pixel from the tophat filter results. The result is then subtracted from the original green channel. Additionally, a small median filter of size 3×3 is applied to the image for noise reduction. The image borders are padded by periodic repetition to allow filter operations on the whole image [16]. 3.2 Soft vessel segmentation Examining the cross-section of the vessel, the gray level is increasing symmetri- callytowardsthevessel’scenterline.Forthedetectionandanalysisofsymmetric and asymmetric structures, Peter Kovesi presented an algorithm, which eval- uates the local symmetry information in gray-level images based on complex logarithmic Gabor functions [21,22] as suggested by Field [23]. An advantage compared to other algorithms is the independency of contrast variations. Furthermore, no segmentation of the objects prior to the symmetry analysisisnecessary.Kovesi’salgorithmuseslogarithmicGaborwaveletstoeval- uate the local phase symmetry, which corresponds to mirror symmetry in the spatial domain. Gabor wavelets are described by complex sine functions multiplied with a Gaussian in the spatial domain gω0(x)= ωσ022e−2ωσ022x2ejω0x (1) wherexistheinputimage’sgrayvalueandσistheGaussianstandarddeviation. In the spectral domain, they are represented by a Gaussian with the sine wave’s frequency w as central point. 0 −σ2 G(ω)=e2ω02 e(ω−ω0)2 (2) As Gabor wavelets are limited in the spatial as well as the spectral domain and provide an optimal compromise between spatial and frequency resolution, theyarewellsuitedaslocalbandpassfilters.Byadaptationofthewaveletscale n,structuresofdifferentsizescanbefound.Forthedetectionofsmallstructures, a low scale n is used, which provides high resolution in the spatial domain. For Retinal blood vessel segmentation with logarithmic Gabor Filters 5 larger structures, a higher scale n is applied. This provides a small bandwidth in the frequency domain, which is essential for noise suppression. Scalingfactorsn=n ·c·mforKovesi’simplementationaredescribedbythe 0 minimum wavelet size n =3, a constant factor c=2.1 and m∈{1...m } 0 scales with m = 3 being the number wavelet scales (detail images: n = 2 and scales 0 m = 4). The parameters are determined empirically. To this end, several scales tests have been performed and results were evaluated by experienced observers. For image processing, however, and especially for the detection of sharper structures, an optimal resolution in the spatial domain with larger bandwidth in the frequency domain is often preferable. A solution is the use of logarithmic Gabor wavelets, which have a Gaussian frequency response over a logarithmic instead of a linear frequency scale. The transfer function on a logarithmic scale is given by G(ω)=e(−log(ωω0)2)/(2(log(ωσ0))2). (3) as recommended by Field [23], where ω is the center frequency. κ is used to 0 keep κ/ω constant while ω changes due to scaling of the filter (denoted by n). 0 0 Compared to normal Gabor wavelets, a larger bandwidth is possible while still maintaining a zero mean. InordertoapplylogarithmicGaborwaveletsindifferentdirections,thefilter is combined with directional filters in six different directions (resolution 30◦). In addition, they are multiplied with a large circular high pass filter. This allows using a rectangular filter without the problem of varying maximum frequencies in different directions. In each direction, filters with multiple scales are applied. All filter results are added up to gain the final filter response. LikeGaborwavelets,logarithmicGaborwaveletsconsistofasymmetricreal andanasymmetricimaginarypart[24].Therefore,byconvolutionwiththegreen channel image, the real part responds to symmetric frequencies, whereas the imaginary part detects asymmetric structures in the image. The preprocessed images as well as the real and imaginary part of the log Gabor filter response are shown in Fig. ??. Real and imaginary part of the filter response with scale n are calculated by convolution of the symmetric and asymmetric Gabor filter pair with the input image’s green channel intensities and represented by e (x) and o (x). n n (cid:112) A = e (x)2+o (x)2 (4) n n n describes the magnitude. Kovesi proposes (cid:80) (cid:98)|e (x)|−|o (x)|−T(cid:99) sym(x)= n n n (5) (cid:80) A (x)+(cid:15) n n to be calculated by subtracting the absolute values of e (x) and o (x) for each n n scale, adding the results together, and normalizing the sum to the magnitude A . n Furthermore, a noise component T is subtracted to suppress symmetry re- sponses caused by noise. T is estimated by the mean squared filter response at 6 S. Gross, M. Klein, D. Schneider thesmallestscale.SinceA isofnon-negativevalue,acomponent(cid:15)>0isadded n topreventdivisionbyzero.Applyingthisapproach,however,thedetectedblood vessel responses appear much too thin in the results. An investigation shows that the responses to the asymmetric filter often eliminatethesmallvesselscompletely.Furthermore,thiscancompletelyprevent the detection of small vessels. Both filter results for overview and detail images can be seen in Fig. 2. Bloodvessels,however,areingeneralwellcontrastedsymmetricalstructures inretinalfundusimages.Sofkapostulatesthatbloodvesselstructuresarepiece- wise symmetric and single edges do not occur [9,10]. Therefore, we propose to ignore the imaginary part and use (cid:80) (cid:98)|e (x)|−T(cid:99) sym(x)= n n (6) (cid:80) A (x)+(cid:15) n n as modified approach. The results of the phase symmetry filter contain the segmentation as well as an orientation matrix containing the direction of the highest filter response for each pixel. In a further step, each pixel above a heuristically selected threshold which is connected to at least seven such pixels is segmented. Additionally, the orientation output of the symmetry filter is used to remove all pixels, which are not connected to a minimum number of pixels with the same orientation. 3.3 Binary classification Forbinaryclassification(binarysegmentation),twodifferentstrategiesarecom- pared: the front-propagation algorithm fast marching [25] and a method using a hysteresis thresholding step. Inthefirstcase,thephasesymmetryfilteristunedtolocateasmanyvessels as possible. This may lead to a reduced accuracy in the detection of the extent of bigger vessels. However, as the filter results sym(x) as presented in Eqn. 5 are skeletonized and used as seed points for the subsequent fast marching step, the full extent of the vessels is detected in the second step. Such seeds for front propagation algorithms form starting regions. These starting regions are then extended into the direction of lowest cost until a heuristically stop threshold is reached. 1 f(x)= (7) 1+e−10x The cost matrix, which is only used in the fast marching step, is created by applying 1 f(x)= (8) 1+e−10x as a scaling function to the green color channel intensity values x. The results show enhanced contrast between blood vessels and background and, therefore, are ideal to calculate marching step costs for each possible fast marching step. For the second approach, the phase symmetry filter is applied twice with different parameters in order to find thick as well as very thin blood vessels. For Retinal blood vessel segmentation with logarithmic Gabor Filters 7 Fig.2.Leftcolumn:fullsizeimage,rightcolumn:detailimage.Toprow:greenchannel image, second row: median and tophat filtered image, third row: absolute value of the real part e (x), bottom row: absolute value of imaginary part o (x). n n 8 S. Gross, M. Klein, D. Schneider therecognitionoflargevessels,manydifferenthighscalesarechosenwhereasfor small vessels, smaller and fewer different scales are used. The resulting segmen- tations are combined by hysteresis tresholding to receive the final segmentation. Theoutcomesofthephasesymmetryfiltersandthefinalsegmentationsforboth algorithms are shown in Fig. 3. Fig.3. Top: fast marching approach: left: phase symmetry output, center: thinned image, right: final segmentation, hysteresis approach: left and center: phase symmetry output for different parameter settings, right: final segmentation. 4 Material The performance of the segmentation algorithms is evaluated on two publicly available databases of retinal fundus images, which have been published for the purpose of comparing approaches on vessel segmentation with identical data. The DRIVE database was published in 2004 by Niemeijer [26] and consists of a test and a training set both containing 20 images. For each image, the field of view (FOV) is specified by a mask. The STARE database was presented by Hoover [12] in 2000 and contains 20 images, 10 of which are pathological. As no predefined FOV is available, it is determined by thresholding prior to the segmentation. Furthermore, STARE doesnotcontainaseparatetrainingset.Ifnotmentionedotherwise,parameters for STARE were determined using the DRIVE training set. Retinal blood vessel segmentation with logarithmic Gabor Filters 9 Forbothdatabasesmanualsegmentationsdevelopedbytrainedobserversare used as gold standard for the evaluation of the segmentation. For the STARE database and the test set of the DRIVE database a second manual segmenta- tion was included by the authors to evaluate inter-observer differences. Intra- observer results are listed in Tab. 1 for reasons of comparison as well. The two algorithms presented in this paper are applied to the DRIVE database [26]. The superior approach is also tested on the STARE [12] database as it offers additional challenges because it contains pathological images. The segmentationsare comparedtothefirstmanualsegmentationsusedas gold standard. To obtain the highest comparability, the evaluation criteria used in previous publications [6–12] are adopted. 5 Results For both algorithms, two parameters T and T may be chosen: in case of the 1 2 fastmarchingapproachthesegmentationthresholdforthephasesymmetryfilter and the stop threshold for the fast marching algorithm, in case of the hysteresis approach the thresholds for the two phase symmetry filters. By modifying T 1 andT ,acompromisebetweencorrectlyclassifiedbloodvesselpixelsandfalsely 2 selected background pixels has to be found. The results of these modifications made during extensive testing can be seen in the curves in Fig. 4. For the evaluation on the DRIVE database, parameters T1 and T2 are ad- justed for the 20 images from the test data set and are subsequently used on test data set containing 20 images as well. The resulting average accuracy for the fast marching approach is 94.74% (sens = 73.07%, spec = 97.99%), for the hysteresis approach 94.92% (sens=71.22%, spec=98.41%). In table 1, the re- sultsofformerpublicationsarelistedalongsideournewresults.Fig.4showsthe resultsofthehysteresisapproachusedontheDRIVEdatabasetestsetforvary- ingparametersT andT .TheaverageaccuracyandtheFalsePositiveRateare 1 2 shownasfunctionoftheTruePositiveRate.Thehighestvalueforthemaximum average accuracy of 94.95% is obtained for T =0.18 and T =0.62. 1 2 Thehysteresisthresholdingapproachshowsthesuperiorclassificationperfor- mance for the DRIVE database and a higher sensitivity. Furthermore, observers acknowledged a better representation of smaller blood vessels in the hystere- sis thresholding results. Therefore, this approach is used in additional tests on the STARE database [12] to evaluate the difference between normal and patho- logical images. The maximum average accuracy for healthy images is 95.81% (sens=74.74%,spec=98.35%).Forpathologicalimages,itdecreasesto95.50% (sens=69.16%,spec=98.28%).Thebestresultfortheoverallaverageaccuracy was 95.65% (sens=71.87%, spec=98.34%). 6 Conclusion and outlook ThepaperevaluatesphasesymmetryfilteringasintroducedbyKovesi[21,22]for the segmentation of vascular patterns on retinal images. The main component 10 S. Gross, M. Klein, D. Schneider Author DRIVE STARE MAA Sens Spec MAA Sens Spec PS+FM 0.9474 0.7307 0.9799 - - - PS+Hysteresis 0.9492 0.7122 0.9841 0.9565 0.7187 0.9834 Budai [7] 0.9490 0.7509 0.9680 0.9380 0.6510 0.9750 Wu [8] - 0.9023 0.9518 - - - Soares [16] 0.9466 - - 0.9480 - - Staal [6] 0.9441 0.7194 0.9773 0.9516 - - Ricci [14] 0.9595 - - 0.9646 - - Mendoca [27] 0.9463 0.7315 0.9781 0.9479 0.7123 0.9758 Niemeijer [26] 0.9416 0.6898 0.9696 - - - Fadzil [28] 0.9316 0.9084 0.9335 - - - Rezatofighi [13] 0.9462 0.7358 0.9767 - - - Garg [29] 0.9361 - - - - - Al-Diri [30] - 0.7282 0.9551 - 0.7521 0.9681 Alonso-Montes [31] 0.9180 - - - - - MacGillivray [32] 0.9400 - - - - - Farzin [33] 0.9370 - - 0.9480 - - Fraz [34] 0.9352 - - 0.9384 - - Hoover [12] - - - 0.9275 - - Lam [35] - - - 0.9474 - - intra-observer [13,26] 0.9473 0.7761 0.9725 0.9349 0.8951 0.9572 only background [26] 0.8727 0.0000 1.0000 0.8958 0.0000 1.0000 Table 1. Comparision of maximal average accuracy (MAA), sensitivity (sens), and specificity (spec) taken from different publications for DRIVE [26] and STARE [12].