Curvature Integration in a 5D Kernel for Extracting Vessel Connections in Retinal Images
Abstract
Treelike structures such as retinal images are widely studied in computeraided diagnosis systems for largescale screening programs. Despite several segmentation and tracking methods proposed in the literature, there still exist several limitations specifically when two or more curvilinear structures cross or bifurcate, or in the presence of interrupted lines or highly curved blood vessels. In this paper, we propose a novel approach based on multiorientation scores augmented with a contextual affinity matrix, which both are inspired by the geometry of the primary visual cortex (V1) and their contextual connections. The connectivity is described with a fivedimensional kernel obtained as the fundamental solution of the FokkerPlanck equation modelling the cortical connectivity in the lifted space of positions, orientations, curvatures and intensity. It is further used in a selftuning spectral clustering step to identify the main perceptual units in the stimuli. The proposed method has been validated on several easy and challenging structures in a set of artificial images and actual retinal patches. Supported by quantitative and qualitative results, the method is capable of overcoming the limitations of current stateoftheart techniques.
I Introduction
Diabetes, as many other systematic, cardiovascular and ophthalmologic diseases, is widespread worldwide, especially in developing countries. Diabetic retinopathy is the progressive damage to the network of tiny blood vessels in the retina and it is due to the diabetes. It is the main cause of blindness and affecting the quality of life of many people in addition to raising healthcare and social costs substantially [1, 2]. Early diagnosis and treatment is essential to stop diseases progression and reduce the financial and emotional costs. Moreover, developing automated computeraided systems facilitate the diagnostic process for a larger population of people in a shorter time and at lower costs.
Based on several studies, the retinal vasculature is one of the easiesttoaccess and highly informative sources of diagnostic information not only for diabetic retinopathy, but also glaucoma, hypertensive retinopathy, and other diseases, considering them being part of the brains vasculature [3, 4]. Geometrical features such as change of vessel width, curvature, branching patterns and fractal dimension are all considered as biomarkers for clinical studies [5, 6, 7, 8]. Quantitative measurement of these features is highly dependent on correct detection and analysis of the morphologic and geometric structure of the retinal vasculature i.e. blood vessels, bifurcations and their connections.
Detection of blood vessels is often done via vessel segmentation, which is a wellstudied topic in retinal image analysis. The segmentation methods basically differentiate between blood vessels and background pixels, but they do not separate individual vessels from each other. Therefore, they have been used often as an initial step in tracking approaches (e.g. in [9, 10, 11, 12, 13, 14, 15]). Tracking the vessels makes it possible to investigate the features along individual vessels separately. The vasculature network constructed during tracking not only needs to be a correct global model, but it also should model the local connections between individual vessel segments at crossings and bifurcations correctly.
Despite all the methods proposed in literature, tracking approaches are still facing some difficulties, which are either coming from imperfect segmentation and their corresponding skeletons, or they arise at crossovers and junction points, where two vessels belonging to separate trees meet or one vessel bifurcates. Any errors created during the segmentation or skeleton extraction steps propagate to the next levels and if the method is not able to resolve these issues, it results in having wrong paths. There have been several efforts in the literature for tackling these difficulties from different points of views. In addition to great performance improvements in stateoftheart segmentation techniques (e.g. [16]), the authors in [11, 12, 13, 14] have used tubularity measurement techniques as the initial step of their method. Moreover, to solve the centreline extraction problems, which is often done via morphological thinning, several modifications and rulebased postprocessing methods have been proposed [15, 12]. On the other hand, in order to find the right connections at junction points, the geometrical configuration of the junctions has been used in local and global costs followed by optimization procedures [17, 13, 14]. In all these works, cost functions were designed in a way to avoid abrupt changes of geometrical properties (most importantly orientation) at junction points and allow a smooth transition from one vessel to the other. The method proposed by [12] is a digraphbased label propagation technique using the matrixforest theorem. The digraph weight matrix is designed in such a way that it does not allow the connected filaments (vertices) to bend too much. In the methods proposed by [11, 18] these constraints on the geometrical properties at junction points were used to find the right connections.
A different research line, aimed to extract geometrical features such as orientation or curvature from natural images, is inspired by the structure of the visual cortex. The works by [19, 20, 21] introduced the mathematical modelling of its functional architecture in terms of differential geometry. The study of [22] described the problem of context perception with psychophysical experiments, introducing the notion of association fields as the information integration along elements of images that satisfy the Gestalt law of good continuation [23, 24]. Models which take into account the structure of the cortex were proposed by [25, 26, 27], where the association fields were described with the FokkerPlanck equations. In these models, different features as orientation and curvature were considered. Then [28] proposed a model of the association fields based on the functional organization of the primary visual cortex (V1), establishing a relation between neural mechanism and image completion further developed by [29]. A similar approach was independently developed by [30, 31] who developed a Lie group theory for image analysis. This cortical method was later adapted to the problem of perceptual unit identification in [32], based on a previous study on neural synchronization [33]. In our previous work [34], we applied an instrument inspired by the geometry of the primary visual cortex (V1) and the approach of [32] to group and separate the blood vessels as individual perceptual units, even though there were informations missing due to poor segmentations. This method could find the right connections between small vessels and their parents, which are usually missing in the literature. However, the bottleneck of this method is that it is not dataadaptive and it does not follow some of the highly curved vessels.
In this work we further develop our previous model in order to solve the limitations found during the tracking of blood vessels, adding several novelties. We develop an special geometrical setting, introducing new theoretical instruments of differential geometry of vector fields. Particularly we work in a 5D space of positions, orientations, intensity and curvatures, and develop a subRiemannian structure in this setting. We introduce the feature of curvature as an important biomarker in retinal image analysis. We will see that detecting the feature of curvature in very noisy retinal images presents technical difficulties, and we will need to extend the previous highly robust singlescale detection technique developed by [6] to a multiscale approach considering that the blood vessels in retinal images have varying widths. Moreover, we apply a selftuning spectral clustering technique proposed by [35] that allows to standardize the method to automatically detect the main groups in the data without any further parameter tuning. Summarizing, the main contributions of this article are:

presenting a novel standardized and fully automatic application of the modeling of the visual cortex to the analysis of medical images;

proposing a new subRiemannian structure in 5D space of position, curvature, orientation and intensity. This allows us to detect salient perceptual units in a selftuning spectral clustering step, despite the presence of different types of variations including scale, orientation, curvature and intensity;

extending the curvature extraction technique to a multiscale approach making it suitable for analysis of multiscale blood vessels;

retrieving the challenging connections among bended, crossed, bifurcated or interrupted vessels in noisy, low quality or diseased retinal images;

demonstrating potential extensions of the method by applying it successfully on several curvilinear synthetic images with various complexities.
The structure of the paper is the following. We will start in section II with introducing the fivedimensional feature space, the new model of the cortical connectivity in this space, and the steps for lifting 2D retinal images. The qualitative and quantitative results of proposed experiments on two different sets are discussed in section III. Finally, the article is concluded in section IV.
Ii Theory and Methodology
In this section, we present: a new algorithm based on a higher dimensional feature space, the model of cortical connectivity in this integrated 5D space and its application in the definition of an affinity matrix. Finally an automatic clustering algorithm, which takes as input our affinity matrix, and a new multiscale curvature extraction technique for retinal images are presented. Additionally, we explain the required preprocessing steps for preparing and lifting the retinal images in 5D space.
Iia Lifting the Image to the Orientation Cortical Space
A regular curve in the twodimensional plane can be represented by . The tangent vector to the curve may be indicated as:
(1) 
Based on previous studies [36, 31], the twodimensional curves can be lifted to the () space of positions and orientations () using the direction of the tangent vector as:
(2) 
where . Therefore, we can write the following:
(3) 
where , are leftinvariant vector fields with respect to the group law of . The neural interpretation of this is that the lifted curves in the cortical space are the integral curves of the two vector fields . The curvature term indicates the rate of change of the orientation and determines the shape of the curve. The Lie algebra is generated by and coincides with , where .
IiB Integrated Model of Position, Orientation, Curvature and Intensity
In this section, we describe a new model of the cortical connectivity in a 5D feature space of position, orientation, intensity and curvature. In our previous model [34] curvature was not present. We will see that this new feature will significantly change the structure of the space, since there is a strong interplay between orientation and curvature. In this way the model will allow strongly curved connectivity patterns.
We assume that the twodimensional curve in cortical plane () is lifted to a 5D space of positions, orientations, intensity and curvature (). Thus, the lifted curve may be written as:
(4) 
Similar to (3), we will have:
(5) 
where and . By defining new vectors in the 5D space as:
(6) 
we are able to write in terms of these vectors:
(7) 
In general, the solution of the following differential equation represents the curves in this lifted domain:
(8) 
The horizontal distribution of planes is . Denoting the directional derivatives in the direction of the vectors , the commutators of these vectors are as follows:
(9) 
Then the Lie algebra generated by is the whole space at every point.
The cortical connectivity can be modeled by a stochastic counterpart of (7). The Markov process that results from the Brownian motion with randomly curved paths has been introduced by [37]. The process is represented by the following differential equations:
(10) 
where is a normally distributed variable with zero mean and variance equal to . If denotes the probability density to find a particle at a point with a certain direction , intensity and curvature , at a specific time , then the FokkerPlanck equation describing the diffusion of the particle density will be:
(11) 
so that and . This partial differential equation means that a particle at a point transports in the direction of in the 5D space. There is no transport in the or direction, but the diffusion in the direction indicates the rate of transport in the direction. The diffusion in direction is independent from the other variables. The nonnegative fundamental solution of (11) satisfies the following equation:
(12) 
A good estimate of this solution is a section of fundamental solutions with fixed ( as a 3D kernel) symmetrized and multiplied to an exponential term which considers the closeness between two points located in the different intensity planes, times an exponential term in the curvature planes. Hence, this new connectivity kernel is presented as:
(13) 
IiC Numerical Approximation of the 5D Kernel
The kernel is numerically estimated using the general Markov Chain Monte Carlo method technique. We say a few words on how to adapt it to the 5D case for reader convenience. The system in equation (10) can be approximated by:
(15) 
where is the number of steps of the random path and is a generator of numbers taken from a normal distribution with mean 0 and standard deviation of . The stochastic path is obtained from the estimate of the kernel as the average of their passages over discrete volume elements, solving this finite difference equation times [32]. The affinity matrix described in (14) is evaluated from this kernel.
Fig. 1 represents two level sets of the 5D kernel by fixing (Fig. 0(a)) and (Fig. 0(b)) dimensions. The 2D projection on by summation over all orientations of five different 4D stochastic kernels having different curvature () values is also presented by Fig. 0(c). The intensity term is kept constant for all figures. As seen in these figures, by increasing the absolute value of curvature, the shape of the kernel also changes and it deviates from the elongated shape.
IiD Lifting Retinal Images in the 5D Space
In this section, we explain how we lift the 2D retinal images to our proposed feature space of positions, orientations, curvatures and intensity in practice.
IiD1 Orientation Score Transform
Being inspired and supported by several electrophysiological experiments, the receptive profiles of simple cells in the primary visual cortex are often interpreted as oriented Gabor filters or directional derivatives of the Gaussian filters [39, 40, 41, 36]. Moreover, [31, 30] introduced the invertible orientation score (OS) transformation for lifting the image to the space. The invertibility property prevents information loss during the transformation, which is guaranteed by certain requirements of the kernel used for the transformation. The cake wavelets introduced by [42, 43] are proper wavelets, which satisfy these criteria. Similar to the Gabor wavelets, they are quadratic anisotropic filters, but unlike these, their summed Fourier transformations cover the entire frequency domain, making them spatially scaleindependent. Therefore, it is possible to disentangle the crossing elongated structures in a 2D image from each other regardless of their scale and without loss of image evidence. In order to construct the orientation score , the image () is convolved with rotated and translated versions of a mother wavelet () as
(16) 
where is the 2D counterclockwise rotation matrix, the overline denotes the complex conjugate and denotes the convolution [31, 30].
Applying this tool for lifting the retinal images, considering the Gaussian profile of the blood vessels and the fact that the blood vessels have a darker intensity compared to the background, the real part of the negative OS response () is used for obtaining the orientation in these images. Because of the quadratic property of the cake wavelets, the real symmetric part of the transformation performs as a vessel ridge detector. Moreover, for each location the response may be different from zero in more than one orientation; therefore, the orientation at which the maximum response is obtained at location is the dominant orientation at that location:
(17) 
so that only takes discrete values.
IiD2 Curvature Extraction Technique
As mentioned before, the curvature determines the rate of orientation change (). If the parametric representation of the curves in the image is available, then by differentiating and once more in (1), we will have:
(18) 
where is the arclength of the curve. So the curvature can be computed as:
(19) 
In computer vision and more specifically in retinal image analysis, several methods have been introduced for measuring the curvatures of curvilinear structures (i.e. blood vessels) [44]. The classical methods that measure the curvatures locally need an initial segmentation, centreline extraction, and separation of segments located between junctions. It is then followed by fitting curves to the segments and by curvature measurement using (19) (e.g. [45]). The drawback of all these methods is their dependency to initial preprocessing and segmentation steps, which may contain errors and missing information. More importantly, the curvature information is not available for junction points because it is not possible to fit a curve to these points where more than one elongated structure meet.
To solve these problems, [6] proposed a local curvature measurement technique by locally fitting exponential curves [46] to the lifted image in . The exponential curves in are interpreted as straight lines considering the curved geometry of and they have constant tangent vectors relative to the rotating frame . The tangent vectors of the exponential curve that best fit the data in the lifted image are obtained by eigensystem analysis of the Gaussian Hessian (expressed in the rotating frame). Then they directly define the curvature value of their spatial projections [47]. This approach makes it possible to assign to each location and orientation in the lifted image a curvature value, without needing explicit curve parameterizations. Such curvature maps (on SE(2)) can be projected on the plane whereby only one value of curvature value is assigned to each spatial location in the image. Finally, these 2D curvature maps can be filtered in a later stage by a vessel confidence map such as , as a Laplacian ridge detector at a single scale () (See ch. 5.3, [47] for detailed explanations).
In order to obtain the curvature value for each location () in our 5D lifted retinal images an altered version of the abovedescribed method is used. The method performs the best when the spatial scale () of the Gaussian derivatives used to calculate the Gaussian Hessian matrix matches the vessel width. By using a single scale, the curvature values are accurate only for the vessels which their width match the used scale. Since the blood vessels in retinal images have different widths, using a multiscale approach helps in covering various vessel widths available in the image and getting more accurate results. Therefore, we altered the method to a multiscale approach. To do so, confidence and curvature maps have been obtained for several scales ( and , ). Then the final curvature map is constructed by assigning to each pixel the curvature value that corresponds to the scale at which the largest confidence value has been obtained, i.e. :
(20) 
IiD3 The Intensity Information
For a retinal image , an appropriate preprocessing step is needed to normalize the luminosity and contrast all over the image, reduce noise and to enhance the blood vessels. We use the same preprocessing technique proposed by [48] on the green and red channels of the original image separately and combine them as , where and are preprocessed red and green channels. This step is essential as it helps in increasing the difference between the intensity of arteries and veins and helps in better discrimination of blood vessels crossing each other. The final value of intensity at position of image is considered as .
IiD4 Set of Interest Positions
As mentioned in section I, vessel segmented images may contain vessel disconnections or wrongly detected vessel pixels. However, we use the segments to get an initial estimation of the location of the blood vessels and only focus on those locations. This helps in reducing the computation complexity as the size of affinity matrix (14) reduces tremendously. The outcome of segmentation can be either a probability map showing the probability for each pixel to be part of a vessel or it can be a deterministic binary map, in which each pixel is labeled as a vessel (label ) or background (label ). In this work, we use the deterministic binary segmentations obtained by the method proposed by [48] for both RGB and SLO retinal images. If the vessel segmented image is called as then the vessel locations are found as:
(21) 
This set provides a set of interest points as we only consider the information at these locations for our connectivity analysis.
To summarize this section, for an image , after preprocessing the image and obtaining , we use the orientation score transform and the segmentation technique on to get and respectively. Then by finding the set of interest points using (21), we create a 5D feature map called as:
(22) 
IiE An Efficient Cortically Inspired Spectral Clustering
In modern data analysis clustering is an important research topic and many algorithms have been proposed in this area (e.g. [49, 50, 51]) where the eigenvectors of the graph Laplacian are used to analyse and reveal the cluster structure of the data. In our previous work [34], the problem of grouping the blood vessels in retinal images (interpreted as perceptual units) was solved by proposing a semiautomatic approach based on spectral analysis of the affinity matrix. So the salient groups were obtained as the first eigenvectors corresponding to the largest eigenvalues. However a fixed threshold was defined for selecting the largest eigenvalues. In order to avoid the manual intervention, we follow the method proposed by [35] in which this problem is solved.
Unlike more basic clustering methods, as kmeans, where the number of clusters and interaction distances of the nodes in the dataset must be assumed a priori, this algorithm automatically and optimally tunes these parameters. In this algorithm, the structure of the eigenvectors is used to determine the number of groups. A cost function is defined and evaluated from the alignment of the eigenvectors. Consecutively, the best number of clusters is considered as the one which minimizes the cost function. Correspondingly, the best clustering quality that has a reverse relation to the alignment cost is obtained (see [35] for detailed explanations). In the final step, the noisy elements that construct small sized groups are removed. From a neural approach, the fact that in the presence of a visual stimulus the perceptual units are defined by the emerging eigenvectors is described by [32].
The full pipeline is presented in Algorithm 1. We assume an image has been given as input. The first step is to lift the image to the 5D feature space (as explained in section IID). Then we need an affinity matrix which is created after modelling the cortical connectivity in the 5D space as we proposed. Finally, the selftuning spectral clustering algorithm returns the final perceptual units in that image , where is the final number of clusters and includes set of points in image which belong to cluster .
Fig. 2 depicts a sample application of the proposed method for clustering the perceptual units in both an artificial image (the first row) and a small patch of retinal images (the second row). The synthetic image includes three crossing circles with different radii and corresponding curvatures, and the retinal patch includes two crossing vessels. For each case, the orientation and curvature of the lifted images are shown in the column (\subreffig:osmap3D) and (\subreffig:CurvMap3D) respectively. The intensity value for the synthetic image is constant all over the circles and for the retinal patch, it is colorcoded in (\subreffig:Orig). A sample level set of corresponding 5D kernels (while keeping two dimensions ( and ) fixed) has been shown in each row (\subreffig:Isosurf). Finally, (\subreffig:Results) shows the three detected groups in the synthetic image and three vessels in the retinal image in different colors.





Iii Experiments
In this section, we present a potential application of the proposed connectivity analysis for solving the aforementioned problems in curvilinear structure tracking methods for retinal vasculature analysis (see section I). After explaining the material used for validating the method, the details of the numerical simulation are described. Then the quantitative and qualitative results of the proposed technique are presented and discussed in detail.
Iiia Material
Two datasets have been used for validating the method. The specifications and the preparation steps of each dataset are explained in detail as follows.
IiiA1 Phantom images
A set of synthetic images () with known orientations and curvature values and constant intensities has been generated to include various rotated, curved and interrupted vessellike structures. Five different groups are created to mimic possible structures that could be present in retinal images, similar to the categories proposed in our previous work [34]. These categories are {enumerate*}[label=()]
crossings;
bifurcations;
close parallel vessels;
bifurcations and crossings; and
vessels with multiple nearby bifurcations . Each of these categories may also include challenging structures. For instance, they may be interrupted or highly curved or include small junction (crossing/ bifurcation) angles. In order to differentiate between the simple and the challenging cases for each category, we name group as if it is challenging. Fig. 4 depicts ten different phantom images in each row of the first column, two per category, together with their colorcoded orientation and curvature maps (second and third columns).
The basic element used for creating these phantom images is a sine wavelike structure, which is generated with several frequencies and amplitudes and it is rotated and located at different positions depending on the target shape. By adjusting the frequency and amplitude of the waves, different curvature values can be created. In addition to the vessellike structures, other challenging structures such as dashing and the Euler spiral have been also used to examine the strength of the method in grouping these curved structures (e.g. Fig. 4, A1).
IiiA2 Retinal images
The public IOSTAR^{1}^{1}1Available at: http://www.retinacheck.org/datasets dataset contains images captured using scanning laser ophthalmoscope (SLO) technology^{2}^{2}2iOptics BV, the Netherlands. These high contrast images have a resolution of with 45 FOV. The blood vessels, junctions and artery/vein labels have been annotated for 24 images and corrected by two different experts in order to decrease the interuser variability. To reduce the computation complexity these images and their annotations are downsampled to half size . Our main purpose is to emphasize the strength of the method in individuation of vessels, specifically at junction points or challenging cases, where most of the proposed methods in the literature face difficulties. Since the images in this set are very noisy and challenging for segmentations, they are good cases for examining our proposed method in detecting missing connections.
Following the procedure described in section IID, we lift each image to the 5D space after preprocessing and segmenting the images. The scales used for obtaining the curvature map for these SLO images are in pixels and . Even though we lift full images, we only consider the information around junction points. In order to crop small patches, we detect the junction points in each image using the hybrid method proposed in [52], and call the set of junction locations as , where is the number of detected junctions in image . If we assume that image is lifted to , then the cropped patch at location () is defined as:
(23) 
where determines the size of the neighborhood considered around each junction. We fix pixels for all patches. At the end, the size of each 5D patch is . In total 272 random patches around junctions points, with resolution of have been selected from this set. Similar to the phantom images, these patches are also categorized in five different groups depending on their structure and complexity. Fig. 3 depicts a sample SLO image of the IOSTAR dataset (\subreffig:colored). It also represents the colorcoded orientation (\subreffig:sloOSMap), confidence (\subreffig:sloConfMap), curvature (\subreffig:sloCurvMap) maps, detected junction points overlaid on the preprocessed image (\subreffig:Junctions), the corresponding neighborhood around each junction location (\subreffig:Patches) and the arteryvein labels of the vessels (\subreffig:vesselGT), which are later used for validation. Note that the depicted confidence and curvature maps are related to one single scale , and the absolute curvature value is shown.
IiiB Validation
IiiB1 Phantom images
To validate the method using the phantom images, the method presented in Algorithm 1 is directly applied on the lifted images, starting from Step 2. The orientation and curvature values of these images are available from the beginning. The intensity is constant for these images, so the intensity term in (13) is equal to 1. The last two columns of Fig. 4 represent two different clustering results per image. For each image, two kernels have been used for obtaining the affinity matrix: the new kernel (adaptive, based on the curvature at each point); and the kernel used in our previous work [34]. This helps in highlighting the importance of including additional contextual information for connectivity analysis. Each color in the final results represents one detected unit. The parameter used in these simulations is 0.01 for the 4D kernel [34] and 0.001 for the 5D one. As seen in this figure, the new method is capable of grouping the elongated, rotated and curved structures despite disconnections, high curvature points or small crossing angles. It not only differentiates well between curved crossing structures, but also groups the bifurcations within the main parent structure so that they construct one unique unit. For other examples we refer to Fig. 1 in the supplementary material.

IiiB2 Retinal patches
To validate the method on retinal image patches, after creating the patches in the lifted domain (23), we perform the analysis starting from Step 2 in Algorithm 1. All 272 patches have been examined semiautomatically. The number of patches processed per group is presented in the second column of Table I. In order to see if each detected unit represent one individual vessel, we create the corresponding 2D patch for each image from the arteryvein labeled images, which are available for this dataset. Existing vessels in each arteryvein labeled patch are compared with detected clusters. In case two vessels with similar labels in the ground truth image (artery or vein) belong to separate vessel trees (parents), they get a different label. In addition, the operator performs a final check to control the final results. Two criteria are used to evaluate the performance. One is the Correct Detection Rate (CDR%) defined as the percentages of correctly grouped patches among all examined cases, similar to [34]. The second criterion is obtained in the spectral clustering step (explained in IIE). In this method, the best number of clusters is the one which minimizes the defined cost function or maximizes the quality (). This criterion represents how well aligned the elements of each group are. These two performance values have been measured for all the patches and presented in the last two columns of Table I for each category separately.
During the experiments the parameters used in numerical simulations and the calculation time of different parts of the experiments have been recorded. Several parameters are involved in creating the kernel. Some of them are determined automatically based on the available information in the data and others need to be set manually. The first set includes the size of the kernel in and dimension ( and respectively). The other one is the number of discrete curvature values (), which for each a 3D section of fundamental solution with fixed is created. Considering the step size as 1 pixel, and are determined by the difference between maximum and minimum coordinates of the vessel locations in and directions. Similarly, considering the step size of for discrete curvature values, is obtained by division of the difference between maximum and minimum of the available curvature values in the patch over the step size. Moreover, the number of steps used in generating the random paths is set automatically as one third of the patch size. The second set of parameters is set manually. The number of discrete orientations (), number of iterations in the MonteCarlo simulation () and (used in (13)) were set to 18, 100000 and 1 respectively and kept constant for all the cases. The other parameters are presented in Table I for each category separately and for all the cases (in mean standard deviation format).
Group  Size  Parameters  Performance  

A  16  0.8125  0.9866  
B  42  0.7143  0.9979  
C  48  0.8542  0.9978  
D  31  0.7742  0.9974  
E  19  0.7368  0.9967  
A1  9  0.6667  0.9756  
B1  23  0.6522  0.9987  
C1  41  0.8537  0.9975  
D1  27  0.8148  0.9974  
E1  18  0.7222  0.9967  
All  274  0.7602  0.9942 
For recording the calculation times the whole process has been divided into four steps: the discretization step before creating the kernel; creating the kernel for several curvature values; creating the affinity matrix, and the spectral clustering step. The times are called , , and , respectively, and they are affected by several parameters including the number of vessel pixels in each patch (), the number of discrete orientations (), curvatures () and the dimension of the kernel in and dimensions ( and ). To consider these effects, the weighted average of calculation times for each image patch is obtained, so that the weight for each timing is defined as the product of the affecting parameters as. Thus we define the weights as , , , where indicates the patch number among patches. It is worth mentioning that since the final number of clusters is determined by comparison among the clustering costs of several cluster sizes (), is additionally affected by the number of examined cluster sizes ( for all the cases), so the final clustering time is divided by . The size of the affinity matrix for patch is . These weighted times and their normal average over all patches per step are presented in Table II.
mean(s)  0.06  43.26  17.86  0.86 
weighted mean (s)  0.06  60.64  17.73  1.06 
A set of sample results for various kinds of patches in challenging groups is depicted in Fig. 5. In this figure, the first column shows the cropped patch from the artery/vein ground truth image. The second column depicts the colorcoded normalized intensity values taken from the preprocessed image. As seen in these figures, the variation of intensity is high even for small children vessels belonging to one parent vessel. The third and fourth columns represent the colorcoded orientations and curvature values. Finally, the last column represents the clusters found in each patch, each shown in an individual color. Other examples are presented in Fig. 2 in the supplementary material.

The presented results, parameters and timings are discussed in the next section.
IiiC Discussion
Based on the results presented in the previous section, the main advantage of the new method is that by including the curvature information as an additional contextual information, the kernel adapts itself naturally according to the available data. If the curvature is high, the kernel rotates as well, otherwise it finds a closer path to the points which are collinear with respect to the reference point. In both datasets, the bifurcations are grouped with the parent vessel, but at crossovers with small crossing angles, despite their similar appearance to junctions, the vessels are totally separated. The main reason is that the curvature at junction points is high (because of sudden change of orientation), while for crossings the orientation for individual vessels changes only slightly (in most of the cases). This is advantageous not only in differentiating between junctions and crossings, but also in separation of arteries from veins or crossing tree structures from each other. Fig. 6 shows the clustering results for two retinal patches obtained using the new 5D kernel and the previously introduced 4D kernel by [34]. This helps in depicting the differences between the two methods visually.

As presented in Fig. 5, the intensity is a less informative feature compared to the geometrical features because of its large variation within a small neighborhood; however, in some cases, it is a good local criterion differentiating arteries from veins. Thus it is also included in the final affinity matrix with a smaller effect (using a relatively large ).
Some examples of the limitations of the method in clustering the phantom and retinal patches are represented in Fig. 7 and 8 respectively. For phantom cases, the presence of very high curvature combined with the cocircularity and colinearity of vessels does not allow to obtain a good clustering result. Considering other features, as the intensity, could be helpful in solving this problem. However, as shown in the top row of Fig. 8, the feature of intensity is not useful for correct clustering. In this image, one of the bifurcations has been assigned as a vessel crossing the other one because it is almost orthogonal to its parent vessel; while the other crossing vessel has been wrongly clustered as a bifurcation. In the bottom row, one of the small bifurcations is totally missing in the segmentation and the other small one is not clustered with its parent vessel because of lack of information.








The statistical analysis on the parameters used during numerical simulations is presented in Table I. Based on these results the curvature diffusion constant parameter () changes in a small range for simple and challenging cases per group, except for category A and A1, considering the fact that the number of available challenging patches in category A is small compared to the others. The parameter has a small variation as well. Therefore, it is reasonable to use the mean value in general for examining new patches in each group. It is worth mentioning that all the patches examined in this work have been selected from retinal images with the same resolution and respective pixel size. If the pixel size increases or decreases, the controlling the size of the 3D kernel needs to be increased or decreased accordingly. It is also important to mention that compared to previous work, the eigensystem analysis is fully automatic in this work, due to the use of selftuning spectral clustering. Therefore, no additional parameters need to be tuned for this step.
Qualitative and quantitative results indicate the better performance of the method on all kinds of retinal patches and challenging structures. Compared to the previous work [34], the performance values have changed for some groups. There are two main reasons. On one hand, in the previous work, the was calculated for 20 patches per group; while in this work each group is categorized into two groups depending on the available structure and also the number of patches per group is different. On the other hand, in this work we evaluated the performance with the assumption that the bifurcations need to be grouped with the main parent vessel, while in previous work, due to use of one elongated kernel, the assumption was to have at least two separate units depending on the bifurcation angle. Another minor difference is that the dataset has changed and the patches have been selected from a different set of retinal images. Therefore, it is not fair to make a onebyone comparison to the previous results.
Last but not least point is about the computation times. The codes are implemented in Matlab and the times are measured on an Apple Macbook Air, Intel Core i7, 1.7 GHz processor and 8GB of memory. The most time consuming step as presented in Table II is the calculation of several 3D kernels. The number of kernels depends on the number of discrete curvature values () existing in the data, and the computation time as mentioned before depend on its size. However, the kernels can be calculated only once for all possible curvature values and then used for analyzing all other patches. The next most time consuming part is the calculation of the affinity matrix which is performed per pair of points. The weighted average times are good indicators of the complexity of each step. Although they are already still relatively small, they can be improved both from hardware and implementation points of view.
Iv Conclusion
In this work a new method for analyzing the connectivities in images containing elongated, rotated and curved structures is proposed. The new connectivity model is inspired by the geometry of the primary visual cortex, where the connectivity between all the lifted points into the fivedimensional space of positions, orientations, intensity and curvatures is represented by a 5D kernel. We use two sets of patches: a retinal and an artificial one, to identify perceptual units in these figures, showing that this can be considered as a good quantitative model for their constitution and in general for the law of good continuation. Including the geometrical contextual properties in addition to local features, such as intensity, makes the method very robust against different kinds of variations and missing information which could exist in clinical images. We analyze different challenging cases that are very informative for clinical studies and we show how the proposed method is successful in solving the problem of grouping, also for blood vessels with high curvature. That was a limitation in the previously proposed method [34] and many of the stateoftheart techniques.
In this work, we only analyzed small patches. However, by replacing the selftuning spectral clustering step with label propagation methods or using the affinity matrix in an optimization framework (e.g. integer programming), the constitution of the entire vasculature network in full retinal images becomes possible. When errors appear in the lifting step due to noise or low contrast regions, and if the measured curvature or orientations are not accurate enough, then the method fails. Therefore, it is essential to validate the curvature and orientation measurement methods in advance.
The proposed method has a great potential in discrimination and separation of arteries from veins in retinal images and, in a general view, separation of all the tree structures crossing each other in the vasculature network. Most of the segmentation or artery/vein separation methods are local, pixel based techniques which do not take into account the global connectivity of the blood vessels in the network. By including this global connectivity criterion, most of the errors and wrong detections will be eliminated and the problem of missing information will be handled appropriately. This technique is also highly suitable for other application areas, such as finding roads in aerial photographs or detecting cracks and faults in materials.
Acknowledgment
We thank Dr. Erik Bekkers who provided us some parts of his codes for measuring the curvatures in retinal images and creating the phantom images.
References
 [1] K. Viswanath and D. M. McGavin, “Diabetic retinopathy: clinical findings and management,” Community Eye Health, vol. 16, no. 46, p. 21, 2003.
 [2] F. B. Hu, “Globalization of diabetes: the role of diet, lifestyle, and genes,” Diabetes Care, vol. 34, no. 6, pp. 1249–1257, May 2011.
 [3] L. D. Hubbard, R. J. Brothers, W. N. King, L. X. Clegg, R. Klein, L. S. Cooper, A. R. Sharrett, M. D. Davis, J. Cai, A. R. in Communities Study Group et al., “Methods for evaluation of retinal microvascular abnormalities associated with hypertension/sclerosis in the atherosclerosis risk in communities study,” Ophthalmology, vol. 106, no. 12, pp. 2269–2280, 1999.
 [4] M. Foracchia, E. Grisan, and A. Ruggeri, “Extraction and quantitative description of vessel features in hypertensive retinopathy fundus images,” in Book Abstracts 2nd International Workshop on Computer Assisted Fundus Image Analysis, E. Grisan, Ed., 2001, p. 6.
 [5] F. Huang, J. Zhang, E. J. Bekkers, B. Dashtbozorg, and B. M. ter Haar Romeny, “Stability analysis of fractal dimension in retinal vasculature,” in Proceedings of the Ophthalmic Medical Image Analysis Second International Workshop, OMIA, Held in Conjunction with MICCAI, ser. Lecture Notes in Computer Science, Munich, Germany, October 9 2015, pp. 1–8.
 [6] E. J. Bekkers, J. Zhang, R. Duits, and B. M. ter Haar Romeny, “Curvature based biomarkers for diabetic retinopathy via exponential curve fits in SE(2),” in Proceedings of the Ophthalmic Medical Image Analysis Second International Workshop, OMIA, Held in Conjunction with MICCAI, ser. Lecture Notes in Computer Science, Munich, Germany, October 9 2015, pp. 113–120.
 [7] N. Chapman, G. Dell’Omo, M. Sartini, N. Witt, A. Hughes, S. Thom, and R. Pedrinelli, “Peripheral vascular disease is associated with abnormal arteriolar diameter relationships at bifurcations in the human retina,” Clinical Science, vol. 103, no. 2, pp. 111–116, 2002.
 [8] M. S. Habib, B. AlDiri, A. Hunter, and D. H. Steel, “The association between retinal vascular geometry changes and diabetic retinopathy and their role in prediction of progressionan exploratory study,” BMC Ophthalmology, vol. 14, no. 1, p. 89, 2014.
 [9] X. Xu, M. Niemeijer, Q. Song, M. Sonka, M. K. Garvin, J. M. Reinhardt, and M. D. Abràmoff, “Vessel boundary delineation on fundus images using graphbased approach,” IEEE Transactions on Medical Imaging, vol. 30, no. 6, pp. 1184–1191, 2011.
 [10] J. De, H. Li, and L. Cheng, “Tracing retinal vessel trees by transductive inference,” BMC Bioinformatics, vol. 15, no. 1, p. 20, 2014.
 [11] J. De, T. Ma, H. Li, M. Dash, and C. Li, Automated tracing of retinal blood vessels using graphical models, ser. Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 2013, vol. 7944, pp. 277–289.
 [12] J. De, L. Cheng, X. Zhang, F. Lin, H. Li, K. Ong, W. Yu, Y. Yu, and S. Ahmed, “A graphtheoretical approach for tracing filamentary structures in neuronal and retinal images.” IEEE Transactions on Medical Imaging, vol. 35, no. 1, p. 257, 2016.
 [13] E. Türetken, G. González, C. Blum, and P. Fua, “Automated reconstruction of dendritic and axonal trees by global optimization with geometric priors,” Neuroinformatics, vol. 9, no. 23, pp. 279–302, 2011.
 [14] E. Türetken, F. Benmansour, and P. Fua, “Automated reconstruction of tree structures using path classifiers and mixed integer programming,” in IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2012, pp. 566–573.
 [15] B. Dashtbozorg, A. M. Mendonça, and A. Campilho, “An automatic graphbased approach for artery/vein classification in retinal images,” IEEE Transactions on Image Processing, vol. 23, no. 3, pp. 1073–1083, 2014.
 [16] G. Azzopardi, N. Strisciuglio, M. Vento, and N. Petkov, “Trainable COSFIRE filters for vessel delineation with application to retinal images,” Medical Image Analysis, vol. 19, no. 1, pp. 46–57, 2015.
 [17] Q. Hu, M. D. Abràmoff, and M. K. Garvin, “Automated construction of arterial and venous trees in retinal images,” Journal of Medical Imaging, vol. 2, no. 4, pp. 044 001–044 001, 2015.
 [18] L. Cheng, J. De, X. Zhang, F. Lin, and H. Li, “Tracing retinal blood vessels by matrixforest theorem of directed graphs,” in Medical Image Computing and ComputerAssisted Intervention, ser. Lecture Notes in Computer Science. Springer, 2014, no. 8673, pp. 626–633.
 [19] J. J. Koenderink and A. J. van Doorn, “Representation of local geometry in the visual system,” Biological Cybernetics, vol. 55, no. 6, pp. 367–375, 1987.
 [20] W. C. Hoffman, “The visual cortex is a contact bundle,” Applied Mathematics and Computation, vol. 32, no. 2, pp. 137–167, 1989.
 [21] ——, “The Lie algebra of visual perception,” Journal of Mathematical Psychology, vol. 3, no. 1, pp. 65–98, 1966.
 [22] D. J. Field, A. Hayes, and R. F. Hess, “Contour integration by the human visual system: Evidence for a local “association field”,” Vision research, vol. 33, no. 2, pp. 173–193, 1993.
 [23] J. Wagemans, J. H. Elder, M. Kubovy, S. E. Palmer, M. A. Peterson, M. Singh, and R. von der Heydt, “A century of Gestalt psychology in visual perception: I. perceptual grouping and figure–ground organization.” Psychological Bulletin, vol. 138, no. 6, p. 1172, 2012.
 [24] M. Wertheimer, “Laws of organization in perceptual forms.” in A Source Book of Gestalt Psychology, W. Ellis, Ed. Routledge and Kegan Paul, 1938, pp. 71–88.
 [25] D. Mumford, Elastica and computer vision. New York, NY: Springer, 1994.
 [26] L. R. Williams and D. W. Jacobs, “Stochastic completion fields: A neural model of illusory contour shape and salience,” Neural Computation, vol. 9, no. 4, pp. 837–858, 1997.
 [27] J. August and S. W. Zucker, “The curve indicator random field: Curve organization via edge correlation,” in Perceptual Organization for Artificial Vision Systems, ser. The Kluwer International Series in Engineering and Computer Science, K. L. Boyer and S. Sarkar, Eds. Boston, MA: Springer, 2000, vol. 546, pp. 265–288.
 [28] G. Citti and A. Sarti, “A cortical based model of perceptual completion in the rototranslation space,” Journal of Mathematical Imaging and Vision, vol. 24, no. 3, pp. 307–326, 2006.
 [29] G. Sanguinetti, G. Citti, and A. Sarti, “Image completion using a diffusion driven mean curvature flow in a subRiemannian space,” in Int. Conf. on Computer Vision Theory and Applications (VISAPP 2008), Funchal, 2008.
 [30] R. Duits, M. Felsberg, G. Granlund, and B. ter Haar Romeny, “Image analysis and reconstruction using a wavelet transform constructed from a reducible representation of the Euclidean motion group,” International Journal of Computer Vision, vol. 72, no. 1, pp. 79–102, 2007.
 [31] R. Duits, M. Duits, M. van Almsick, and B. ter Haar Romeny, “Invertible orientation scores as an application of generalized wavelet theory,” Pattern Recognition and Image Analysis, vol. 17, no. 1, pp. 42–75, 2007.
 [32] A. Sarti and G. Citti, “The constitution of visual perceptual units in the functional architecture of V1,” Journal of Computational Neuroscience, vol. 38, no. 2, pp. 285–300, 2015.
 [33] A. Sarti, G. Citti, and M. Manfredini, “From neural oscillations to variational problems in the visual cortex,” Journal of PhysiologyParis, vol. 97, no. 2, pp. 379–385, 2003.
 [34] M. Favali, S. AbbasiSureshjani, B. ter Haar Romeny, and A. Sarti, “Analysis of vessel connectivities in retinal images by cortically inspired spectral clustering,” Journal of Mathematical Imaging and Vision, vol. 56, no. 1, pp. 158–172, 2016.
 [35] L. ZelnikManor and P. Perona, “Selftuning spectral clustering,” in Advances in Neural Information Processing Systems, 2004, pp. 1601–1608.
 [36] G. Citti and A. Sarti, Eds., Neuromathematics of Vision, Lecture Notes in Morphogenesis. Springer, 2014.
 [37] J. August and S. W. Zucker, “Sketches with curvature: The curve indicator random field and Markov processes,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 25, no. 4, pp. 387–400, 2003.
 [38] M. Favali, G. Citti, and A. Sarti, “Local and global gestalt laws: A neurally based spectral approach,” Neural Computation, vol. 29, no. 2, pp. 394–422, 2017.
 [39] J. G. Daugman, “Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters,” Journal of the Optical Society of America A, vol. 2, no. 7, pp. 1160–1169, 1985.
 [40] L. M. Florack, B. M. ter Haar Romeny, J. J. Koenderink, and M. A. Viergever, “Scale and the differential structure of images,” Image and Vision Computing, vol. 10, no. 6, pp. 376–388, 1992.
 [41] T. S. Lee, “Image representation using 2D Gabor wavelets,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 18, no. 10, pp. 959–971, 1996.
 [42] R. Duits, “Perceptual organization in image analysis,” Ph.D. dissertation, Eindhoven University of Technology, Department of Biomedical Engineering, The Netherlands, 2005.
 [43] R. Duits and E. Franken, “Line enhancement and completion via linear left invariant scale spaces on se (2),” Scale space and variational methods in computer vision, pp. 795–807, 2009.
 [44] A. A. Kalitzeos, G. Y. Lip, and R. Heitmar, “Retinal vessel tortuosity measures and their applications,” Experimental Eye Research, vol. 106, pp. 40–46, 2013.
 [45] R. Annunziata, A. Kheirkhah, S. Aggarwal, P. Hamrah, and E. Trucco, “A fully automated tortuosity quantification system with application to corneal nerve fibres in confocal microscopy images,” Medical Image Analysis, vol. 32, pp. 216–232, August 2016.
 [46] R. Duits, M. Janssen, J. Hannink, and G. Sanguinetti, “Locally adaptive frames in the rototranslation group and their applications in medical imaging,” Journal of Mathematical Imaging and Vision, pp. 1–36, 2016.
 [47] E. Franken and R. Duits, “Crossingpreserving coherenceenhancing diffusion on invertible orientation scores,” International Journal of Computer Vision, vol. 85, no. 3, pp. 253–278, 2009.
 [48] S. AbbasiSureshjani, I. SmitOckeloen, J. Zhang, and B. ter Haar Romeny, “Biologicallyinspired supervised vasculature segmentation in SLO retinal fundus images,” in Image Analysis and Recognition, ser. Lecture Notes in Computer Science. Springer, 2015, vol. 9164, pp. 325–334.
 [49] J. Shi and J. Malik, “Normalized cuts and image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 22, no. 8, pp. 888–905, 2000.
 [50] P. Perona and W. Freeman, “A factorization approach to grouping,” in Computer Vision—ECCV’98. Springer, 1998, pp. 655–670.
 [51] U. Von Luxburg, “A tutorial on spectral clustering,” Statistics and Computing, vol. 17, no. 4, pp. 395–416, 2007.
 [52] S. AbbasiSureshjani, I. SmitOckeloen, E. Bekkers, B. Dashtbozorg, and B. ter Haar Romeny, “Automatic detection of vascular bifurcations and crossings in retinal images using orientation scores,” in IEEE International Symposium on Biomedical Imaging, April 2016, pp. 189–192.