Analysis of Vessel Connectivities in Retinal Images by Cortically Inspired Spectral Clustering
Marta Favali, Samaneh Abbasi-Sureshjani, Bart ter Haar Romeny, Alessandro Sarti
JJ Math Imaging Vis manuscript No. (will be inserted by the editor)
Analysis of Vessel Connectivities in Retinal Imagesby Cortically Inspired Spectral Clustering
Marta Favali · Samaneh Abbasi-Sureshjani · Bart ter Haar Romeny · Alessandro Sarti
Received: 23 October 2015 / Accepted: 08 February 2016
Abstract
Retinal images provide early signs of dia-betic retinopathy, glaucoma and hypertension. Thesesigns can be investigated based on microaneurysms orsmaller vessels. The diagnostic biomarkers are the changeof vessel widths and angles especially at junctions, whichare investigated using the vessel segmentation or track-ing. Vessel paths may also be interrupted; crossings andbifurcations may be disconnected. This paper addressesa novel contextual method based on the geometry ofthe primary visual cortex (V1) to study these difficul-ties. We have analysed the specific problems at junc-tions with a connectivity kernel obtained as the funda-mental solution of the Fokker-Planck equation, whichis usually used to represent the geometrical structureof multi-orientation cortical connectivity. By using thespectral clustering on a large local affinity matrix con-structed by both the connectivity kernel and the featureof intensity, the vessels are identified successfully in ahierarchical topology each representing an individualperceptual unit.
M. Favali and S. Abbasi-Sureshjani contributed equally tothis work.M. Favali · A.SartiCentre d’Analyse et de Math´ematique Sociales, Ecole desHautes ´Etudes en Sciences Sociales, 190-198, Avenue deFrance 75244 Paris Cedex 13 FranceE-mail: { marta.favali, alessandro.sarti } @ehess.frS. Abbasi-SureshjaniDepartment of Biomedical Engineering, Eindhoven Univer-sity of Technology, P.O. Box 513, 5600 MB Eindhoven, TheNetherlandsE-mail: [email protected]. ter Haar RomenyDepartment of Biomedical and Information Engineering,Northeastern University, 500 Zhihui Street, Shenyang, 110167ChinaE-mail: [email protected] Keywords
Retinal image analysis · Fokker-Planckequation · Cortical connectivity · Spectral clustering · Gestalt
Clinical importance of retinal blood vessels:
Epidemicgrowth of systemic, cardiovascular and ophthalmologicdiseases such as diabetes, hypertension, glaucoma, andarteriosclerosis [38, 48, 67], their high impact on thequality of life, and the substantial need for increase inhealth care resources [43, 68] indicate the importance ofconducting large screening programs for early diagnosisand treatment of such diseases. This is impossible with-out using automated computer-aided systems becauseof the large population involved.The retina is the only location where blood vesselscan be directly monitored non-invasively in vivo [77].Retinal vessels are connected and form a treelike struc-ture. The local orientation and intensity of vesselschange gradually along their lengths; however, theselocal properties may vary highly for different vessels.Studies show that the quantitative measurement ofmorphological and geometrical attributes of retinal vas-culature, such as vessel diameter, tortuosity, arteriove-nous ratio and branching pattern and angles are veryinformative in early diagnosis and prognosis of severaldiseases [11, 28, 35, 40, 46, 65, 71]. More specifically,since two blood vessels with the same type never crosseach other and arteries are more likely to cross overveins, studying the behaviour of blood vessels at cross-ings for detecting the branch retinal vein occlusion (oc-clusion of the vein) and arteriovenous nicking helpsin diagnosis of hypertension (increased arterial bloodpressure), arteriosclerosis (hardening of arteries) and a r X i v : . [ c s . C V ] M a y Marta Favali et al. stroke [32, 75, 76]. The other clinically highly promis-ing but still underrated information is based on study-ing the smaller vessels, because it is expected that thesigns of diseases such as diabetic retinopathy appearsin smaller vessels earlier than in larger ones.
Vessel extraction and its difficulties:
The vasculaturecan be extracted by means of either pixel classificationor vessel tracking. Several segmentation and trackingmethods have been proposed in the literature [9, 26,31]. In pixel classification approaches image pixels arelabeled either as vessel or non-vessel pixels. Therefore,a vessel likelihood (soft segmentation) or binary map(hard segmentation) is created for the retinal image.Although the vessel locations are estimated in theseapproaches, they do not provide any information aboutvessel connectivities. On the contrary, in tracking basedapproaches, several seed points are selected and the bestconnecting paths between them are found [2, 10, 12, 16–18, 34, 57, 58, 78]. The main benefit of vessel trackingapproaches is that they work at the level of a singlevessel rather than a single pixel and they try to find thebest path that matches the vessel profile. Therefore, theinformation extracted from each vessel segment (e.g.diameter and tortuosity) is more accurate and reliable.There are several difficulties for both vessel segmen-tation and tracking approaches. Depending on imag-ing technology and conditions, these images could beaffected by noise in several degrees. Moreover, non-uniform luminosity, drift in image intensity, low con-trast regions and also central vessel reflex make thevessel detection and tracking complicated. Several im-age enhancement, normalization and denoising tech-niques have been developed to tackle these complica-tions (e.g. [1, 29, 50]).The tracking methods are often performed exploit-ing the skeleton of the segmented images. Thus, non-perfect segmentation or wrong skeleton extraction re-sults in topological tracing errors e.g. disconnectionsand non-complete subtrees as discussed in several meth-ods proposed in the literature [2, 16, 17, 34, 44]. Typicalnon-perfections include missing small vessels, wronglymerged parallel vessels, disconnected or broken up ves-sel segments and the presence of spur branches in thin-ning. Moreover, the greater difficulty arises at junc-tions and crossovers: small arteriovenous crossing an-gles, complex junctions when several junctions are closetogether, or presence of a bifurcation next to a cross-ing makes the centerline extraction and tracing chal-lenging. These difficulties are mentioned as the trackinglimitations in the literature. Some of these challengingcases are depicted in Fig. 1 with their correspondingartery/vein ground truth labels. Arteries and veins are annotated in red and blue colors respectively. The greencolor represents the crossing and the types of the whitevessels are not known.
Gestalt theory and cortically-inspired spectral cluster-ing:
Visual tasks like image segmentation and groupingcan be explained with the theory of the Berliner Gestaltpsychology, that proposed local and global laws to de-scribe the properties of a visual stimulus [70, 73]. Inparticular, the laws of good continuation, closure andproximity have a central role in the individuation ofperceptual units in the visual space, see Fig. 2. In [54] (a) (b) (c)
Fig. 2: Examples of (a) good continuation, (b) closureand (c) proximity Gestalt laws.perceptual grouping was considered to study the prob-lem of finding curves in biomedical images. In order tostudy the property of good continuation, Field, Hayesand Hess introduced in [27] the concept of an associ-ation field, that defines which properties the elementsof the stimuli should have to be associated to the sameperceptual unit, such as co-linearity and co-circularity.In [7] Bosking showed how the rules of association fieldsare implemented in the primary visual cortex (V1),where neurons with similar orientation are connectedwith long-range horizontal connectivity. A geometricmodel of the association fields based on the functionalorganization of V1 has been proposed in [13]. This ge-ometric approach is part of the research line proposedby [37, 45, 49, 56, 63, 80] and applications to imageprocessing can be found in [6, 23, 24].In this work, a novel mathematical model based onthis geometry has been applied to the analysis of retinalimages to overcome the above mentioned connectivityproblems in vessel tracking. The proposed method rep-resents an engineering application of segmenting andrepresenting blood vessels inspired by the modeling ofthe visual cortex. This shows how these models can beapplied to the analysis of medical images and how thesetwo fields can be reciprocally used to better understandand reinforce each other. ortically Inspired Connectivity Analysis in Retinal Images 3
C1 C2 C5 C6C3 C4 C1 C2C3 C4 C5 C6C2C1 C3 C4 C5 C6
Fig. 1: A sample image selected from the DRIVE dataset [66] with centrelines extracted by applying the morpho-logical skeletonization on the segmented image (top left), and its corresponding artery/vein ground truth fromthe RITE dataset [39] (top right). Several difficult cases are highlighted in this figure. C1: complex junction withthe presence of a bifurcation and a crossing with a narrow crossing angle; C2: interrupted lines and missing smallvessels; C3: high curvature vessel; C4: complex junction, a crossing and bifurcation, wrong thinning at crossing; C5:two nearby parallel vessels merged as one group; C6: missing small vessel, merged parallel vessels and interruptedsegment.This method, which is not dependent on centrelineextraction, is based on the fact that in arteriovenouscrossings there is a continuity in orientation and in-tensity of the artery and vein respectively, i.e. the lo-cal variation of orientation and intensity of individualvessels is very low. The proposed method models theconnectivity as the fundamental solution of the Fokker-Planck equation, which matches the statistical distri-bution of edge co-occurrence in natural images and is agood model of the cortical connectivity [60].Starting from this connectivity kernel and consider-ing the Euclidean distance between intensities of bloodvessels, we build the normalized affinity matrix. Sinceat crossings the vessels have different intensities (types),including this feature in the construction of the affin-ity matrix adds more discriminative information. Thisspectral approach, first proposed for image processing,is inspired by [15, 55, 64, 72]. Moreover, recent resultsof [62] show how the spectral analysis could actually beimplemented by the neural population dynamics of theprimary visual cortex (V1). Finally, we use a spectralclustering algorithm to find and group the eigenvectorslinked to the highest eigenvalues of the affinity matrix.We will describe how these groups represent different perceptual units (vessels) in retinal images. Originally,the spectral clustering was exploited for good contin-uation, closure and proximity, see Fig. 2. It excelledin finding connections, e.g. broken vessel segments. Inthis paper, we go further by solving many challenges atvessel crossings and bifurcations.
Paper structure:
The remainder of the article is orga-nized as follows. In Sect. 2, after describing the geome-try of the functional architecture of the primary visualcortex, it is explained in detail how to lift the stimulusin cortical space using a specific wavelet transform. Thislifting converts the image from the space of positions( R ) to the joint space of positions and orientations( R × S ). Then the connectivity kernel and the con-struction of the affinity matrix based on this kernel andthe intensity is described. Next step is the spectral clus-tering algorithm, that is used to extract the groupinginformation of the stimuli, explained in Sect. 3. The ex-periments applied on retinal images are explained stepby step and the results are presented in Sect. 4. Finally,the paper is concluded in Sect. 5 by briefly summariz-ing the proposed method, discussing its strengths inpreserving the connectivities in the retinal vasculature Marta Favali et al. network and proposing potential improvements as fu-ture work. θ . In other words, simplecells extract the orientation information at all locationsand send a multi-orientation field to higher levels in thebrain. Also, it is well known that objects with differentorientations can be identified by the brain even whenthey are partly occluded, noisy or interrupted [41].Motivated by these findings, a new transformationwas proposed [21, 22],to lift all elongated structures in2D images to a new space of positions and orientations( R × S ) using elongated and oriented wavelets. Bylifting the stimulus, multiple orientations per positioncould be detected. Thus, crossing and bifurcating linesare disentangled into separate layers corresponding totheir orientations.Constructing this higher dimension structure sim-plifies the higher level analysis. However, such repre-sentation is constrained and the invertibility of thetransformation needs to be guaranteed using the rightwavelet. Cake wavelets as a class of proper wavelets [21,22, 33] satisfy this constraint and avoid loss of informa-tion. Cake wavelets are directional wavelets similar toGabor wavelets and have a high response on orientedand elongated structures. Moreover, similar to Gaborwavelets, they have the quadratic property; so the realpart contains information about the symmetric struc-tures, e.g. ridges, and the imaginary part contains infor-mation about the antisymmetric structures, e.g. edges.Although blood vessels can have several scales, by us-ing the cake wavelets multi-scale analysis is not needed,because they capture the information at all scales. Avisual comparison between Gabor and cake wavelets ispresented in Fig. 3. As seen in this figure, by using thecake wavelet in all orientations, the entire frequency do-main is covered; while the Gabor wavelets, dependingon their scale cover a limited portion of Fourier space.Therefore, they are scale dependent. The reader is re-ferred to [5] for more detail.For lifting the stimulus and constructing the higherdimension representation, U f : SE (2) → C , the inputimage f ( x, y ) is correlated with the anisotropic cake wavelet ψ [21, 22, 30]: U f ( x , θ ) = ( W ψ [ f ])( x , θ ) = ( R θ ( ψ ) (cid:63) f )( x )= (cid:90) R ψ ( R − θ ( y − x )) f ( y )d y (1)where R θ is the 2D counter-clockwise rotation matrixand the overline denotes the complex conjugate. Usingthis transformation and considering only one orienta-tion per position (the orientation with highest trans-formation response), the points of a curve γ = ( x, y )are lifted to new cortical curves and are described inthe space ( x, y, θ ):( x, y ) → ( x, y, θ ) . These curves have been modelled in [13] as integralcurves of suitable vector fields: X = (cos θ, sin θ, , X = (0 , , . (2)The points of the lifted curves are connected by integralcurves of these two vector fields such that: γ : R → SE (2) , γ ( s ) = ( x ( s ) , y ( s ) , θ ( s )) (3) γ (cid:48) ( s ) = ( k ( s ) X + k ( s ) X )( γ ( s )) , γ (0) = 0 . These curves projected on the 2D cortical plane repre-sent a good model of the association fields, as describedin [13] (see Fig. 4).In order to include the intensity term, we use theEuclidean distance between the intensities of two cor-responding points. If f ( x, y ) represents the image in-tensity at position ( x, y ) the stimulus is lifted to theextended 4-dimensional feature space:( x, y, θ ) → ( x, y, θ, f ( x, y )) . An admissible curve in this space is defined as the so-lution of the following differential equation: γ (cid:48) ( s ) = ( k ( s ) X + k ( s ) X + k ( s ) X )( γ ( s )) (4) γ (0) = ( x , y , θ , f ) , γ (1) = ( x , y , θ , f )where the vector fields are: X = (cos θ, sin θ, , , X = (0 , , , , (5) X = (0 , , , k and k represent a distance in the( x, y, θ ) domain and k is a Euclidean distance. Start-ing from these vector fields we can model the corticalconnectivity. ortically Inspired Connectivity Analysis in Retinal Images 5(a) (b) (c) (d) Fig. 3: Visual comparison between cake and Gabor wavelets (top and bottom rows respectively). (a) and (b) realand imaginary parts of the wavelets in spatial domain, (c) the wavelet in Fourier domain, and (d) the illustrationof Fourier space coverage by wavelets in 36 orientations. (a) Integral curves (b) Distribution of integral curves
Fig. 4: a) The integral curves of the vector fields X and X in the (x,y, θ ) space. In blue the projections of theintegral curves on the xy plane. b) The distribution of the integral curves, modeling the connectivity betweenpoints.2.2 The Connectivity KernelsThe cortical connectivity can be modelled as the prob-ability of connecting two points in the cortex and isrepresented by the stochastic counterpart of the curvesin Eq. 3:( x (cid:48) , y (cid:48) , θ (cid:48) ) = X + N (0 , σ ) X (6)where N (0 , σ ) is a normally distributed variable withzero mean and variance equal to σ .This process, first described in [49], is discussed in[3, 4, 61, 74]. We denote v the probability density to finda particle at the point ( x, y ) considering that it startedfrom a given location ( x (cid:48) , y (cid:48) ) and that it is moving withsome known velocity. This probability density satisfiesa deterministic equation known in literature as the Kol- mogorov forward equation or Fokker Planck equation: ∂ t v = X v + σ X v where X is the directional derivative cos( θ ) ∂ x +sin( θ ) ∂ y and X = ∂ θθ is the second order derivative.This equation has been largely used in differentfields [3, 4, 23, 74]. In [61] its stationary counterpartwas proposed to model the probability of co-occurrenceof contours in natural images. The Fokker Planck op-erator has a nonnegative fundamental solution Γ thatsatisfies: X Γ (( x, y, θ ) , ( x (cid:48) , y (cid:48) , θ (cid:48) ))+ σ X Γ (( x, y, θ ) , ( x (cid:48) , y (cid:48) , θ (cid:48) ))= δ ( x, y, θ )which is not symmetric. The connectivity kernel ω ob-tained by symmetrization of the Fokker Planck funda- Marta Favali et al. mental solution is: ω (( x, y, θ ) , ( x (cid:48) , y (cid:48) , θ (cid:48) ))= 12 ( Γ (( x, y, θ ) , ( x (cid:48) , y (cid:48) , θ (cid:48) )) + Γ (( x (cid:48) , y (cid:48) , θ (cid:48) ) , ( x, y, θ )) . In order to measure the distances between intensitieswe introduce the kernel ω (( x i , y i ) , ( x j , y j )). This newintensity kernel is obtained as: ω ( f i , f j ) = e ( − ( fi − fjσ )) (7)The final connectivity kernel can be written as theproduct (as these are probabilities) of the two compo-nents: ω f (( x i , y i , θ i , f i ) , ( x j , y j , θ j , f j ))= ω (( x i , y i , θ i ) , ( x j , y j , θ j )) ω ( f i , f j ) . (8)2.3 Affinity MatrixStarting from the connectivity kernel defined previ-ously, it is possible to extract perceptual units fromimages by means of spectral analysis of suitable affinitymatrices. The eigenvectors with the highest eigenval-ues are linked to the most salient objects in the scene[55]. The connectivity is represented by a real symmet-ric matrix A i,j : A i,j = ω f (( x i , y i , θ i , f i ) , ( x j , y j , θ j , f j )) (9)that contains the connectivity information between allthe lifted points. The eigenvectors of the affinity matrixare interpreted as perceptual units [62]. The goal of clustering is to divide the data points intoseveral groups such that points in the same group aresimilar and points in different groups are dissimilar toeach other. The cognitive task of visual grouping canbe considered as a form of clustering, with which it ispossible to separate points in different groups accordingto their similarities. In order to perform visual group-ing, we will use the spectral clustering algorithm. Tra-ditional clustering algorithms, such as K-means, are notable to resolve this problem [51]. In recent years, dif-ferent techniques have been presented to overcome theperformance of the traditional algorithms, in particularspectral analysis techniques. It is widely known thatthese techniques can be used for data partitioning andimage segmentation [47, 55, 64, 72] and they outper-form the traditional approaches. Above that, they aresimple to implement and can be solved efficiently bystandard linear algebra methods [69]. In the next sec-tion we will describe the spectral clustering algorithmused in the numerical simulations of this paper. 3.1 Spectral Clustering TechniqueDifferent algorithms based on the theory of graphs havebeen proposed to perform clustering. In [55] it has beenshown how the edge weights { a ij } i,j =1 ,...n of a weightedgraph describe an affinity matrix A . This matrix con-tains information about the correct segmentation andwill identify perceptual units in the scene, where thesalient objects will correspond to the eigenvectors withthe highest eigenvalues. Even though it works success-fully in many examples, in [72] it has been demon-strated that this algorithm also can lead to clusteringerrors. In [69] and [47] the algorithm is improved con-sidering the normalized affinity matrix. In particularwe will use the normalization described in [47]. Defin-ing the diagonal matrix D as formed by the sum of theedge weights (representing the degrees of the nodes, d i = (cid:80) nj =1 a ij ), the normalized affinity matrix is ob-tained as: P = D − A (10)This stochastic matrix P represents the transition prob-ability of a random walk in a graph. It has real eigen-values { λ j } j =1 ,...n where 0 ≤ λ j ≤
1, and its eigenvec-tors { u i } i =1 ,...K , related to the K largest eigenvalues λ ≥ λ ≥ ... ≥ λ K , represent a solution of the cluster-ing problem [69]. The value of K determines the numberof eigenvalues and eigenvectors considered informative.The important step is selecting the best value of K,which can be done by defining an a-priori significancethreshold (cid:15) for the decreasingly ordered eigenvalues λ i ,so that λ i > − (cid:15), ∀ ≤ i ≤ K . However, selecting thebest (cid:15) value is not always trivial, and the clustering re-sults get very sensitive to this parameter in many cases.Hence, considering the diffusion map approach of [15]and following the idea of [14], using an auxiliary dif-fusion parameter ( τ , big positive integer value) to ob-tain the exponentiated spectrum { λ τi } i =1 ,...n , the gapbetween exponentiated eigenvalues increases and sensi-tivity to the threshold value decreases very much. Usingthis new spectrum, yields to the stochastic matrix P τ ,that represents the transition matrix of a random walkin defined τ steps. The difference between thresholdingthe eigenvalues directly or the exponentiated spectrumis shown in an example in Fig. 5. As seen in this figure,selecting the best discriminative threshold value for theeigenvectors (Fig. 5c) is not easy, while with the expo-nentiated spectrum (Fig. 5d) the threshold value canbe selected in a wide range (e.g. 0 . ≤ − (cid:15) ≤ . τ need to be selected as a large positiveinteger number (e.g. 150).After selecting the value of K, the number of clustersis automatically determined using Algorithm 1. ortically Inspired Connectivity Analysis in Retinal Images 7(a) Image patch
20 40 60 8020406080 (b) Affinity matrix
10 20 3000.20.40.60.81 (c) Eigenvalues
10 20 3000.20.40.60.81 (d) Eigenvalues, τ = 150 Fig. 5: (a) A sample image patch at a crossing, (b) its affinity matrix ( A ) built upon a connectivity measure,(c) the eigenvalues ( λ i , i = 1 , ...n ) of the normalized affinity matrix ( P ) and the threshold value (1 − (cid:15) = 0 . λ τi , i = 1 , ...n ) with τ = 150 and threshold value of 0 . Algorithm 1
Spectral clustering algorithm
1: Define the affinity matrix A i,j from the connectivity ker-nel.2: Evaluate the normalized affinity matrix: P = D − A.
3: Solve the eigenvalue problem
P u i = λ i u i , where the orderof i is such that λ i is decreasing.4: Define the thresholds (cid:15) , τ and evaluate the largest integerK such that λ τi > − (cid:15) , i = 1 , . . . , K .5: Let U be the matrix containing the vectors u , . . . , u K ascolumns.6: Define the clusters k = arg max j { u j ( i ) } with j ∈{ , . . . , K } and i = 1 , . . . , n .7: Find and remove the clusters that contain less than aminimum cluster size elements. Possible neural implementations of the algorithmare discussed in [14]. Particularly, in [8, 25] an im-plementation of the spectral analysis is described as amean-field neural computation. Principal eigenvectorsemerge as symmetry breaking of the stationary solu-tions of mean field equations. In addition, in [62] itis shown that in the presence of a visual stimulus theemerging eigenvectors are linked to visual perceptualunits, obtained from a spectral clustering on excitedconnectivity kernels. In the next section the applica-tion of this algorithm in obtaining the vessel clusteringin retinal images will be presented.
In this section, the steps proposed for analyzing theconnectivities of blood vessels in retinal images and val-idating the method are described. In addition, the pa-rameter settings and the obtained results are discussedin detail. 4.1 Proposed TechniqueIn order to prove the reliability of the method in re-trieving the connectivity information in 2D retinal im-ages, several challenging and problematic image patchesaround junctions were selected. First step before de-tecting the junctions and selecting the image patchesaround them, is to apply preconditioning on the greenchannel ( I ) of a color fundus retinal image. The greenchannel provides a higher contrast between vessels andbackground and it is widely used in retinal image analy-sis. The preconditioning includes: a ) removing the non-uniform luminosity and contrast variability using themethod proposed by [29]; b ) removing the high fre-quency contents; and c ) denoising using the non-linearenhancement in SE(2) as proposed by [1] . A samplecolor image before and after preconditioning ( I enh ) areshown in Fig. 6a and 6b respectively.In next step, soft ( I soft ) and hard ( I hard ) segmen-tations are obtained using the BIMSO (biologically-inspired multi-scale and multi-orientation) method forsegmenting I enh as proposed by [1]. These images areshown in Fig. 6c and 6d respectively. The hard segmen-tation is used for detecting the junctions and selectingseveral patches with different sizes around them; whilesoft segmentation is used later in connectivity analysis.In order to find the junction locations automatically,the skeleton of I hard is produced using the morpholog-ical skeletonization technique. Then the method pro-posed by [52] is applied on this skeleton and the junc-tion locations are determined as shown in Fig. 6e. Usingthe determined locations, several image patches withsimilar sizes ( s = 10 pixels) are selected at first stage.However, as seen in Fig. 6e, some of the junctions arevery close to each other and their distances are smallerthan s/
2. For these junctions, a new patch includingboth nearby junctions (with a size equal to three times
Marta Favali et al.(a) (b) (c)(d) (e) (f)
Fig. 6: The different steps applied for selecting several image patches around junctions, (a) original RGB image,(b) enhanced image ( I enh ), (c) soft segmentation ( I soft ), (d) hard segmentation ( I hard ), (e) detected junctions andthe skeleton of the segmentation overlaid on color image, (f) selected patches overlaid on artery/vein ground truththe distance between them) is considered, and its centreis used for finding the distance of this new patch withthe other ones. These steps are repeated until no moremerging is possible or the patch size reaches the max-imum possible size (we assumed 100 as the maximumpossible value). Thus, all nearby junctions are groupedin order to decrease the number of patches that over-lap in a great extent. This results in having differentpatch sizes (0 ≤ s p i ≤ , ≤ i ≤ m ) that could in-clude more than one junction all over the image. Fig. 6fshows the junction locations and the corresponding se-lected patches overlaying on artery/vein ground truth.In order to analyze the vessel connectivities for eachimage patch ( I p i ), we need to extract the location ( x, y ),orientation ( θ ) and intensity ( f ( x, y )) of vessel pixels inthese patches. Hence for each group of junctions ( i ) withthe size s i , two patches from I enh and I soft are selected,called I enh,p i and I soft,p i respectively. Then I soft,p i isthresholded locally to obtain a new hard segmentedimage patch (called I hard,p i ). This new segmented im- age patch is different from selecting the correspond-ing patch from I hard , because I hard was obtained bythresholding the entire I soft using one global thresh-old value, but this is not appropriate at all regions. Ifthere are regions with very small vessels with a lowcontrast (often they get a very low probability of beingvessel pixels), they are normally removed in the globalthresholding approach. Accordingly, wrong threshold-ing leads into wrong tracking results e.g. C1, C2, C6in Fig. 1 are some instance patches with missing smallvessels. In this work, we selected one threshold valuefor each patch specifically using Otsu’s method [53], tokeep more information and cover a wider range of vesselpixels. Consequently, thicker vessels will be created in I hard,p i and the results will be more accurate.By knowing the vessel locations ( x, y ) other infor-mation could be extracted for these locations using I enh,p i . So f ( x, y ) equals the intensity value in I enh,p i at location ( x, y ). Moreover, by lifting I enh,p i using cakewavelets (see Eq. 1), at each location the angle corre- ortically Inspired Connectivity Analysis in Retinal Images 9 sponding to the maximum of the negative orientationresponse (real part) in the lifted domain is consideredas the dominant orientation ( θ d ) as Eq. 11. The nega-tive response is considered because the blood vessels inretinal images are darker than background. θ d = arg max θ ∈ [0 ,π ] Re ( − U f ( x, y, θ )) (11)Next step is approximating the connectivity kernels.The first kernel ( ω ), was calculated numerically. Sothe fundamental solution Γ was estimated using theMarkov Chain Monte Carlo method [59] by developingrandom paths based on the numerical solution of Eq. 6.This solution can be approximated by: x s + ∆s − x s = ∆s cos( θ ) y s + ∆s − y s = ∆s sin( θ ) θ s + ∆s − θ s = ∆sN (0 , σ ) , s ∈ , . . . , H where H is the number of steps of the random path and σ is the diffusion constant (the propagation variance inthe θ direction). This finite difference equation is solvedfor n (typically 10 ) times, so n paths are created. Thenthe estimated kernel is obtained by averaging all thesolutions [36, 62]. An overview of different possible nu-merical methods to compute the kernel is explained in[79], where comparisons are done with the exact solu-tions derived in [19, 20, 23]. From these comparisonsit follows that the stochastic Monte-Carlo implementa-tion is a fair and accurate method. The intensity-basedkernel ( ω ), the final connectivity kernel ( ω f ) and theaffinity matrix ( A ), were calculated using Eq. 7, 8 and9 respectively. Finally, by applying the proposed spec-tral clustering step in Sect. 3, the final perceptual units(individual vessels) were obtained for each patch.The above-mentioned steps for a sample crossing ina 21 ×
21 image patch are presented in Fig. 7. Afterenhancing the image (Fig. 7a), obtaining soft segmen-tation (Fig. 7b) and thresholding it locally (Fig. 7c),the vessel locations, intensity and orientation have beenextracted. As shown in Fig. 7d arteries and veins havedifferent intensities and this difference helps in discrim-inating between them. Though, orientation informa-tion is the most discriminative one. The lifted image in SE (2) using the π -periodic cake wavelets in 24 differentorientations is shown in Fig. 7h. The disentanglement oftwo crossing vessels at the junction point can be seenclearly in this figure. The dominant orientations ( θ d )for the vessel pixels are also depicted in Fig. 7e, usingline segments oriented according to the correspondingorientation at each pixel.In the next step, this contextual information (inten-sity and orientation) is used for calculating the connec-tivity kernel (Fig. 7i) and the affinity matrix (Fig. 7j) as mentioned in Sect. 2. For this numerical simulation, H , n , σ and σ have been set to 7, 100000, 0 .
05 and0 . (cid:15) and τ as0 . ×
584 ( ∼ µm/px ) and a 45 ◦ field of view. The selected patches from each image weremanually categorized into the following groups: simplecrossing (category A), simple bifurcation (category B),nearby parallel vessels with bifurcation (category C),bifurcation next to a crossing (category D) and multi-ple bifurcations (category E), and each category nar-rowed down to 20 image patches. These patches havedifferent complexities, number of junctions and sizesand they could contain broken lines, missing small ves-sels and vessels with high curvature. The parametersused in the numerical simulation of the affinity matrixand spectral clustering step (including σ , H , n , σ , σ , (cid:15) and τ ) are chosen for each patch differently, with theaim of achieving the optimal results for each case. Au-tomatic parameter selection remains a challenging taskand will be investigated in future work.Some sample figures of these cases are depicted inFig. 8. For each example, the original gray scale en-hanced image, hard segmentation (locally thresholded),orientation and intensity information, and finally theclustering result together with artery/vein labels aredepicted (Fig. 8a to 8f respectively). Although the com-plexity of these patches is quite different in all cases, thesalient groups are detected successfully. All the vesselpixels grouped as one unit have similarity in their ori-entations and intensities, and they follow the law ofgood continuation. Therefore, at each bifurcation orcrossover point, two groups have been detected.In this figure, G1 is a good example of a crossingwith a small angle. The method not only differenti-ates well between vessels crossing each other even witha small crossing angle, but it also determines the or-der of vessels, being at the bottom or passing over incrossover regions. The image patch in G2 is a good ex-ample showing the strength of the method in detecting (cid:1)(cid:1)(cid:1) y (cid:2) (cid:1) (cid:3) x (cid:4) (cid:1) (cid:3) (cid:4)(cid:1) (cid:1) (cid:4)(cid:1) (cid:2) (cid:1)(cid:1) (cid:1) (h) (cid:1)(cid:2)(cid:3)(cid:2) y (cid:4) (cid:1) (cid:5)(cid:1) (cid:1) (cid:5) (cid:3)(cid:4) (cid:1) (cid:2)(cid:2) (cid:2)(cid:5)(cid:4) x (cid:2)(cid:4) (i)
40 80 120 160 20040 80 120160200 (j) (k)
Fig. 7: A sample 21 ×
21 image patch at a crossing: (a) I enh,p i , (b) I soft,p i , (c) I hard,p i , (d) the differences inintensities are shown in color, (e) each oriented line represents the orientation at its position, (f) final perceptualunits shown in different colors (g) the ground truth artery and vein labels, (h) lifted image in SE (2), (i) connectivitykernel ( ω ), (j) affinity matrix ( A ) obtained using both orientation and intensity information, (k) thresholding theeigenvalues of the normalized affinity matrixsmall vessels. The detected small vessel in this imageis even not annotated in the artery/vein ground truth.However, this detection is highly dependent on the softsegmented image and the threshold value used for ob-taining the hard segmentation. If the small vessel is notdetected in the soft segmentation or if a hight thresh-old value is selected, then it also will not be availablein the final result. Other cases in this figure are goodrepresentations of the robustness of the method againstthe presence of a central vessel reflex (as in G3), inter-rupted lines (as in G10) or even noise (as in G9). InG9, noisy pixels are detected as individual units whichare not similar to the other groups. They can be differ-entiated from others based on their sizes. If there arevery few pixels in one group, then it can be consideredas noise and removed. There are also several cases withcomplex junctions in this figure. Presence of multiplebifurcations in one image patch, or presence of severalbifurcations close to the crossing points does not leadto wrong grouping results (as seen in G5, G6, G7, G8and G10).The parameters used during the numerical simula-tions of the image patches shown in Fig. 8 and theircorresponding sizes are presented in Table 1. For all ex-periments the values of n , (cid:15) and τ were set to 100000,0 . H , σ and σ . H and σ determine the shape Table 1: The parameters used in numerical simulationof the image patches shown in Fig. 8 and their corre-sponding sizes Name size
H σ σ G1 21 ×
21 7 0.02 0.3G2 21 ×
21 8 0.03 0.3G3 41 ×
41 10 0.03 0.1G4 39 ×
39 9 0.03 0.3G5 33 ×
33 8 0.03 0.3G6 51 ×
51 20 0.03 0.3G7 71 ×
71 17 0.07 0.3G8 73 ×
73 24 0.03 0.3G9 89 ×
89 30 0.03 0.3G10 97 ×
97 24 0.03 0.3 of the kernel. Based on the experiments, the appropri-ate value for the number of steps of the random pathgeneration is approximately 1 / σ and σ which determine the propagation variance in the θ direction and the effect of the intensity-based similar-ity term do not have a large sensitivity to variation.To quantify this, the mean and variance of these twoparameters for each of the above-mentioned categoriesare calculated and presented in Table 2. Since the se-lected patches have varying sizes and H is dependenton that, this parameter is not presented in this table.Moreover, to evaluate the performance of the method, ortically Inspired Connectivity Analysis in Retinal Images 11 G , A G , B G , C G , E G , E G , D G , D G , D G , C G , D (a) (b) (c) (d) (e) (f) Fig. 8: Sample image patches selected from the DRIVE dataset. Columns from left to right present the image patchat the green channel, segmented image, extracted orientation and intensity, clustering result and the artery/veinlabels. we introduced the correct detection rate (
CDR ) as thepercentage of correctly grouped image patches for eachcategory. These values are presented in Table 2. By con-sidering higher number of image patches per categorythe
CDR values will be more realistic.Table 2: The correct detection rate and the mean andvariance of σ and σ used in numerical simulation foreach category Category CDR% σ σ mean variance mean varianceA 85 0.032 0.0001 0.28 0.0039B 95 0.033 (cid:39) (cid:39)
0C 85 0.0269 (cid:39) (cid:39) (cid:39) If there are some high curvature vessels, then de-pending on their curvature increasing σ might help inpreserving the continuity of the vessel. As an example,G4 in Fig. 8 is relatively more curved compared to theother cases, but the clustering works perfectly in thiscase. However, for some cases it does not solve the prob-lem totally, and other kernels need to be considered forpreserving the continuity. An example 49 ×
49 imagepatch with a highly curved vessel is shown in Fig 9,where the method fails in clustering the vessels cor-rectly. The parameters used for this case are H = 16, σ = 0 .
03 and σ = 0 . σ , the distance between intensities getsa higher value and it helps in differentiating better be-tween the groups. Fig. 10 represents a sample 67 × σ from 0 . H = 24, σ = 0 .
02 and n = 100000). The other impor-tant difference between these two results is that thenoisy pixels close to the thicker vessel have been totallyremoved in the correct result. Although they seem tobe oriented with the thick vessel their intensities aretotally different. Therefore, by increasing the effect ofintensity, they are clustered as several small groups andremoved in the final step of the spectral clustering al-gorithm. In this work, we have presented a novel semi-automatictechnique inspired by the geometry of the primary vi-sual cortex to find and group different perceptual unitsin retinal images using spectral methods. Computingeigenvectors of affinity matrices, which are formed usingthe connectivity kernel, leads to the final grouping. Theconnectivity kernel represents the connectivity betweenall lifted points to the 4-dimensional feature space ofpositions, orientations and intensities, and it presents agood model for the Gestalt law of good continuation.Thus, the perceptual units in retinal images are theindividual blood vessels having low variation in theirorientations and intensities.The proposed method allows finding accurate junc-tion positions, which is the position where two groupsmeet or cross each other. The main application of theseconnectivity analyses would be in modelling the retinalvasculature as a set of tree networks. The main graphconstructed by these trees would be very informativein analyzing the topological behaviour of retinal vas-culature which is useful in diagnosis and prognosis ofseveral diseases especially in automated application inlarge scale screening programs.The detection of small vessels highly depends onthe quality of the soft segmentation, not the hard seg-mentation. These vessels could easily be differentiatedfrom noise based on the size of the group. Noisy pix-els have random orientations and intensities and theybuild smaller groups. Our method represents some lim-itations at blood vessels with high curvature. One pos-sible solution is to merge the two detected perceptualunits and form one unique unit, if there are no junctionsat these locations. The other stronger extension is touse other kernels that take into account the curvatureof structures in addition to positions and orientations.Moreover, it is also possible to enrich the affinity ma-trix with other terms e.g. the principle curvature of themultiscale Hessian (ridgeness or vesselness similarity).All these solutions will be investigated in the future.With this model we have analyzed many challeng-ing cases, such as bifurcations, crossovers, small and dis-connected vessels in retinal vessel segmentations. Thesecases not only have been reported to create tracing er-rors in the state-of-the-art techniques, but also are veryinformative for the clinical studies. Based on the resultsshown in the numerical simulations, the method is suc-cessful in detecting the salient groups in retinal images,and robust against noise, central vessel reflex, interrup-tions in vessel segments, presence of multiple junctionsin a small area and presence of nearby parallel vessels.For this reason, this can be considered as an excellent ortically Inspired Connectivity Analysis in Retinal Images 13(a) (b) (c) (d) (e) (f)
Fig. 9: Failure of clustering in presence of highly curved vessels. Columns from left to right: enhanced image; itssegmentation; orientation and intensity information; clustering result and the artery/vein labels. (a) (b) (c) (d) (e) (f) (g)
Fig. 10: The effect of including intensity term in calculating the connectivity kernel. Columns from left to right:enhanced image; its segmentation; orientation and intensity information; correct and wrong clustering results; andthe artery/vein labels. σ = 0 . Acknowledgements
This project has received funding fromthe European Union’s Seventh Framework Programme, MarieCurie Actions- Initial Training Network, under grant agree-ment n o n o . . References
1. Abbasi-Sureshjani, S., Smit-Ockeloen, I., Zhang, J.,ter Haar Romeny, B.: Biologically-inspired super-vised vasculature segmentation in SLO retinal fun-dus images. In: Image Analysis and Recognition,vol. 9164, pp. 325–334. Springer (2015)2. Al-Diri, B., Hunter, A., Steel, D.: An active contourmodel for segmenting and measuring retinal vessels.IEEE T. Med. Imaging (9), 1488–1497 (2009)3. August, J., Zucker, S.W.: The curve indicator ran-dom field: Curve organization via edge correlation.In: Perceptual organization for artificial vision sys-tems, pp. 265–288. Springer (2000)4. August, J., Zucker, S.W.: Sketches with curvature:The curve indicator random field and Markov pro- cesses. IEEE T. Pattern. Anal. (4), 387–400(2003)5. Bekkers, E., Duits, R., Berendschot, T., terHaar Romeny, B.: A multi-orientation analysis ap-proach to retinal vessel tracking. J. Math. ImagingVis. (3), 583–610 (2014)6. Boscain, U., Duplaix, J., Gauthier, J.P., Rossi, F.:Anthropomorphic image reconstruction via hypoel-liptic diffusion. SIAM J. Control Optim. (3),1309–1336 (2012)7. Bosking, W.H., Zhang, Y., Schofield, B., Fitz-patrick, D.: Orientation selectivity and the arrange-ment of horizontal connections in tree shrew striatecortex. J. Neurosci. (6), 2112–2127 (1997)8. Bressloff, P.C., Cowan, J.D., Golubitsky, M.,Thomas, P.J., Wiener, M.C.: What geometric vi-sual hallucinations tell us about the visual cortex.Neural. Comput. (3), 473–491 (2002)9. B¨uhler, K., Felkel, P., La Cruz, A.: Geo-metric methods for vessel visualization andquantification—a survey. In: Geometric Modelingfor Scientific Visualization, Math. Visual., pp. 399–419. Springer Berlin Heidelberg (2004)10. Can, A., Shen, H., Turner, J.N., Tanenbaum, H.L.,Roysam, B.: Rapid automated tracing and featureextraction from retinal fundus images using directexploratory algorithms. IEEE T. Inf. Technol. B. (2), 125–138 (1999)11. Chapman, N., Dell’Omo, G., Sartini, M., Witt, N.,Hughes, A., Thom, S., Pedrinelli, R.: Peripheralvascular disease is associated with abnormal arte- riolar diameter relationships at bifurcations in thehuman retina. Clin. Sci. (2), 111–116 (2002)12. Chutatape, O., Zheng, L., Krishnan, S.: Retinalblood vessel detection and tracking by matchedGaussian and Kalman filters. In: Engineering inMedicine and Biology Society, 1998. Proceedings ofthe 20th Annual International Conference of theIEEE, vol. 6, pp. 3144–3149. IEEE (1998)13. Citti, G., Sarti, A.: A cortical based model of per-ceptual completion in the roto-translation space. J.Math. Imaging Vis. (3), 307–326 (2006)14. Cocci, G., Barbieri, D., Citti, G., Sarti, A.: Corti-cal spatio-temporal dimensionality reduction for vi-sual grouping. Neural. Comput. (6), 1252–1293(2015)15. Coifman, R.R., Lafon, S.: Diffusion maps. Appl.Comput. Harmon. A. (1), 5–30 (2006)16. De, J., Li, H., Cheng, L.: Tracing retinal vesseltrees by transductive inference. BMC Bioinformat-ics (1), 20 (2014)17. De, J., Ma, T., Li, H., Dash, M., Li, C.: AutomatedTracing of Retinal Blood Vessels Using GraphicalModels, LNCS , vol. 7944, pp. 277–289. SpringerBerlin Heidelberg (2013)18. Delibasis, K.K., Kechriniotis, A.I., Tsonos, C., As-simakis, N.: Automatic model-based tracing algo-rithm for vessel segmentation and diameter estima-tion. Comput. Meth. Prog. Bio. (2), 108–122(2010)19. Duits, R., van Almsick, M.: The explicit solutions oflinear left-invariant second order stochastic evolu-tion equations on the 2D-Euclidean motion group.Technical Report CASA-report 43, Eindhoven Uni-versity of Technology Dep. of mathematics andcomputer science (2005)20. Duits, R., van Almsick, M.: The explicit solutions oflinear left-invariant second order stochastic evolu-tion equations on the 2D Euclidean motion group.Q. Appl. Math. (1), 27–68 (2008)21. Duits, R., Duits, M., van Almsick, M., terHaar Romeny, B.: Invertible orientation scores asan application of generalized wavelet theory. S.Mach. Perc. (1), 42–75 (2007)22. Duits, R., Felsberg, M., Granlund, G., terHaar Romeny, B.: Image analysis and reconstruc-tion using a wavelet transform constructed from areducible representation of the Euclidean motiongroup. Int. J. Comput. Vision (1), 79–102 (2007)23. Duits, R., Franken, E.: Left-invariant parabolic evo-lutions on SE(2) and contour enhancement viainvertible orientation scores—part i: Linear left-invariant diffusion equations on SE(2). Q. Appl.Math. (2), 255–292 (2010) 24. Duits, R., F¨uhr, H., Janssen, B., Bruurmijn, M.,Florack, L., van Assen, H.: Evolution equationson gabor transforms and their applications. Appl.Comput. Harmon. A. (3), 483–526 (2013)25. Faugeras, O., Veltz, R., Grimbert, F.: Persistentneural states: stationary localized activity pat-terns in nonlinear continuous n-population, q-dimensional neural networks. Neural. Comput. (1), 147–187 (2009)26. Felkel, P., Wegenkittl, R., Kanitsar, A.: Vesseltracking in peripheral CTA datasets-an overview.In: Computer Graphics, Spring Conference on,2001., pp. 232–239. IEEE (2001)27. Field, D.J., Hayes, A., Hess, R.F.: Contour integra-tion by the human visual system: Evidence for alocal “association field”. Vision Res. (2), 173–193 (1993)28. Foracchia, M., Grisan, E., Ruggeri, A.: Extractionand quantitative description of vessel features inhypertensive retinopathy fundus images. In: BookAbstracts 2nd International Workshop on Com-puter Assisted Fundus Image Analysis, p. 6 (2001)29. Foracchia, M., Grisan, E., Ruggeri, A.: Luminosityand contrast normalization in retinal images. Med.Image. Anal. (3), 179 – 190 (2005)30. Franken, E., Duits, R.: Crossing-preservingcoherence-enhancing diffusion on invertible ori-entation scores. Int. J. Comput. Vision (3),253–278 (2009)31. Fraz, M., Remagnino, P., Hoppe, A., Uyyanonvara,B., Rudnicka, A., Owen, C., Barman, S.: Blood ves-sel segmentation methodologies in retinal images –a survey. Comput. Meth. Prog. Bio. (1), 407 –433 (2012)32. Fruttiger, M.: Development of the retinal vascula-ture. Angiogenesis (2), 77–88 (2007)33. F¨uhr, H.: Abstract harmonic analysis of continuouswavelet transforms. Springer (2005)34. Gonz´alez, G., T¨uretken, E., Fleuret, F., Fua, P.:Delineating trees in noisy 2D images and 3D image-stacks. In: IEEE Conference on Computer Visionand Pattern Recognition (CVPR), pp. 2799–2806(2010)35. Habib, M.S., Al-Diri, B., Hunter, A., Steel, D.H.:The association between retinal vascular geome-try changes and diabetic retinopathy and their rolein prediction of progression-an exploratory study.BMC Ophthalmol. (1), 89 (2014)36. Higham, D.J.: An algorithmic introduction to nu-merical simulation of stochastic differential equa-tions. SIAM Rev. (3), 525–546 (2001)37. Hoffman, W.C.: The visual cortex is a contact bun-dle. Appl. Math. Comput. (2), 137–167 (1989) ortically Inspired Connectivity Analysis in Retinal Images 15
38. Hu, F.B.: Globalization of diabetes the role of diet,lifestyle, and genes. Diabetes Care (6), 1249–1257 (2011)39. Hu, Q., Abr`amoff, M.D., Garvin, M.K.: AutomatedSeparation of Binary Overlapping Trees in Low-Contrast Color Retinal Images, LNCS , vol. 8150,pp. 436–443. Springer Berlin Heidelberg (2013)40. Hubbard, L.D., Brothers, R.J., King, W.N., Clegg,L.X., Klein, R., Cooper, L.S., Sharrett, A.R.,Davis, M.D., Cai, J., in Communities Study Group,A.R., et al.: Methods for evaluation of retinal mi-crovascular abnormalities associated with hyper-tension/sclerosis in the atherosclerosis risk in com-munities study. Ophthalmology (12), 2269–2280 (1999)41. Hubel, D.H., Wiesel, T.N.: Receptive fields, binoc-ular interaction and functional architecture in thecat’s visual cortex. J. Physiol. (1), 106 (1962)42. Hubel, D.H., Wiesel, T.N.: Ferrier lecture: Func-tional architecture of macaque monkey visual cor-tex. P. Roy. Soc. Lond. B Bio. (1130), 1–59(1977)43. The International Council of Ophthalmology: ICOGuidelines for Diabetic Eye Care (2014)44. Joshi, V.S., Garvin, M.K., Reinhardt, J.M.,Abramoff, M.D.: Automated method for the iden-tification and analysis of vascular tree structuresin retinal vessel network. In: SPIE Medical Imag-ing, pp. 79,630I–79,630I. International Society forOptics and Photonics (2011)45. Koenderink, J.J., van Doorn, A.J.: Representationof local geometry in the visual system. Biol. Cy-bern. (6), 367–375 (1987)46. Lowell, J., Hunter, A., Steel, D., Basu, A., Ryder,R., Kennedy, R.L.: Measurement of retinal vesselwidths from fundus images based on 2-D modeling.IEEE T. Med. Imaging (10), 1196–1204 (2004)47. Meila, M., Shi, J.: A random walks view of spec-tral segmentation. In: Artif. Int. Stat. (AISTATS)(2001)48. Mozaffarian, D., Benjamin, E.J., Go, A.S., Arnett,D.K., Blaha, M.J., Cushman, M., de Ferranti, S.,Despres, J.P., Fullerton, H.J., Howard, V.J., et al.:Heart disease and stroke statistics-2015 update: areport from the American Heart Association. Cir-culation (4), e29 (2015)49. Mumford, D.: Elastica and computer vision, pp.491–506. Springer New York (1994)50. Narasimha-Iyer, H., Mahadevan, V., Beach, J.M.,Roysam, B.: Improved detection of the central re-flex in retinal vessels using a generalized dual-Gaussian model and robust hypothesis testing.IEEE T. Inf. Technol. B. (3), 406–410 (2008) 51. Ng, A.Y., Jordan, M.I., Weiss, Y., et al.: On spec-tral clustering: Analysis and an algorithm. Adv.Neur. In. , 849–856 (2002)52. Olsen, M.A., Hartung, D., Busch, C., Larsen, R.:Convolution approach for feature detection in topo-logical skeletons obtained from vascular patterns.In: IEEE Workshop on Computational Intelligencein Biometrics and Identity Management (CIBIM),pp. 163–167 (2011)53. Otsu, N.: A threshold selection method from gray-level histograms. IEEE T. Syst. Man. Cyb. (1),62–66 (1979)54. Parent, P., Zucker, S.W.: Trace inference, curvatureconsistency, and curve detection. IEEE T. Pattern.Anal. (8), 823–839 (1989)55. Perona, P., Freeman, W.: A factorization ap-proach to grouping, LNCS , vol. 1406, pp. 655–670.Springer Berlin Heidelberg (1998)56. Petitot, J., Tondut, Y.: Vers une neurog´eom´etrie.fibrations corticales, structures de contact et con-tours subjectifs modaux. Math. Inform. Sci. Hu-maines , 5–101 (1999)57. Poon, K., Hamarneh, G., Abugharbieh, R.: Live-vessel: Extending livewire for simultaneous extrac-tion of optimal medial and boundary paths invascular images,
LNCS , vol. 4792, pp. 444–451.Springer Berlin Heidelberg (2007)58. Quek, F.K., Kirbas, C.: Vessel extraction in medicalimages by wave-propagation and traceback. IEEET. Med. Imaging (2), 117–131 (2001)59. Robert, C., Casella, G.: Monte Carlo statisticalmethods. Springer Science & Business Media(2013)60. Sanguinetti, G., Citti, G., Sarti, A.: Image comple-tion using a diffusion driven mean curvature flowin a sub-Riemannian space. In: Proceedings of theThird Int. Conf. on Computer Vision Theory andApplications (VISIGRAPP 2008), pp. 46–53 (2008)61. Sanguinetti, G., Citti, G., Sarti, A.: A model ofnatural image edge co-occurrence in the rototrans-lation group. J. Vision (14), 37 (2010)62. Sarti, A., Citti, G.: The constitution of visual per-ceptual units in the functional architecture of V1.J. Comput. Neurosci. (2), 285–300 (2015)63. Sarti, A., Citti, G., Petitot, J.: The symplecticstructure of the primary visual cortex. Biol. Cy-bern. (1), 33–48 (2008)64. Shi, J., Malik, J.: Normalized cuts and image seg-mentation. IEEE T. Pattern. Anal. (8), 888–905(2000)65. Smith, W., Wang, J.J., Wong, T.Y., Rochtchina,E., Klein, R., Leeder, S.R., Mitchell, P.: Retinalarteriolar narrowing is associated with 5-year inci- dent severe hypertension. The Blue Mountains EyeStudy. Hypertension (4), 442–447 (2004)66. Staal, J., Abr`amoff, M.D., Niemeijer, M.,Viergever, M.A., van Ginneken, B.: Ridge-based vessel segmentation in color images of theretina. IEEE T. Med. Imaging (4), 501–509(2004)67. Tabish, S.A.: Is diabetes becoming the biggest epi-demic of the twenty-first century? Int. J. HealthSci. (2), V (2007)68. Viswanath, K., McGavin, D.M.: Diabetic retinopa-thy: clinical findings and management. CommunityEye Health (46), 21 (2003)69. Von Luxburg, U.: A tutorial on spectral clustering.Stat. Comput. (4), 395–416 (2007)70. Wagemans, J., Elder, J.H., Kubovy, M., Palmer,S.E., Peterson, M.A., Singh, M., von der Heydt,R.: A century of Gestalt psychology in visual per-ception: I. perceptual grouping and figure–groundorganization. Psychol. Bull. (6), 1172 (2012)71. Wasan, B., Cerutti, A., Ford, S., Marsh, R.: Vascu-lar network changes in the retina with age and hy-pertension. J. Hypertens (12), 1724–1728 (1995)72. Weiss, Y.: Segmentation using eigenvectors: a uni-fying view. In: The proceedings of the SeventhIEEE International Conference on Computer vi-sion, vol. 2, pp. 975–982 (1999)73. Wertheimer, M.: Laws of organization in perceptualforms. In: W. Ellis (ed.) A Source Book of GestaltPsychology, pp. 71–88. Routledge and Kegan Paul(1938)74. Williams, L.R., Jacobs, D.W.: Stochastic comple-tion fields: A neural model of illusory contour shapeand salience. Neural. Comput. (4), 837–858 (1997)75. Wong, T., Mitchell, P.: The eye in hypertension.Lancet (9559), 425–435 (2007)76. Wong, T.Y., Klein, R., Couper, D.J., Cooper, L.S.,Shahar, E., Hubbard, L.D., Wofford, M.R., Shar-rett, A.R.: Retinal microvascular abnormalities andincident stroke: the atherosclerosis risk in commu-nities study. Lancet (9288), 1134–1140 (2001)77. Xu, L., Luo, S.: A novel method for blood vesseldetection from retinal images. Biomed. Eng. Online (1), 14 (2010)78. Xu, X., Niemeijer, M., Song, Q., Sonka, M.,Garvin, M.K., Reinhardt, J.M., Abr`amoff, M.D.:Vessel boundary delineation on fundus images us-ing graph-based approach. IEEE T. Med. Imaging30