Image Anal Stereol 2014;33:189-200 doi:10.5566/ias.1101 Original Research Paper AUTOMATIC DETECTION AND CLASSIFICATION OF RETINAL VASCULAR LANDMARKS HADI HAMADC,1,2, DOMENICO TEGOLO2,3 AND CESARE VALENTI2 1Department of Mathematics, An-Najah National University, Palestine; 2Dipartimento di Matematica e Informatica, Universit` a degli Studi di Palermo, Italy; 3Centro Interdipartimentale di Tecnologie della Conoscenza, Universit`a degli Studi di Palermo, Italy e-mail: hadihamad@najah.edu, domenico.tegolo@unipa.it, cesare.valenti@unipa.it (Received December 13, 2013; revised April 9, 2014; accepted May 7, 2014) ABSTRACT The main contribution of this paper is introducing a method to distinguish between different landmarks of the retina: bifurcations and crossings. The methodology may help in differentiating between arteries and veins and is useful in identifying diseases and other special pathologies, too. The method does not need any special skills, thus it can be assimilated to an automatic way for pinpointing landmarks; moreover it gives good responses for very small vessels. A skeletonized representation, taken out from the segmented binary image (obtained through a preprocessing step), is used to identify pixels with three or more neighbors. Then, the junction points are classi.ed into bifurcations or crossovers depending on their geometrical and topological properties such as width, direction and connectivity of the surrounding segments. The proposed approach is applied to the public-domain DRIVE and STARE datasets and compared with the state-of-the-art methods using proper validation parameters. The method was successful in identifying the majority of the landmarks; the average correctly identi.ed bifurcations in both DRIVE and STARE datasets for the recall and precision values are: 95.4% and 87.1% respectively; also for the crossovers, the recall and precision values are: 87.6% and 90.5% respectively; thus outperforming other studies. Keywords: retinal vessel landmark points, retinal vessel structure classi.cation. INTRODUCTION this is due to the deep and comprehensive knowledge she/he has. What drives many researchers to focus on studying Vascular landmarks or junction points in the retina and extracting diverse features in retinal images is are divided into two main types: their vital role in understanding, diagnosing, and predicting not only different pathologies of the eye, – branching junctions which divide the vessel into but also estimations of a variety of more general two sub-vessels; diseases in the body. Pathologies of the eye make – crossing junctions, where two separated vessels changes in the properties and structure of the vascular overcross each other. tree, while diabetic retinopathy usually leads to the formation of new vessels and branches, hypertension Branching junctions might be additionally divided into two types: Y-junction or bifurcation junction in reduces the caliper of arteries, and vessel occlusion which a vessel splits into two vessels of comparable extends their length (Fathi et al., 2013). Moreover, measures, and T-junction where a minor (of smaller the identi.cation of vascular landmarks can help width) vessel protrudes from a major vessel. a lot in pinpointing network connectivity and in constructing retinal vascular trees with anatomical In our paper we use the terms branching and reality, which would be used in many applications such bifurcations interchangeably since for the purpose of as personal identi.cation, mosaicing, and biometric the study there is no need to distinguish between them. security features applications. The landmarks of the retina, especially the For a correct understanding of a vascular tree, crossovers, have a huge number of shapes and the landmarks (bifurcations and crossings) have to structures which makes it dif.cult to construct a be identi.ed. These vascular landmarks are milestone direct technique to locate them. Spherical shape of features which can be used in personal authentication the eye, inconsistent re.ectance of light in the vessels, systems, and also in image registration processes to the existence of pathologies, a big variety of vessel follow up a patient’s status (Fathi et al., 2013). The widths and directions do not let the fundus image be perceptual recognition of the landmarks is much easier good enough to obtain the perfect segmentation which for a physician compared with technical approaches; would be used to identify these landmarks. In the last decade two main techniques were designed to pinpoint the vascular landmarks: geometrical-features based (adopted by this study) and model based approaches (Azzopardi et al., 2011). The geometrical-features based technique depends on image processing manipulation including segmentation and thinning followed by pixel processing and junction-point analysis. This is a robust technique for localizing bifurcations but needs extensive pixel processing. The model based approach is more adaptive and needs less computation, but it suffers from being less generalized since not all features can be modeled easily (Azzopardi et al., 2011). Many of the important features of vessels, such as width and angles between them, have to be measured in the process of classifying landmarks; these features provide other diagnosing parameters which can be calculated directly, e.g., the arterial/venous width ratio (AV-Ratio) which in certain cases is an independent risk factor for stroke, as well as some eye diseases (Xu et al., 2011). In the classi.cation process, on the contrary of our method, many techniques do not take these features into consideration, which may result in increasing the false positive validation parameter (FP). RELATED WORK This .eld of research involved a number of academics who worked on retinal segmentation, but few of them paid enough attention to identifying the landmarks of the vascular tree of the retina. Al-diri et al. (2010) presented an algorithm in which a retinal vessel graph is formed by analyzing the potential connectivity of segmented retinal vessels. Self-organizing feature maps have been adopted to model implicit cost functions for the junction geometry, and the network connectivity is identi.ed by resolving the con.guration of local sets of segment ends. Another paper done by Calvo et al. (2011) divided the method into two successive steps. The former, applies imaging techniques (.lters and morphologic operations) to obtain the base structure for vessels detection. The latter, classi.es (crossover or bifurcation) by analyzing the feature point’s environment. Fathi et al. (2013) suggested also a proper local vessel pattern operator, applied to the centerline of the already thinned segmentation. This operator has a circular structure to extract and classify vessel-point features. Nayebifar et al. (2013) used an approach based on particle .ltering to specify and track the vessel paths in retina. After a median .ltering stage is applied, the product of the green and blue channels, in the RGB retinal image, is used as the probability density function in tracking and describing blood vessels. Initiating the starting point set around the optic disc, the propagating test decides whether a pixel is inside or outside a vessel, and then clusters the pixels by thresholding. The tracking proceeds towards a bifurcation or ends of the vessels. In their paper, Lin et al. (2012) intended to identify the vascular tree by vessel segmentation using the Kalman .lter depending on the continuities in curvature, width and intensity changes at the bifurcation and the crossing points. Tracing errors are resolved by applying a minimum-cost matching algorithm at junctions. Their system was trained on 20 images of the DRIVE dataset, and tested on the remaining 20 images. Another method given by Azzopardi et al. (2011) was used to identify the landmarks by an initial training phase to con.gure a bank of Gabor .lters; yet, that method did not proceed to distinguish between bifurcations and crossovers. Aibinu et al. (2010) used a new hybrid approach called the combined cross-point number method to detect the bifurcation and intersection points in the fundus images using a window of 5×5 pixels. Dashtbozorg et al. (2014) proposed an automatic graph-based approach to classify the vessels into arteries and veins. In order to identify each graph node as being a bifurcation or a crossing, the following steps were applied: the centerlines are extracted from the segmented image, inferring a graph representation of the centerlines, and different modi.cations and enhancements are applied on the junction points in the graph in order to classify the landmarks. We already developed methodologies to automatically segment retinal images by using AdaBoost (Lupas¸cu et al., 2010), self-organizing maps and K-means (Lupas¸cu et al., 2011a), Fuzzy C-Means clustering (Lupas¸cu et al., 2011b). All these approaches were validated by using the datasets Digital Retinal Images for Vessel Extraction (DRIVE) (Staal et al., 2004). The method proposed here introduces a new technique which uses the graph representation of the segmented images and added extra features round the junction points to achieve better reliable classi.cation performances. Both DRIVE and STARE were chosen in this work since they are popular publicly available datasets and a lot of previous and current works are using them (Aibinu et al., 2010; Azzopardi et al., 2011; Fathi et al., 2013), thus more comprehensive comparison between our method and those of other researchers will be achieved. Image Anal Stereol 2014;33:189-200 MATERIAL The 40 DRIVE photographs are randomly selected ones out of 400 for diabetic subjects aged between 25 and 90 years by using a diabetic retinopathy screening program in The Netherlands; only 7 out of 40 show signs of mild early diabetic retinopathy. Each image was captured using 8 bits per color plane at 768 by 584 cropped into 565 by 584 pixels saved in a JPG compressed format. The set is divided into training and test sets, where the .rst part contains 20 images with their corresponding manual binary segmentation, while the latter part (test set) which contains also 20 images, has 2 independent manual binary segmentations, the 1st is used as the ground truth while the 2nd is used in general for comparative evaluation. Observers that manually segmented the vasculature were instructed and trained by experienced ophthalmologists. They were asked to mark as belonging to vessels those pixels of which they were at least 70% certain. Detailed information can be seen later on (Niemeijer et al., 2004). Regarding the STARE dataset (Hoover et al., 2000), it also consists of twenty colored retinal images; ten of them have pathological features. The size of each image is 700 by 605 pixels, and 8 bits per color channel and segmentation of these images was done by two observers (Fraz et al., 2012). METHODOLOGY Two main phases are presented in the method of this study: preprocessing and classi.cation. During the preprocessing stage, the binary segmented image is analyzed specifying the segments, junctions, widths, directions... While during the second stage, each junction is classi.ed into a bifurcation or a crossover. PREPROCESSING STAGE A colored fundus image, with its binary segmentation, is the starting point for the preprocessing stage. Thinning of the segmented images (Lam et al., 1992) was adopted in this contribution instead of skeletonization (Whan et al., 1993), owing to the fact that the thinned image presents smoother segments and contains less spikes which may give a false impression of branches and cusps. Different examples are shown in Fig. 1. Thinning is a binary morphological iterative operation which preserves the original topological information of the structures and returns a thin representation of the centerline of a given vessel by assuming eight-connectivity con.guration (Bhuiyan et al., 2007). Fig. 1. Sample segmentations (a); skeleton (b); thinning (c). Main differences are highlighted. Due to the low resolution of the segmented images, one or two pixel gaps confuse the thinning operation, since that will create an edge connecting the two sides while this edge might not be existing in the actual segmented image. Therefore, our methodology enhances the segmentation automatically by applying a set of three rules, in the following order: – should two components with less than 200 black pixels be weakly-connected one to another, then this connection must be reinforced to become a strong one by introducing a pixel as belonging to background (see Fig. 2a – red arrow); – should a black pixel be weakly connected to a component with more than 1000 black pixels, then the former pixel must be labeled as belonging to a vessel (see Fig. 2a – green arrow); – all strongly connected components with no more than two black pixels must be assigned to the surrounding vessel (see Fig. 2a – blue arrow). Fig. 2. Original (a) and modi.ed (b) segmentations with the corresponding thinning representations. The meaning of the arrow color is described in the text. A portion of a segmented image (Fig. 2a) is processed (Fig. 2b) to clarify the effect of modi.cations on their corresponding thinning images Fig. 2 left respectively. The thinned resulting graph will represent the vascular network (Lupas¸cu et al., 2013a), in which links and nodes identify vessels and intersection points (bifurcations and crossings) respectively; each link is connected at most to two nodes, while each node is connected at least to three links. Thus the graph is composed by two sets: the .rst set contains the links and the other one the intersection points, all reserved in image form. In order to identify all those elements, a labeling process is applied; as a result, each element in the image will be identi.ed by a different label. The essential structures of the thinned form are de.ned as follows: a vertex is a pixel in the binary image connected to at least three components. In particular, a bifurcation is a vertex connected to exactly three components while a direct crossing is de.ned as a vertex which has four connected components. When we have two wide segments crossing each other, usually their thinning is represented by two vertices with a common edge in addition to four outer segments (Fig. 3). This type of crossing is called a non-direct crossing. This last structure makes it dif.cult to distinguish between a crossing and two consecutive bifurcations. It was experimentally veri.ed that 25 pixels are suitable as the maximum allowable edge length, for considering the structure as a non-direct crossing. For the sake of simplicity, we will consider both direct and non-direct crossings just as crossings. We need now to de.ne and sort two sets of vertices: those of probable direct crossings and bifurcations; each of these two sets is used in the classi.cation stage. As far as we know, those direct crossings, connected to 4 surrounding segments, are plausible ones and in need of con.rmation, while the non-direct crossings consist of all pairs of adjacent bifurcations such that the distance between the pair of vertices is less than 25 pixels. a b c Fig. 3. A colored crossing (a), its segmented version (b) and the centerline (c) with two vertices (red) and the edge between (yellow). Extracting the real crossings from the probable ones is really challenging owing to the complexity of the vascular structure; variant possible shapes, tortuosity, in addition to distance between the two vertices in a non-direct crossing and the number of segments, had to be taken into consideration while designing the needed technique for specifying the segments connected to the crossing. The proposed technique runs as follows: .rst, proposing the minimum rectangle that includes the vertex/vertices followed by expanding the sides of the rectangle by one pixel in each direction then detecting those involved segments in the proposed crossing. Each segment is labeled with a distinctive number as can be seen in Fig.4. In this way, we specify the segments which are directly connected to this/these vertex/vertices without interference with other non-intended segments. The segments connected to the vertices are counted in order to propose the plausible landmark type and their anticlockwise ordering round the vertex/vertices is determined and reserved for further investigation; since the features of each segment such as caliper, direction, and connectivity are involved in the criteria for de.ning the real crossings. Counting the number of segments round a vertex is a different approach with respect to the circular probe introduced by Calvo et al. (2011), where changed radius is critical in detecting the neighboring segments and consequently in the overall performance. In the example of Fig. 4, bifurcations are represented in blue, a probable direct crossing is highlighted in red and a probable non-direct crossing with an edge between its two successive bifurcations colored in yellow. Fig. 4. Detecting the landmarks: the segmented image (a), labeled thinned representation of the vessels (b) and rectangles including the vertex/vertices (c,d,e) for direct crossing, bifurcation and a non-direct crossing respectively. Measuring Vessel Direction The vessel direction is one of the most important segment features on which many methods depend, in classifying junctions in the vascular tree (Al-diri et al., 2010; Lin et al., 2010; Fathi et al., 2013; Dashtbozorg et al., 2014). Segment length and direction can be Image Anal Stereol 2014;33:189-200 assessed using its thinned form; length is calculated by counting the number of pixels in the corresponding centerline, but more effort is needed to measure the direction. To identify the slope, we used the singular value decomposition (SVD) algorithm which produces a diagonal matrix S and unitary matrices U and V such that X = U × S ×V *, where S has the same dimension of the input matrix X with non-negative diagonal elements in decreasing order, called the singular values of S. The slope angle of the given segment can be calculated from the two matrices V * and S as the .nal output of a special program named SVD perpendicular direction (SVDPD). Since vessels are usually considered smooth with slight curvature (Fang et al., 2005), the direction of a vessel round a junction is approximated by measuring the direction of the nearest portion of the vessel to the junction. As seen in Fig. 5, the direction of each segment round the junction point is measured in the range between 0 and 180 degrees, which does not immediately determine the angles between the segments around the junction points or the difference in direction between the opposite segments. Fig. 5. To compute the directions, each segment is measured by considering a portion of the centerline round a junction and then appropriate values are added to let the angles range between 0 and 360 degrees. An additive value is estimated to de.ne the segment direction round the vertex as de.ned in the usual Cartesian plane, with the modi.ed measure of the angles lying in the range between 0 and 360 degrees; thinning may have a segment with non-linear pixels so that different quadrants are involved by the segment, especially for the non-direct crossing where there are two vertices. The suggested form to de.ne the additive value uses a Cartesian plane rotated 45 degrees, as shown in Fig. 6, where the yellow region is a buffering zone. Thus, a modi.cation of the measure of the direction of each segment is obtained by adding an additive value depending on the position of the segment relative to the vertex and the measure of the direction being less than or greater than 90 degrees. An illustration of this is given in Fig. 6. Fig. 6. The additive value .ad is de.ned to modify the direction of the segment regarding the vertex. Measuring Vessel Width Several techniques are presented in the literature to measure the width of a vessel in the retina (Niemeijer et al., 2011; Fathi et al., 2013; Lupas¸cu et al., 2013b) each with different level of accuracy, and since an approximated width is suf.cient for the method due to the used resolution, two essential techniques were used both of which rely on the following principle: a perpendicular straight line to the centerline is de.ned through a pixel of the centerline to minimize the error during measurement. The width of the vessels at each pixel in the segment (Niemeijer et al., 2011) is measured by .nding the number of pixels common to the segment and the perpendicular line to it including the pixel of the centerline itself. This process is applied to each pixel in the centerline round the vertices within a9×9 window. Shortness of length of the thinned segments and the presence of several connected vessels to the segment, limit the effectiveness of the calculation process. Hence, two different techniques are used to extract segment’s width: for those segments of length with 6 pixels or more, we used the perpendicular line method; since this will be ef.cient and reliable in this case. But, when the length of the segment is less than 6 pixels a rotational technique is employed. These two techniques are de.ned below. The perpendicular line method. The same method we de.ned in detecting the direction of the centerline is used here; the width is measured through a perpendicular straight line whose direction is de.ned at each pixel and not as detecting an overall direction of the centerline round a vertex by using the segment as a whole. To do so, a gradient box with the size of order two by two pixels is used. In Fig. 7a a perpendicular line extracted from the thinned segment is superimposed on the segmented image to calculate the width which is the Euclidean distance between the two edge pixels intercepted by the edges of the segment and the perpendicular line. The rotational line method. For small segments, a set of twelve rotated straight lines, each spaced 15 degrees from the next one, are positioned around the intended centerline pixel; the number of pixels common to each line together with the centerline pixel is found and the minimum number of the twelve cases is taken to be the width of the segment (See Fig. 7b) (Hatanaka et al., 2011). Fig. 7. Perpendicular (a) and rotational (b) methods to compute the width of the vessel. Both approaches return a measure equal to 3 pixels. CLASSIFICATION STAGE In a crossing, at least one segment of each pair of the two crossing vessels should be connected to some other vessels. Otherwise, this vessel (pair of segments) lacks the mechanism of supply of blood that is characteristic of the vascular tree; this necessitate its being 2 branches of a main vessel. In order to make the decision more reliable, a threshold segment length for such cases was used, such that if at least one of the two non-connected segments is longer than the threshold, then this candidate crossing will be checked through the crossing identi.cation process (Fig. 8). Fig. 8. The two short non-connected segments show an apparent crossing, yet they are classi.ed as 2 branches of a main vessel. Identi.cation of Direct Crossings In order to con.rm or reject each candidate crossing de.ned in the preprocessing stage, we have to take into account further distinctive features such as width and direction of the segments; each two continuing segments generally have nearly equal widths and directions (see Fig. 9). This assumption might be violated in case of a disease or a non-precise segmentation, and moreover the low resolution of the images obliges a permissible threshold parameters. Rules about the width and the direction are de.ned as follows: Width restrictions for direct crossings. Suppose the widths of the .rst pair of segments of a vessel are w1,w3 and the second pair has the widths w2,w4 (see Fig. 9a). Two width conditions are checked for each two segments of the same pair: the width ratio (.) and the width difference (.). For each pair, we compare their width percentages, and then the minimum percentage min{w1,w3} min{w2,w4} is chosen, the function 1 gives its mathematical expression. . = min,. (1) max {w1, w3}max{w2,w4} This ratio should not be less than 0.15 (empirically determined), otherwise, the widths of the two segments of a pair are not considered as the same vessel and hence the assumption of a crossing is not true. The width difference (.) is calculated by the equation: . = max{|w1 - w3|,|w2 - w4|} , (2) which represents the maximum difference between the width of the two segments of each pair, the difference value (empirically determined) was chosen to be less than 6.2 pixels; otherwise, the two segments are not considered to belong to the same vessel; hence the assumption of a crossing is not true. Direction restrictions for direct crossings. To identify the direction of each surrounding segment to a vertex, the centerlines of the segments are used. All directions are reported lying in the range 0–360 degrees (see Fig. 9), with respect to the vertex taken as the center of the crossing. Fig. 9. Direction and width of each segment as . and w respectively (a), angles between the segments (b) and the supplementary angles (c). Given .1 + .2 . ß1 + ß2 (Fig. 9b), the following de.nitions are adopted: .1: The smallest angle of the pair .1 and .2 ß1: The smallest angle of the pair ß1 and ß2. Image Anal Stereol 2014;33:189-200 The following conditions were empirically determined: 0 . .1 . .2 . 173. , 35.. ß1 . ß2 . 180. . .1 and .2: The supplementary angles of the difference in direction between each segment and its continuation (Fig.9c). The values of .1 and .2 have positive measure less than or equal 90., so as to avoid allowing crossings of such high turning vessels. If the above conditions for direction of segments in a candidate direct crossing do not hold then the assumption of a crossing is not true. If a candidate direct crossing satis.es both conditions of width and direction then it is regarded as a real direct crossing. The con.rmed direct crossings are not used any more in further classi.cation procedures, while still candidate direct crossings have to be checked again in next identi.cations. Identi.cation of non-direct crossings Not any two successive adjacent bifurcations will constitute a non-direct crossing; this depends on the distance between their vertices, and the directions of the opposite segments in addition to their width. As in the direct crossings, width and direction are the main features used to classify the non-direct crossings in addition to connectivity, but with more restrictive conditions. The main difference between the images of the two types of crossings is the non-actual edge imposed between the two vertices in the thinned form. The vertices are taken in pairs under the following condition: The edge (d) between two connected vertices must not exceed 25 pixels (Fig. 10); otherwise, they are considered as two bifurcations. All of the candidate non-direct crossings are sorted in ascending order regarding the edge length. Starting from the shortest edge, each pair is checked if it satis.es the non-direct crossing criteria depending on width and direction. Width restrictions for non-direct crossings. In the non-direct crossing, the same two width conditions used in direct crossing are also checked for the two corresponding segments: the width ratio (.) and the width difference (.). The . ratio is de.ned as in Eq. 1, it is not to be less than 0.25. If . does not satisfy this condition, the two segments will not be considered as a continuation of the same vessel. In the same way, as in direct crossing, the width difference . is calculated through Eq. 2, where . should not exceed 3 pixels. Direction restrictions for non-direct crossings. To identify the direction of each surrounding segment to a vertex, approximately the same technique was used as the one in the direct crossing. The following de.nitions are adopted (Fig. 10): – .1: The smaller angle of the two opposite angles connected by the segment “d”. – .2: The greater of the two opposite angles. – ß1 < ß2. The restrictions on these conditions are as follows: 0 . .1 . .2 . 166. , .1 + .2 . 260. , .2-.1 .135. , ß2-ß1 .113. – .1 and .2: The difference in direction between the two segments of each pair in the proposed crossing vessels, where 88.. .1 . .2 . 267. The non-direct crossings candidates satisfying these conditions are considered non-direct crossings, and their vertices will not be involved in any further investigations. Non-standard crossings Although crossings with more than four segments are quite rare in fundus images (Dashtbozorg et al., 2014), the misclassi.cation, resulting from the complexity of such crossings, might destroy the retinal vascular tree reconstruction. Only few researchers tried to identify this kind of crossing in their work (Saez et al., 2012; Dashtbozorg et al., 2014). The proposed approach in this study is limited to the case involving 5 segments constituted of two main vessels crossing each other and a daughter branch from one of them. RESULTS The main performance measures which are used in evaluating the landmark classi.cation, depend on the comparison between the description of the found landmarks and those manually provided by expert and known as ground truth. The following indices have been used in evaluating our methods classi.cation performance: – True Positive (TP) indicates the amount of real feature points detected by the system. – False Negative (FN) is the amount of real feature points not detected by the system. – False Positive (FP) indicates the number spurious feature points detected by the system – Recall re.ects the ability of the algorithms to detect the landmarks, RE = TP/(TP + FN) – Precision is the rate of true positive to all positive results, PR = TP/(TP + FP) – Jaccard Index which measures the similarity between two sets A and B is de.ned by the equation: J(A,B)= |A . B|/|A . B|. The Jaccard index was calculated by using the formula J = TP/(TP + FP + FN). The suggested unsupervised approach was validated on the 20 test images of both the DRIVE and STARE datasets, using their available .rst manual segmentation (Fraz et al., 2012) to identify the structural components of the vascular tree. The optic disk, which can be located by using the Harris detector (Bellavia et al., 2011; Dehghani et al., 2012), was excluded from the evaluation of the method due to the complexity, and multi intertwining of its vessels (Dashtbozorg et al., 2014). This area usually contains many vessels and the corresponding graph is not reliable; detection percentages for both bifurcations and crossings are shown in Table 1. Table 1. The detection percentages of bifurcations and crossings for DRIVE and STARE datasets. DRIVE TP % FN % FP % Bifurcations 94.3 5.7 17.7 Crossings 86.5 13.5 13.1 STARE TP % FN % FP % Bifurcations 95.3 4.7 13.0 Crossings 86.8 13.2 7.6 Moreover, the method took into consideration all vessels, without excluding the smallest ones. The proposed methodology asserts its being more robust and powerful than most of the other methodologies (Bhuiyan et al., 2007; Al-diri et al., 2010; Aibinu et al., 2010; Fathi et al., 2013), which is con.rmed by the results. The method works properly with different datasets resolutions (Table 2); the recall and precision values for both the DRIVE and STARE datasets show comparable results for their mean, standard deviation, minimum, and maximum values. This gives a certi.cation of robustness to the methodology. Table 2. Performance for classi.cation of bifurcations and crossings in DRIVE (DR) and STARE (ST). Bifurcations % Crossings % RE PR RE PR mean DR 94.9 84.5 86.8 87.6 ST 95.9 89.6 88.4 93.3 stdDR 3.45.16.87.7 ST 3.37.48.07.0 min DR 86.8 71.8 76.0 75.9 ST 89.2 77.4 71.4 79.0 max DR 100. 93.8 96.2 100. ST 100. 100. 100. 100. Since the recall and precision values do not show the in.uence of FN and FP on the performance of the method simultaneously, a Jaccard index was introduced and the results are shown in Table 3 where the effect of FP values reduced the performance in comparison to recall and precision values. Table 3. Detection performances percentages of bifurcations and crossings of different studies on the DRIVE and STARE. Best results are highlighted. DRIVE Bifurcations % RE PR CO J RE Crossings % PR CO J Al-Diri et al. 71. - - - 48. 69. 58.5 - Aibinu et al. Fathi et al. Our method 94.8 80.8 94.9 51.9 88.4 84.5 73.3 84.6 89.7 --80.2 4.6 85.9 86.8 85.3 78.9 87.6 44.9 82.4 87.2 --76.5 STARE Bifurcations % RE PR CO J RE Crossings % PR CO J Aibinu et al. 97.2 52.8 75. -6.1 87.5 46.8 - Fathi et al. 88.2 83.8 86. -82.9 85.8 84.3 ­ Our method 95.9 89.6 92.7 84.3 88.4 93.3 90.8 80.7 On the other hand, our methodology is stable in classifying the landmark points, with respect to both Image Anal Stereol 2014;33:189-200 precision and recall: the plots in Fig. 11 show that most classi.cations exibit a value of at least 80% for PR and RE measurements. In order to compare and to evaluate the performance of our methodology, a manual detection of landmark points was carried out by an expert ophthalmologist on 20 test images for both the datasets, where he identi.ed the locations and types of the junction points as bifurcations or crossings. In particular, the DRIVE test images were manually identi.ed to have 2167 bifurcations and 520 crossings, while the algorithm located 2427 bifurcations, 192 direct crossings and 326 non-direct crossings (a total number of 518 crossings). Regarding the STARE dataset, it was also manually identi.ed to have 1146 bifurcations and 370 crossings, while the algorithm located 1241 bifurcations, 97 direct crossings and 252 non-direct crossings (a total number of 349 crossings). points by our methodology compared with those classi.ed by the physician; – Two segments may intertwine more than once (Figs. 12c, d) which is not easy to analyze in the thinned form; – Wrongly identi.ed bifurcations or crossings due to non-veri.ed segments in the binary image (Figs. 12e, f); – Non usual vessel structure (Figs. 12g, h). – Crossings not correctly identi.ed are considered as bifurcations (Lin et al., 2010), this leads to a reduced number of automatically detected crossings with respect to those manually identi.ed ones and to more bifurcations. RE RE RE RE PR PR PR PR Fig. 11. Precision vs recall for the DRIVE (blue) and STARE (red) datasets classi.ed into bifurcations (triangles) and crossings (circles). We observed a difference in the classi.cation ability among different images regarding existing pathologies and the quality of their respective segmentations. In classifying particular junction points, dif.culties may arise due mainly to low resolution segmentation (Al-diri et al., 2010). For instance, here are some of these complications: – Two segments may pass very near to each other but do not really cross, (Figs. 12a, b) which may lead to a higher number of detected bifurcation Image Anal Stereol 2014;33:189-200 The explanation of extra identi.ed false bifurcations with respect to crossings, as shown in Table 1, might be a result of high detection ability even for very small segments, i.e., protrusions, and also to non-suf.cient separation between near segments; again, low resolution in the segmented image limits the ability to have a proper separation for such near segments (Figs. 12a, b). Moreover, the non identi.ed crossings are regarded as bifurcations. The difference between the best and the worst results are evident in Fig. 13. Comparing the recall and precision for both bifurcations and crossings in Table 2, we notice that, bifurcations recall is higher than precision. In the crossing case, precision is generally higher than recall, which indicates that the algorithm detects real crossings in case of standard ones but it is liable to get confused in case of nonstandard crossings. To the best of our knowledge, the method excels the present mentioned results in literature for detecting landmarks of the retina on the used datasets in this study. To evaluate the performance of our methodology regarding other works, a comparison is presented in Table 3 between 3 previous studies and ours. Where correctness (CO) was de.ned as the average of the recall and the precision values. The recall and the precision rates in some previous works, are evaluated with big differences, such as the results in Aibinu et al. and the crossings correctness value did not exceed 58.5% in Al-Diri et al. on the DRIVE dataset. However, comparing the previous two results to Fathi et al. his work got more reasonable performance measurements for both datasets. Our method shows approximately 6% improvement in overall performance in comparison to the best results which were presented by Fathi et al. (2013). DISCUSSION We introduced here a novel non-supervised method to detect landmark points in binary images coming from retinal photographs. These segmented images can be obtained by using a variety of approaches already present in the literature and validated on both DRIVE and STARE datasets which are considered as a standard corpus to make fair comparisons (Niemeijer et al., 2004). The twenty training images of the DRIVE dataset were used to .ne-tune and modify the performance of our method before applying the algorithm to the two testing set images. The methodology on which we built our work was mathematical modeling of the way the physician understands and analyzes the segmented image. It must be stressed that both segmentation and classi.cation of the landmark points in the ground truth are subjective and prone to errors. Our approach showed a high performance in locating and classifying most junction points. The .nal results of the study were compared with state-of­the-art works, showing generally better overall recall and precision of the study in both bifurcations and crossings. The average recall value for bifurcations and crossings detection in the DRIVE and STARE datasets is 91.5%, and the average precision value is 88.8%. In the case of ambiguous junctions, our approach tried to introduce the best expected solution of that particular con.guration. Moreover, further investigation is needed due to the presence of eventual pathologies, the low resolution of the images we used and the subjective segmentations provided by the physicians, which added an extra burden to the work. Building on the study an ef.cient algorithm may be developed for classifying the vasculature of the retina into arteries and veins which is important in identifying diseases; the approach adapted in this paper can also be used in several practical applications such as personal authentication, registration, etc. ACKNOWLEDGMENT The authors would like to thank Dr. Luigi Di Rosa from the Department of Ophthalmology at the University of Palermo in Italy for his precious time either in identifying the landmarks and also for his valuable notes and recommendations. REFERENCES Aibinu AM, Iqbal MI, Sha.e AA, Salami MJE, Nilsson M (2010). Vascular intersection detection in retina fundus images using a new hybrid approach. Comput Biol Med 40(1):81–9. Al-Diri B, Hunter A, Steel D, Habib M (2010). Automated analysis of retinal vascular network connectivity. Comput Med Imaging Grap 34(6):462–70. Azzopardi G, Petkov N (2011). Detection of retinal vascular bifurcations by trainable V4-like .lters. In: Real P, Dias-Pernil D, Molina-Abril H, Berciano A, Kropatsch W, eds. Proc 14th Int Conf CAIP. Seville, Spain. August 29–31. Lect Not Comput Sci 6854:451–9. Bellavia F, Tegolo D, Valenti C (2011). Improving Harris corner selection strategy. IET Comput Vis 5(2):87–96. Bhuiyan A, Nath B, Chua J, Ramamohanarao K (2007). Automatic detection of vascular bifurcations and crossovers from color retinal fundus images. In: Y´ etongnon K, Chbeir R, Dipanda A, eds. Proc 3rd Int IEEE Conf SITIS. Shanghai, China, December 16–18. 711–8. Calvo D, Ortega M, Penedo MG, Rouco J (2011). Automatic detection and characterisation of retinal vessel tree bifurcations and crossovers in eye fundus images. Comput Meth Prog Biomed 103(1):28–38. Dashtbozorg B, Mendonca AM, Campilho A (2014). An automatic graph-based approach for artery/vein classi.cation in retinal images. IEEE T Image Process 23(3):1073–83. Dehghani A, Moin MS, Sagha. M (2012). Localization of the optic disc center in retinal images based on the Harris corner detector. Biomed Eng Lett 2:198–206. Fang B, You X, Tang YY, Chen WSH (2005). Morphological structure reconstruction of retinal vessels in fundus images. Int J Pattern Recogn 19(7):937–48. Fathi A, Naghsh-Nilchi AR, Mohammadi FA (2013). Automatic vessel network features quanti.cation using local vessel pattern operator. Comput Biol Med 43:587– 93. Fraz MM, Remagnino P, Hoppe A, Uyyanonvara B, Rudnicka AR, Owen CG, Barman SA (2012). Blood vessel segmentation methodologies in retinal images – A survey. Comput Meth Programs Biomed 108:407–33. Hatanaka Y, Muramatsu C, Hara T, Fujita H (2011). Automatic arteriovenous crossing phenomenon detection on retinal fundus images. In: Summers RM, van Ginneken B, eds. Medical Imaging 2011: Computer-Aided Diagnosis. Lake Buena Vista, USA. February 15–17. Proc SPIE 7963:79633V Hoover A, Kouznetsova V, Goldbaum M (2000). Locating blood vessels in retinal images by piecewise threshold probing of a matched .lter response. IEEE T Med Imaging 19(3):203–10. Lam L, Lee SW, Suen CY (1992). Thinning methodologies-a comprehensive survey. IEEE T Pattern Anal 14(9):869–85. Lin KSh, Yeh PY, Tsai ChL, Chen ShJ, Lin WY (2010). Retinal vascular tree construction with multimodal .uorescein angiogram sequence. Biomed End-App Bas C 22(2):101–10. Lin KS, Tsai CL, Tsai CH, Sofka M, Chen SJ, Lin WY (2012). Retinal vascular tree reconstruction with anatomical realism. IEEE T Biomed Eng 59(12):3337– 47. Lupas¸cu CA, Tegolo D, Trucco E (2010). FABC: Retinal vessel segmentation using AdaBoost. IEEE T Inf Technol B 14(5):1267–74. Lupas¸cu CA, Tegolo D (2011a). Automatic unsupervised segmentation of retinal vessels using self-organizing maps and K-means clustering. In: Rizzo R, Lisboa PJG, eds. Proc 7th Int Meet CIBB. Palermo, Italy. September 16–18. Lect Not Comput Sci 6685:263–74. Lupas¸cu CA, Tegolo D (2011b). Stable automatic unsupervised segmentation of retinal vessels using self-organizing maps and a modi.ed Fuzzy C-Means clustering. In: Fanelli AM, Pedrycz W, Petrosino A, eds. Proc 9th Int Worksh WILF. Trani, Italy. August 29–31. Lect Not Comput Sci 6857:244–52. Lupas¸cu CA, Tegolo D, Bellavia F, Valenti C (2013a). Semi-automatic registration of retinal images based on line matching approach. In: Proc 26th IEEE Int Symp CBMS. Porto, Portugal. June 20–22. Proc IEEE 26:453–6. Lupas¸cu CA, Tegolo D, Trucco E (2013b). Accurate estimation of retinal vessel width using bagged decision trees and an extended multiresolution Hermite model. Med Image Anal 17(8):1164–80. Nayebifar B, Moghaddam HA (2013). A novel method for retinal vessel tracking using particle .lters. Comput Biol Med 43(5):541–8. Niemeijer M, Staal J, van Ginneken B, Loog M, Abramoff MD (2004). Comparative study of retinal vessel segmentation methods on a new publicly available database. In: Fitzpatrick JM, Sonka M, eds. Medical Imaging 2004: Image Processing. San Diego, USA. February 16–19. Proc SPIE 5370:648–56. Niemeijer M, Xiayu X, Dumitrescu AV, Gupta P, Ginneken Bvan, Folk JC, Abramoff MD (2011). Automated measurement of the arteriolar-to-venular width ratio in digital color fundus photographs. IEEE T Med Imaging 30(11):1941–50. Saez M, Gonz´azquez S, Gonz´oalez-V´alez-Penedo M, Barcel´MA, Pena-Seijo M, Tuero GCDe, Pose-Reino A (2012). Development of an automated system to classify retinal vessels into arteries and veins. Comput Meth Programs Biomed 108(1):367–76. Staal J, Abramoff MD, Niemeijer M, Viergever MA, Ginneken Bvan (2004). Ridge-based vessel segmentation in color images of the retina. IEEE T Med Imaging 23(4):501–9. http://www.isi.uu.nl/Research/Databases/DRIVE/ Whan LS, Louisa L, Ching YS (1993). A systematic evaluation of skeletonization algorithms. Int J Pattern Recogn 7(5):1203–25. Xu X, Niemeijer M, Song Q, Garvin MK, Reinhardt JM, Abramoff MD (2011). Retinal vessel width measurements based on a graph-theoretic method. Proc 8th IEEE Int Symp Biomed Imaging: From Nano to Macro. Chicago, USA (U.S.A.). March 30–April 2. 641–4.