Retrieving challenging vessel connections in retinal images by line co-occurrence statistics

Natural images contain often curvilinear structures, which might be disconnected, or partly occluded. Recovering the missing connection of disconnected structures is an open issue and needs appropriate geometric reasoning. We propose to find line co-occurrence statistics from the centerlines of blood vessels in retinal images and show its remarkable similarity to a well-known probabilistic model for the connectivity pattern in the primary visual cortex. Furthermore, the probabilistic model is trained from the data via statistics and used for automated grouping of interrupted vessels in a spectral clustering based approach. Several challenging image patches are investigated around junction points, where successful results indicate the perfect match of the trained model to the profiles of blood vessels in retinal images. Also, comparisons among several statistical models obtained from different datasets reveal their high similarity, i.e., they are independent of the dataset. On top of that the best approximation of the statistical model with the symmetrized extension of the probabilistic model on the projective line bundle is found with a least square error smaller than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\%$$\end{document}2%. Apparently, the direction process on the projective line bundle is a good continuation model for vessels in retinal images.


Introduction
Tracking curvilinear structures Tree like structures such as the retinal vasculature, corneal nerve fibers, and roads from aerial photographs for cartography are widely studied both in quantitative computer-aided diagnosis systems in large-scale screening programs, and high-volume industrial settings.Delineation of curvilinear structures in these images is essential for investigating their characteristics.For instance, several studies highlighted the importance of using quantitative measurements of morphological and geometrical properties of blood vessels in retinal images for early diagnosis and prognosis of several diseases such as hypertension and diabetic retinopathy (e.g.Chapman et al, 2002;Smith et al, 2004).Despite all improvements in the segmentation of curvilinear structures in two-dimensional images, the proposed methods often present limitations when two or more structures branch or cross, or when there are areas with missing information or interruptions (Fraz et al, 2012).Consequently, several tracking-based techniques provided solutions for preserving the connections in tree-shaped networks (e.g.Türetken et al, 2012;Cheng et al, 2014;Bekkers et al, 2014;Estrada et al, 2015;Hu et al, 2015;De et al, 2016).One of the common approaches has been to manually design cost functions, which penalized abrupt changes of the contextual features such as orientation, width and color.These costs were used in later stages in optimization or graph theory based techniques for constructing the full retinal vasculature network.In these methods, not only the cost functions were designed manually and depended on existing topological structures, but also tracing errors were often created due to the use of imperfect pixelbased vessel segmentations and skeletons that do not guarantee the connections among pixels belonging to one vessel.
Geometry of primary visual cortex The human visual system is capable of interpreting visual scenes and of completing disconnected contours among interrupted segments, following the Gestalt law of good continuation (Wertheimer, 1938).Fig. 1 represents a sample interrupted phantom image (Fig. 1a) and the units detected by our perceptual system in different colors (Fig. 1b).Inspired by this capability a new method was proposed by Favali et al (2016) for resolving the missing and complex connections among blood vessels at junction points in retinal images.In this approach, first the image is lifted to the coupled space of positions and orientations using the so-called orientation score transformation introduced by Duits et al (2007).
The multi-orientation score is augmented with a contextual affinity matrix inspired by the long-range contextual connections between the multi-orientation pinwheel structures discovered in the primary visual cortex (V1) (Hubel and Wiesel, 1962;Bosking et al, 1997).
The affinities in this matrix are augmented by another similarity measurement of the feature of intensity and further processed in a spectral clustering step, which resulted in separate groups each representing one individual blood vessel.The cortical connectivity representing the contextual connections in V1 can be modeled as the fundamental solution of the time-independent Fokker-Planck (FP) equation for Mumford's direction process (Mumford, 1994;Williams and Jacobs, 1997;August and Zucker, 2003;Citti and Sarti, 2006).Another closely related model for perceptual grouping of local orientations is a hypo-elliptic Brownian motion, whose Fokker-Planck equation describes hypo-elliptic diffusion without convection for the generator (Citti and Sarti, 2006;Duits and Franken, 2009;Agrachev et al, 2009).These models explain well the notion of association fields introduced by Field et al (1993) as a model for contour integration by the human visual system for image perception.There exist various numerical approximations, and exact solutions.For a recent overview and detailed comparison of all the solutions, see Zhang et al (2016b) and references therein.Based on this study, the Fourier based technique (Duits and van Almsick, 2008) is the best approximation of the exact solution, having the smallest error both in the spatial and Fourier domains.The second best approximation is provided by the stochastic method based on the Monte Carlo simulation (Robert and Casella, 2005).The stochastic solution was used by Favali et al (2016).
Edge statistics in natural images Second and higher order edge statistics are commonly used for representing the mutual relation between connected edges in natural images.Several studies (e.g.August and Zucker, 2000;Geisler, 2008;Sanguinetti et al, 2010;Perrinet and Bednar, 2015) investigated the statistics of the edges in our surrounding environment and their relation to the adaptation of connectivity patterns in our perceptual system.These studies showed that the individual edges are dependent on each other, and the strongest characteristic determining their connections is their co-circularity relation.In the work by August and Zucker (2000), the edge statistics for several types of image categories were measured and it was shown that the resulting patterns were dependent on the structures available in the images.In a recent study by Perrinet and Bednar (2015), these low-level edge statistics were used for a high level judgmental task successfully.Moreover, Sanguinetti et al (2010) showed that there is a close relation between the statistics of edge cooccurrence and the probabilistic, geometric model of the cortical connectivity in V1.
Our proposed method We propose to train the probabilistic connectivity kernel using the line cooccurrence statistics on the lifted space of positions and (π−periodic) orientations (R 2 × P 1 ) extracted directly from the retinal images.To this end, we make an adaptation of the direction process to the projective line bundles in R 2 × P 1 instead of R 2 × S 1 (the space of positions and 2π−periodic orientations), as this extension is necessary for comparison of the probabilistic model to the statistical co-occurrences.By comparing the statistical kernel to the symmetrized probability kernel, its best approximation resulting in the least square error is found.In fact, we show the relation between the probabilistic model of cortical connectivity and the edge statistics in our retinal imaging application is even closer when including both symmetrization and a projective line bundle structure.The dependency of the parameters of the statistical kernel with respect to the dataset is also investigated using different retinal image datasets, varying their resolutions and pixel sizes.
Finally, we show the application of this trained model in retrieving the vessel connections at locations with complex structures in retinal images.To this end an affinity matrix is created based on this statistical model and the similarities among vessel intensities and is analyzed in a self-tuning spectral clustering technique (Zelnik-Manor and Perona, 2004), which does not need any parameter tuning and manual thresholding of the eigenvalues.It automatically determines the number of salient groups in the image by rotating the eigenvectors to create the maximally sparse representation and by minimizing a clustering cost defined accordingly.
Summarizing, we demonstrate the following points in this article: (a) the statistical line co-occurrence kernel learned from retinal images matches remarkably well our symmetrized extension of the probabilistic model on the projective line bundle; (b) the statistical kernels do not change significantly over different datasets and are reproducible; (c) the low-level line statistics are successfully used to perform the high level task of grouping of the interrupted blood vessels in the retinal images automatically; (d) Mumford's direction process is a very good stochastic model for connecting interrupted vessels in segmented retinal images.
Paper structure The rest of the article is structured as follows.In Sect.2, the steps for deriving the line co-occurrences from retinal images and the theoretical details about modeling the cortical connectivity are described.In Sect.3, after introducing the datasets, the resulting line co-occurrences are presented and compared against each other quantitatively and qualitatively.The best probability model approximating each kernel is presented afterwards.Application of the statistical kernel in retinal image analysis is presented at the end of the section.Finally, the results are discussed and the paper is concluded in Sect. 4.

Methodology
In this section, in addition to introducing the steps for extracting the line co-occurrences from retinal images, a numerical model of the connectivity kernel is proposed.

Line co-occurrence
In order to extract the line co-occurrences from retinal images, a similar approach as the method of Sanguinetti et al ( 2010) is used.However, there are some differences.We only use retinal images, which include multiple elongated structures: the vessels; the vessel centerlines have been used instead of the edges (the resulting kernel is called line co-occurrences rather than edge cooccurrences) and no line polarity has been taken into account; the orientation score transformation has been used to find the orientation information at each point.
Recall that the projective circle P 1 is obtained from the normal circle S 1 by identifying antipodal points.The orientation score (OS) transform R 2 → R 2 × P 1 is obtained by correlating the input image f with rotated isotropic (bi-directional) cake wavelets ψ (Bekkers et al, 2014) in where R θ is the 2D counter-clockwise rotation matrix, the overline denotes the complex conjugate and denotes the correlation (Duits et al, 2007).
The proposed method for finding the line statistics of the images of a dataset S = {I 1 , I 2 , . . ., I n }, where I i ∈ R 2 is the i th retinal image, is explained step by step in Algorithm 1.The initial step is to create a set of interest vessel positions and orientations for each image.To obtain the vessel pixel locations the vessel ground truth is used, and the binary vessel centerlines (I c,i , i = 1, . . ., n) are extracted in a standard morphological thinning approach (Lam et al, 1992).So if a pixel at location (x, y) belongs to a centerline, then I c,i (x, y) = 1, otherwise I c,i (x, y) = 0.If the ground truth is not available, then the vessel segmentation is obtained using one of the state-of-the-art techniques (e.g.Zhang et al, 2016a).The orientations at interest centerline positions are obtained by lifting the image using the OS transform (Step 2), and finding the angles with the maximum response at these locations (θ mi (x) in Step 3).It is worth mentioning that only the real part of the OS has been considered which acts as a ridge detector on the Gaussian profiles of blood vessels.Besides, the blood vessels in retinal images are darker compared to the background.As a result, they get negative responses (large absolute values) in this transformation.The negative sign used at Step 3 compensates for that.Using these locations and orientations, the set of interest points L i is created in Step 4 for each image.
In the next step, pairs of interest points located at less than a certain distance (d) from each other are used to create a difference set S d i considering the translationinvariance property (See Step 5 and Fig. 2).In order to make the set rotation-invariant, the relative positions are rotated with respect to the relative orientations and the shift-twist difference set Q d i is created in Step 6.
(x p , y p , θ p ) (x q , y q , θ q ) x q − x p y q − y p θ q − θ p Fig. 2: A sample pair of edges with positions and orientations of (x p , y p , θ p ) and (x q , y q , θ q ).The relative positions and orientations are depicted.Adapted from (Sanguinetti et al, 2010, Fig. 1).
By counting the number of iterations of relative positions and orientations in Step 7 the statistical kernel per image is created.The statistical kernel of the entire dataset is obtained by accumulating the individual kernels of all the images in set S. Finally, the kernel is l 1 -normalizedand and called the data-driven or statistical kernel (Step 8).The final dimension of this kernel is (2d + 1) × (2d + 1) × n θ .The parameter d is selected heuristically during experiments.
In another approach, the artery/vein (AV) labels of vessels and the fact that arteries are not connected to veins are taken into account.By knowing these labels, the vessel profiles for arteries and veins are separated and their line co-occurrences are calculated individually (K stat i,A and K stat i,V ), using the similar steps as described before.Finally, the artery and vein histograms are added to each other to find the final histogram for each retinal image (K ).This is a more accurate assumption about the connections among vessel centerlines; however, it is only possible to use this setup if the AV labels are available.More details about datasets, parameter settings and results are given in Sect.3.
2.2 Cortical connectivity in R 2 × P 1 Considering Mumford's direction process in the differential structure of the sub-Riemannian SE(2) group, the fundamental solution of the FP equation represents the probability of having a contour at a certain position and orientation, starting from a reference position and orientation.In order to model the cortical connectivity kernel in R 2 × P 1 (projective line bundle), we propose to create the connectivity kernel by adding the solutions of the FP equation in forward and backward directions in R 2 × S 1 (for symmetrization) and the π−shifted solutions (for taking into account the final periodicity).Therefore, the connection probability between two points in R 2 × P 1 is obtained by: where Γ is the fundamental solution of the time integrated FP equation centered around (x , θ ) represented as: with the resolvent kernel R D α obtained by integrating Green's function K D t : SE(2) → R + as: which is the solution of the following PDE: Note that this PDE is defined on R 2 × S 1 and not on R 2 × P 1 as the first order part flips when applying θ → θ + π.Therefore, in Eq. 2, besides a π−shift, we need inversion invariance yielding a double symmetric kernel (see Fig. 3).In this work the numerical solution has been created using the Fourier-based technique (Duits and van Almsick, 2008), because it is not only the best approximation to the exact solution, but also computationally the least expensive one compared to the other solutions (Zhang et al, 2016b).Referring to (Fig. 13 Zhang et al, 2016b), the slight spatial blurring with 0 < s 1 Algorithm 1 The proposed steps for obtaining the line co-occurrence statistics from a set S of retinal images 1: Obtain the binary vessel centerlines (I c,i ) from the vessel ground truth of each image I i ∈ S, (i = 1, . . ., n) by a standard morphological thinning approach.2: Lift the original image (I i ∈ R 2 , ∀i = 1, . . ., n) to the rototranslation group (U I i (x, θ) ∈ R 2 × P 1 ) using isotropic cake wavelets rotated in n θ directions in Eq. 1 so that θ ∈ {−π/2, . . ., −π/2 + (π(n θ − 1)/n θ )}.3: Find the orientation with the highest real value of negative orientation score at each position x of I i , (i = 1, . . ., n) as: 4: Create a set of interest points for each image defined as where θm i (x) is the dominant orientation at position x.5: Create a difference set defines as: 6: Create the shift-twist difference set as: 7: Obtain K stat i so that 8: Calculate the total statistical kernel for the entire dataset as: and normalize it as: . corresponds to a one-pixel bin used in statistical kernel.So the exact probability kernel (with small s) is considered with the same resolution as the statistical kernel.
The final probability kernel and the statistical datadriven kernels are compared in the spatial domain.The key parameters in creating the probability kernel are α and D 33 = σ 2 /2, which determine the expected life time of the resolvent kernel (E(T ) = 1/α) and the diffusion matrix (D = diag{0, 0, D 33 } ) respectively (see Zhang et al (2016b) for more details).To have uniform notations in the rest of the article, we assume the following relations hold for both probabilistic and statistical kernels: where g, h ∈ L i for i = 1, . . ., n.
As a side note, the exact non-symmetrized kernel on R 2 P 1 is given by: where RD,a,∞ α is expressed in two Mathieu functions in (Eq.5.5, Zhang et al, 2016b).Alternatively, one may take the simplified exact solution on R 2 × S 1 (Eq.5.11 Zhang et al, 2016b) and consider only the first and third terms in Eq. 2.

Experiments
In this section, first the datasets and the parameters used for finding the data-driven kernels are introduced.Then the results of the comparison of the data-driven kernels against each other, and comparison of these statistical kernels with the probability kernels are presented.At the end, the application of the data-driven kernel in retrieving vessel connections is explained.

Material
Two retinal datasets have been used in this study.The public DRIVE dataset (Staal et al, 2004) including 40 color fondus images taken with a Canon CR5 nonmydriatic 3CCD camera, with a resolution of 565×584,   Two different visualizations of the final statistical kernels are shown in Fig. 5 and 6.The rows in Fig. 5 from top to bottom represent the K stat DR , K stat DR−AV , K stat IO , and K stat IO−AV respectively.Each square has the dimension of 131 × 131 and it depicts the value of the kernel at fixed relative orientation (K stat (x, θ c ), θ c ∈ {±π/8, ±π/16, 0}).The kernel values of only five orientations are depicted as the information at other angles is very small.As seen in these figures, the maximum values of the statistical kernels occur at small orientation differences i.e., two aligned lines are more probable to appear in the images.This probability decreases when the orientation differences increase.Comparing these four histograms qualitatively, the statistical kernels of the DRIVE dataset are a bit less elongated compared to the ones of the IOSTAR datasets.Moreover, the separation of the lines of arteries and veins from each other results in less noisy histograms for both datasets; however, the difference is very small.Fig. 6 visualizes the isosurfaces of these four data-driven kernels, which shows their high similarity.
These data-driven kernels are compared with each other quantitatively and their mutual differences are are the two statistical kernels in comparison.These least square errors are presented in Table 1 for each pair of kernels.Based on these quantitative results, their differences are very small.Therefore, it is possible to use them interchangeably, regardless of the dataset or the ground truth used for obtaining them.

Comparison to the probability model
In this section we find the best approximations of the statistical kernels by comparing them against various Table 1: The least square errors obtained by comparing each pair of the statistical kernels probability kernels (obtained using Eq. 2) with different parameters.The kernels that result in the least square errors are considered as the best approximations.In our experiments, the parameter α takes 50 different values from 0.00001 to 0.01 and the parameter D 33 takes 100 values from 0.000001 to 0.005.These ranges are determined heuristically.The minimum errors and the corresponding parameters for each kernel are presented in Table 2. Based on these results, the errors are very small for all the kernels.The largest error is related to the K stat IO−AV , which is close to 1%.A sample qualitative comparison of these two types of kernels is shown in Fig. 7.In addition to very similar profiles of the two kernels in Fig. 7a, b and d, the summations of the line distribution over the spatial dimension also matches the information density of the probability kernel at every θ layer (Fig. 7c).In both kernels, the maximum value appears at θ = 0 as expected.

Application in retinal image analysis
In this section, the application of the statistical model of the cortical connectivity pattern in identifying vessel connections in retinal images is explained.
In this method the vessel connections are retrieved from image segmentations (not necessarily centerlines).So an initial segmentation of the image I (using the aforementioned segmentation techniques proposed in the literature) is required.The binary image representing the segmentation is called I s .Repeating the steps as the ones proposed in Algorithm 1 (Steps 2, 3 and 4) the image is lifted (U I ), dominant orientations (θ m ) are obtained and the set of interest points is created as: Considering m = |L|, an m × m affinity matrix (A) is created to take into account the connection probability between each pair of points in set L and is later aggregated with another affinity matrix ( Ã) representing the similarities of the corresponding vessel intensities.So for each pair of (x i , θ m (x i )) and (x j , θ m (x j )) ∈ L: where G σint (x) is the normalized Gaussian kernel with the standard deviation of σ int and I n is the image intensity in normalized green channel.Here we normalized luminosity and contrast using the method by Foracchia et al (2005).The final affinity matrix is analyzed using the self-tuning spectral clustering technique, which identifies the salient groups in the image automatically.
As discussed before, the differences between the kernels obtained from the same datasets are very small and they can be used interchangeably.Hence, for analyzing the retinal image patches taken from the DRIVE and IOSTAR datasets we only use the K stat DR−AV and the K stat IO−AV kernels respectively.For testing the method, several patches with the sizes of 51 × 51 and 101 × 101 have been selected around junctions from the DRIVE and IOSTAR datasets respectively.A vessel segmentation by Abbasi-Sureshjani et al ( 2015) developed for both color RGB images and SLO images is applied on these patches.After segmenting the images, the positions, orientations and normalized intensities for the vessel pixels are extracted.
Fig. 8 shows five image patches for each dataset.The first five patches (D 1 to D 5 ) belong to the DRIVE and the second five patches (I 1 to I 5 ) are selected from the IOSTAR datasets.For each patch, from top to bottom the vessel intensities, the color-coded orientation maps (representing the values of θ m (x)), the AV labels, and finally the clustering results have been presented.As seen in this figure, the detected final perceptual units (shown in different colors) correspond to the separate blood vessels existing in the image patch.These patches have been selected in a way to represent several complex topological structures of the blood vessels, varying the number of the vessel bifurcations and crossings and existence of parallel or curved vessels.For all the experiments the parameter σ int was set to 0.2.Despite the complexity of the structures, using orientation and intensity features for determining the contextual connections among pixels, individuates well the blood vessels from each other.Using the feature of intensity helps in the cases where there is an abrupt change in the orientation but not the intensity e.g. in D 3 and D 4 .Moreover, the presence of disconnections (e.g. in I 5 ) does not affect the correct grouping.In addition to these patches, 20 more patches per dataset have been analyzed.For all these cases, the method is successful in correctly grouping the vessel pixels.The limitation arises when both the feature of intensity and orientation of a vessel are very noisy or change suddenly.In these cases, the vessel splits into smaller clusters.Involving additional contextual information such as curvature or scale may help in resolving this problem as proposed by Abbasi-Sureshjani et al (2016).

Conclusion
In this work, we exploit the relation between the statistical co-occurrences of line elements in natural images and the high level task of contour grouping in our brain, and use it in the retinal image processing for contour completion.Firstly, the steps for obtaining the rotation and translation-invariant line statistics from different retinal image datasets are explained.Secondly, their relation to the symmetrized probability kernel of the direction process on the projective line bundle modeling the cortical connectivity is investigated.The results reveal their remarkable high similarity.However, in edge co-occurrences (rather than line co-occurrences) of natural images (e.g.Sanguinetti et al, 2010) the shapes seem to resemble the shape of the hypo-elliptic heat kernels (without convection) in R 2 × P 1 , but the relation needs further investigation.In addition, all the statistical kernels obtained from different datasets are compared with each other quantitatively and qualitatively.The obtained results indicate their high similarity and reproducibility despite the differences in the datasets and the setups used.
Furthermore, the statistical kernels were used directly to retrieve the individual vessels from segmented retinal images.The successful results show that Mumford's direction process is a very good model for centerlines of vessels, and together with the Lie group theory, the proposed connectivity analysis technique is useful for retinal image analysis.Using the data-driven kernel (that does not need parameter tuning) in addition to adding the automatic self-tuning spectral clustering technique, forms a robust and fully automatic connectivity analysis technique.This provides an effective solution for challenging situations in which most of the methods fail (but our visual perception succeeds) because of non-perfect imaging conditions, interruptions, or occlusions.
The method presented here for extracting line cooccurrences from retinal images and the improved connectivity analysis approach can be extended to a rich number of other application areas which contain curvilinear structures such as corneal never fibers, plant roots and road networks.For each dataset, it is possible to learn the connectivity kernel and use it directly for similar images.
Since the visual cortex deploys additional contextual information and its receptive fields are not only sensitive to orientation, but also other information (such as scale and curvature), one potential extension of the method is to use this additional information in deriving the line statistics and creating higher order kernels.

Fig. 1 :
Fig. 1: (a) A sample image with interrupted segments and (b) the salient units identified by our perceptual system

Fig. 4 :
Fig. 4: Two sample images from the IOSTAR (top row) and the DRIVE (bottom row) datasets.(a) The original images, (b) the vessel and (c) the AV ground truth images (artery in red and vein in blue), (d) the vessel skeleton and (e) the color-coded orientation maps

Fig. 5 :Fig. 6 :
Fig. 5: The statistical kernels obtained in Step 8 of Algorithm 1 for each dataset.The value of θ is shown for each column and the values of all figures are clipped between 0 to 0.2 of the maximum value of the K stat DR

Fig. 7 :
Fig. 7: Comparison between the K stat DR and its best approximating kernel K prob DR : the level sets at 0.012 of the maximum value of (a) the statistical and (b) the probability kernels, (c) the summation of the values of the two kernels over the spatial dimension at different angles, and (d) the line distribution at five different orientations θ ∈ {±π/8, ±π/16, 0} in the data-driven (top row) and the probability kernel (bottom row) Fig. 8: Some example retinal patches selected from the DRIVE (D 1 to D 5 ) and the IOSTAR (I 1 to I 5 ).The first three rows from top to bottom indicate the intensity, main orientations and the AV labels of the vessels.The last row represents the final clustering results after removing the small groups (noise).Detected clusters are shown in different colors

Table 2 :
The least square errors and the corresponding parameters resulting in these errors between the statistical and probabaility kernels