A Multi-Orientation Analysis Approach to Retinal Vessel Tracking

This paper presents a method for retinal vasculature extraction based on biologically inspired multi-orientation analysis. We apply multi-orientation analysis via so-called invertible orientation scores, modeling the cortical columns in the visual system of higher mammals. This allows us to generically deal with many hitherto complex problems inherent to vessel tracking, such as crossings, bifurcations, parallel vessels, vessels of varying widths and vessels with high curvature. Our approach applies tracking in invertible orientation scores via a novel geometrical principle for curve optimization in the Euclidean motion group SE(2). The method runs fully automatically and provides a detailed model of the retinal vasculature, which is crucial as a sound basis for further quantitative analysis of the retina, especially in screening applications.

only suitable for investigation of ocular diseases such as glaucoma and age-related macular degeneration (AMD), but also for systemic diseases such as diabetes, hypertension and arteriosclerosis. This makes the retinal vasculature a rewarding and well researched subject and and a growing number of image processing techniques are developed to segment and analyze the retinal vasculature [1,47].
Retinal vessel tracking Typically there are two types of methods for vessel extraction: pixel classification methods [42,49,13,46] and vessel tracking methods [2,14,65,37,25,50,16]. The first type of method classifies pixels as either being part of a vessel or background, resulting in a pixel map in which white pixels represent blood vessels. Of the pixel classification methods, the approach by Krause et al. [42] is most similar to ours as both methods rely on a transformation to a higher dimensional domain. In their work they applied vessel detection based on the local Radon transform, of which we will show later in this article that this is a special case of an orientation score transform based on cake wavelets. The other type of method, vessel tracking, is based on recursively expanding a model of the vasculature from a set of seed points. One advantage of vessel tracking over pixel classification is that it guarantees connectedness of vessel segments, whereas in pixel classification methods this is not necessarily the case. For further quantitative analysis of the vasculature, tracking algorithms are preferred because they intrinsically provide geometrical and topological information. For example, vessel widths, curvatures, segment lengths, bifurcation density and other features can relatively easily be extracted from the generated vessel models.
Several different approaches to vessel tracking can be found in literature. There are methods based on active contours [2,25], matched filters [14,37,16], and probabilistic models [65] among others [65,50]. The majority of papers on vessel tracking report limitations regarding tracking blood vessels through crossings, bifurcations and/or more complex situations. In this paper we aim at solving these problems by means of orientation analysis via so called orientation scores, which are objects in which image information is neatly organized based on position and orientation [17]. We propose two new tracking algorithms that act directly on the domain of an orientation score, and we show that these methods are highly capable of dealing with the aforementioned problems. Afterwards, we will extend one of the orientation score based algorithms to a vasculature tracking algorithm, which is capable of constructing models of the complete retinal vasculature.
Orientation scores Inspired by the cortical orientation columns in the primary visual cortext [39], Duits et al. developed a mathematical framework for image processing and analysis based on 2D orientation scores [20]. Similar to the perceptual organization of orientation in the visual cortex, a 2D orientation score is an object that maps 2D positions and orientation an-gles (x, y, θ) to complex scalars. Instead of assigning an orientation to each position and thereby extending the codomain, we extend the domain of the image where our modeling can deal with multiple orientations per position. When constructing an orientation score it is crucial one does not tamper the data evidence before tracking takes place. Therefore we consider invertible orientation scores that provide a full comprehensive overview of how the image is decomposed out of local (multiple) orientations. In invertible orientation scores, all orientated structures are disentangled, see Fig. 1.
Paper structure The article is structured as follows: First, in Section 2, theory about orientation scores is provided. In Section 3, two vessel tracking approaches based on orientation scores are described: the ETOS-algorithm: an all-scale approach based on a new class of wavelets, the so-called cake wavelets the CTOS-algorithm: a multi-scale approach based on the classical Gabor wavelets Both tracking methods rely on a novel generic geometrical principle for curve optimization within the Euclidean motion group, which is explained and mathematically underpinned in Appendix A. We will show that ETOS generally works with different types of orientation scores, however with best performance on invertible orientation scores based on cake wavelets [17,22] (in comparison to non-invertible orientation scores based on Gabor wavelets). The second approach requires a multi-scale and orientation decomposition. The two approaches are described in Section 3.1, and evaluated in Section 3.2. It will turn out that ETOS based on cake wavelets has several advantages over CTOS based on Gabor wavelets. We have validated ETOS more extensively by comparing it to the state of the art in retinal vessel tracking [2,6,64] using the publicly availabe REVIEW database [2]. In Section 4.1 we describe our vasculature tracking algorithm, composed of proper initialization, junction detection and junction resolver algorithms. In Section 4.2 the correctness of topology of the models is evaluated using images of the HRFIdatabase [12]. General conclusions can be found in Section 5.

Orientation scores
Orientation detection and encoding is a common subject in image processing. E.g., Frangi et al. [30] detect retinal blood vessels by calculating a vesselness value, obtained by eigenvalue analysis of the Hessian matrix.
In [55] this is done by taking the maximum modulus over a set of oriented Gabor wavelet responses. Besides the presence of local orientations, the orientation value(s) may be relevant as well. For instance in vessel tracking and certain other image enhancement techniques like coherence enhanced diffusion [61] and orientation channel smoothing [26] The most commonly used methods to detect orientations are capable of detecting only one orientation per position. However, by using oriented wavelets and steerable filters [33,53,48,35] orientation confidence measures can be extracted for any given orientation, thus allowing for the detection of multiple orientations per position. Oriented wavelets allow for a transformation from an image to an orientation score, where each locally present combination of position and orientation is mapped to a single value [41,20,22], see Fig. 1.
In his pioneering paper, Kalitzin [41] proposed a specific wavelet, given by that, by approximation, guaranteed invertibility from the orientation scores back to the original image, without loss of information. These wavelets belong to a specific class of so called proper wavelets (wavelets that allow well-posed reconstruction [20,22]), and are found by expansion in eigenfunctions of the harmonic oscillator. The advantage of this expansion is that this steerable basis is Fourier invariant, allowing to control the wavelet shape in both the spatial and Fourier domain.
The disadvantages of such kernels ψ are however that 1) their series do not converge in L 2 and truncation of the pointwise converging series heavily affects the shape and induces undesirable oscillations [17, p.140-142]; 2) the wavelets ψ explode in the radial direction along its orientation, e.g. for the case in (1) it explodes with 3) does not allow approximate reconstruction by integration over S 1 only.
In Kalitzin's paper the well-posedness of the reconstruction was not quantified. Analysis of the wellposedness was done by Duits using the function M ψ , which will be explained in more detail in Section 2.1. The function M ψ indicates how well spatial frequencies in an image are preserved after a transformation, and it can be seen as a measure for stability of the inverse transformation [22,27,17,20]. Within these papers a new class of proper wavelets called cake wavelets are presented. These are oriented wavelets, able to capture all image scales without any bias to a specific scale. This property is crucial in our vessel edge tracking approach since it reduces the need for a multi-scale approach, as will become more clear in the next sections.
In contrast, the Gabor wavelet is another oriented wavelet, which is widely used in the field of image processing because of its capability to detect oriented features at a certain scale. An example of segmenting crossing lines using Gabor wavelet based orientation scores is given in [15]. The property to tune the Gabor wavelet to capture features at a specific scale can be very useful, but the scale selection also implies exclusion of other scales. The single-scale Gabor wavelet transformation causes information from the original image to be lost and the transformation is therefore non-invertible. In order to introduce invertibility, one has to use a multiscale approach [28], which is also computationally more expensive.
Since no information should be lost during the orientation score transformation W ψ : f → U f (see Fig. 1e), the notion of invertibility of an orientation score transformation is essential. The invertibility allows us to relate operators on images to operators on orientation scores, and vice versa. Using invertible orientation scores, one can employ the automatic disentanglement of local orientations involved in a crossing (cf. Fig. 1f). For line/contour enhancement, this has lead to a generic crossing preserving diffusion method [32,21,31] which outperformed related diffusions acting directly on the image domain. For tracking of crossing blood vessels a similar advantage can be employed, as we will show in this article.

Construction of orientation scores
Consider a 2D image f as a function f : R 2 → R, with compact support on the image domain Ω = [0, X] × [0, Y ], with image dimensions X, Y ∈ R + , and which we assume to be square integrable, i.e. f ∈ L 2 (R 2 ). An orientation score, constructed from image f , is defined as a function U f : R 2 × S 1 → C and depends on two variables (x, θ), where x = (x 1 , x 2 ) ∈ R 2 denotes position and θ ∈ [0, 2π] denotes the orientation variable.
An orientation score U f := W ψ f of a function f can be constructed by means of convolution with some anisotropic wavelet ψ via where ψ ∈ L 1 (R 2 ) L 2 (R 2 ) is the convolution kernel with orientation θ = 0, i.e. aligned with the vertical axis in our convention, and where W ψ denotes the transformation between image f and orientation score U f . The overline denotes complex conjugate,ψ θ (x) = ψ θ (−x) and the rotation matrix R θ is given by see Fig. 1e. Exact reconstruction 1 from the orientation scores constructed by (2) is given by where F is the unitary Fourier transform on R 2 , where W * ψ denotes the adjoint wavelet transformation (see [17] for details), and M ψ : The function M ψ provides a stability measure of the inverse transformation. Theoretically, reconstruction is well-posed, as long as where δ is arbitrarily small, since then the condition number of W ψ is bounded by M δ −1 , see [17,Thm. 20].
If we do not restrict ourselves to band-limited/disklimited functions, this requirement (Eq. (6)) bites the assumption ψ ∈ L 1 (R 2 ) L 2 (R 2 ) since it implies Fψ is a continuous function vanishing at infinity (see e.g. [52]) and so is M ψ . In that case we have to resort to distributional wavelet transforms 2 whose closure is again a unitary map from L 2 (R 2 ) into a reproducing kernel subspace of L 2 (R 2 × S 1 ), for details see Appendix B.
In practice, to prevent numerical problems, it is best to aim at M ψ (ω) ≈ 1 for ω < , where is the Nyquist frequency of the discretely sampled image, meaning that all relevant frequency components within a ball of radius are preserved in the same way. Because of the discontinuity at ω = , which causes practical problems with the discrete inverse Fourier transform, we will use wavelets ψ, where t denotes a scale parameter. To fix the inflection point close to the Nyquist frequency, say at ρ = γ with 0 γ < 1, we set t = 2(γ ) 2 1 + 2N (to fix the bending point: Fig. 2). The function M N basically is a Gaussian function at scale t, multiplied with the Taylor series of its inverse up to a finite order 2N to ensure a slower decay. The function M N smoothly approximates 1 on the domain ρ ∈ [0, ], see Fig. 2. A wavelet ψ : R 2 → C with such a M ψ will be called a proper wavelet. The methods presented in this paper are best described by a moving frame of reference that lives in the tangent bundle (T (R 2 × S 1 )) above the domain R 2 × S 1 of an orientation score. This moving frame of reference is given by the map with e x = (1, 0, 0) and e y = (0, 1, 0), e θ = (0, 0, 1) tangent vectors (formally attached at g). To simplify the notation we introduce coordinates ξ = x cos(θ) + y sin(θ), η = −x sin(θ) + y cos(θ) and tangent vectors e ξ = cos(θ)e x + sin(θ)e y = (cos(θ), sin(θ), 0) e η = − sin(θ)e x + cos(θ)e y = (− sin(θ), cos(θ), 0) e θ = (0, 0, 1).
As a result we have that at a given point in the orientation score (x, y, θ), the tangent vector e ξ points in the spatial direction aligned with the orientation score kernel used at layer θ, see Fig. 1b and 1d. We will often rely on the notation in Eq. (9) in the remainder of this article.

The domain of orientation scores: SE(2)
The domain of an orientation score is the set R 2 × S 1 . However, from Fig. 1d one can recognize a curved geometry on the domain of orientation scores. This is reflected in the fact that e ξ and e η vary with θ. This is modeled by imposing a group structure on the set R 2 × S 1 . This group structure comes from rigid body Note that (x, θ) = (x, R θ )(0, 0) which allows us to uniquely identify i.e. to identify the space of positions and orientations with the rigid body motion group SE(2) = R 2 SO(2).
As the combination of two rigid body motions is again a rigid body motion, SE(2) is equipped with the group product: which is consistent with Eq. (10). The moving frame of reference (9) corresponds to the so-called left-invariant vector fileds in SE (2). For details see [31, Fig. 2.5a].

Cake wavelets
Cake wavelets are constructed from the Fourier domain. By using polar coordinates, the Fourier domain can be uniformly divided into N o equally wide samples ("pieces of cake") in the angular direction, see Fig. 3. The spatial wavelet is given by where G σs is a Gaussian window, with 0 < 1 σ s , that is used to avoid long tails in the spatial domain. Note that multiplication with a large window in the spatial domain corresponds to a convolution with a small window in the Fourier domain, such that M ψ is hardly affected with σ s sufficiently large. Functionψ cake is given bỹ with ω = (ρ cos ϕ, ρ sin ϕ) and where s θ = 2πN −1 o is the angular resolution in radians. The function M n specifies the radial function in the Fourier domain given by (7). B k denotes the kth order B-spline given by Orientation scores constructed from an image f using cake wavelets are denoted by U cake f . Fig. 4c demonstrates a typical cake wavelet based orientation score for a certain orientation.
The approach of constructing wavelets directly from the Fourier domain allows indirect control over the spatial shape of the filter, and it can easily be adapted. For example, the number of orientations N o specifies the angular resolution s θ : If N o is large, the resolution in the orientation dimension is large and the filters become very narrow. This is illustrated in Fig. 5. The cut-off frequency (at the inflection point) of the function M n , which is usually set as the Nyquist-frequency, could be lowered to filter out high-frequency noise components. Moreover, because B-splines and the function M n are used to sample the Fourier domain, the sum of all cake wavelets is approximately 1 over the entire Fourier domain (within a ball of radius γ ), see Fig. 3 and 5. Thus the cake wavelets indeed are proper wavelets, allowing stable reconstruction via Eq. (4). In Eq. (4)   Remark If N o → ∞ thenψ cake → δ 0 in the distributional sense. In this case the wavelet transform converges to the localized Radon transform, which has been proposed for effective retinal vessel detection in [42]. The advantage however of taking N 0 ∞ is that we obtain well-posed non singular kernels in the spatial domain while allowing a stable reconstruction for all apriori set N o ∈ N.
Cake wavelets are quadrature filters, meaning that the real part contains information about the locally even (symmetric) structures, e.g. ridges, and the imaginary part contains information about the locally odd (antisymmetric) structures, e.g. edges. That is, the real and imaginary part of the filter ψ θ are related to each other by the Hilbert transform in the direction perpendicular to the wavelets orientation, which is defined by (15) where e η specifies the direction in which the Hilbert transform is performed, recall Eq. (9).
The quadrature property is useful in our vessel tracking approach, since it allows us to directly detect vessel edges from the imaginary part of the orientation scores, without having to calculate first-order derivatives perpendicular to the vessel orientation. In our applications we did remove the DC-component for the real part to avoid responses on locally constant images. Fig. 5 shows the real and imaginary part of the cake wavelet in the spatial domain, as well as the coverage of the wavelet in the Fourier domain.
A final remark on cake wavelets: Since the cake wavelets uniformly cover the Fourier domain (M N ≈ 1), they allow us to use a fast approximate reconstruction scheme, which is given by integration of the orientation scores over the angles only: for details see [17].

Gabor wavelets
Gabor wavelets are directional wavelets, and can be tuned to specific spatial frequencies (and inherently scales).
In the field of retinal image processing they are used for vessel detection in various studies [55,43]. We can exploit the tuning of the wavelet to specific spatial frequencies to match differently sized blood vessels. The 2D Gabor wavelet is a Gaussian, modulated by a complex exponential, and is defined as: where A = diag[ −1/2 , 1] with ≥ 1 is a 2 × 2 diagonal matrix that defines the anisotropy of the wavelet. The vector k 0 defines the spatial frequency of the complex exponential and C ψ normalizes the wavelet to unity. In our method we use = 4, which makes the filter elongated in the x-direction and we choose k 0 = (0, 3), which causes oscillations perpendicular to the orientation of the wavelet. We can dilate the filter by a scaling parameter a > 0: Orientation scores constructed from an image f using Gabor wavelets at scale a are denoted by U Gabor f,a . Fig. 4d-e demonstrate typical Gabor wavelet based orientation scores at a small and large scale respectively, for a certain orientation.
The real and imaginary part of the wavelet in the spatial domain, as well as the coverage of the wavelet in the Fourier domain are shown in Fig. 5. Similar to the cake wavelets, Gabor wavelets also have the quadrature property orthogonal to their orientation.
In the Fourier domain, the Gabor filters are represented as Gaussian functions shifted from the origin by k 0 . The set of all rotated Gabor functions at a certain scale a covers therefore a certain annulus in the Fourier domain, and a single Gabor wavelet can thus be regarded as an oriented band-pass filter. This is also clearly depicted by the outlines of the frequency responses as shown in Fig. 5.

Double-sided vs single-sided wavelets, orientation vs direction
The cake and Gabor wavelets are double-sided wavelets which do not distinguish between a forward or backward direction (they are symmetric with respect to the y-axis). In order to distinguish between π symmetries and 2π symmetries (see Fig. 6), and to be able to handle bifurcations, we decompose the orientation scores into a forward and backward direction 3 , denoted by a + and − symbol respectively: where and where with Note that by using the error function, we have ψ = . It is thus possible to choose one of the single-sided wavelets to construct a directional orientation score, while still being able to access the original (double-sided) orientation score. Fig. 6 demonstrates the advantage of using single-sided wavelets over double-sided wavelets in the case of direction estimation based on the orientation column of a score U f (x, y, ·).

Vessel tracking in orientation scores via optimization in transversal tangent planes V
In orientation scores, information from the original image is both maintained and neatly organized in different orientations, leading to the disentanglement of crossing structures. Moreover, because of the quadrature property of the wavelets used in the construction of orientation scores, important edge information is well represented in the imaginary part of the score. We therefore propose tracking algorithms that directly act on the orientation score: 1. The ETOS algorithm: Edge Tracking based on Orientation Scores (Section 3.1.1). 2. The CTOS algorithm: Centerline Tracking based on multi-scale Orientation Scores (Section 3.1.2).
(a) Vessels in OS and detection plane In both tracking algorithms we rely on a fundamental geometric principle to extract the most probable paths in orientation scores. Consider to this end Fig. 7a, where spanned by ∂ θ and ∂ η = − sin θ∂ x + cos θ∂ y in tangentbundle where T g (SE(2)) denotes the tangent space at g = (x, y, θ) ∈ SE (2). For more details on this principle we refer to Appendix A.

ETOS: Edge tracking in orientation scores
The ETOS algorithm tracks both vessel edges simultaneously through an orientation score. The method iteratively expands a blood vessel model by detecting, at each forward step k, the optimal edge locations (u k , θ k ), (v k , θ k ) ∈ SE(2) from the orientation score. Here u k and v k denote the 2D left and right edge position respectively, θ k denotes the orientation of the blood vessel. At each iteration the vessel center point c k and the vessel width w k are defined as follows: To describe our method we will rely on a moving frame of reference with basis vectors e η k , e ξ k and e θ k , which are described by the orientation parameter θ k and Eq. (9). In our method the edge positions g u k and g v k are detected in the orientation score from the tangent plane V (yellow plane in Fig. 7a). An edge can be detected as a local optimum from the imaginary part of this plane; a local minimum and maximum for the left and right edge respectively (as indicated by the two red dots in Fig. 7c). A schematic overview of the tracking process, including the symbols used in this section, is presented in Fig. 8.
For the sake of speed and simplicity, we follow a 2-step approach where the process of detecting the optimal edge positions is separated into two 1D optimization tasks which simply involve the detection of local minima (left edges) and maxima (right edge); in step 1 the edge locations u k and v k are optimized in the e η direction (η-optimization), in step 2 the corresponding orientation is optimized in the e θ direction (θ-optimization), see Fig. 7b-d. By considering the continuous properties of the blood vessels (e.g. continuous vessel widths), we use a paired edge tracking approach where the left and right edges are detected simultaneously. This approach has the advantage that even if one of the edges is less apparent in the image (e.g. at crossing points, parallel vessels and bifurcations), both edges and their orientation can still be tracked. A possible disadvantage of this approach is however that abrupt changes in vessel width (e.g. at stenoses and aneurysms) may become unnoticed or detected with less detail.
Step 1: Based on a-priori knowledge about the previous vessel orientation, edges and center point, stored in where we set M = 10, a new vessel center pointc k is calculated as where λ (typically in the order of 2 pixels) is the tracking step-size. New edge points are selected from a set of points p k (η) on a line L k going through the estimated center pointc k and perpendicular to the vessel orientation θ k−1 : with where η is a parameter describing the distance to the estimated vessel center point and η max > u k − v k denotes the maximum distance to the estimated the vessel center point. Note that orientation θ k−1 of the previous iteration is used as the new orientation is yet to be detected. An intensity profile I η k (η) can then be obtained from the orientation scores according to see Fig. 7d. New edge points can now easily be found by detecting the local optima on the imaginary part of this profile. However, the detection of vessel edges is made more robust by taking into account that the vessel wall is a continuous structure, and that the width of a blood vessel gradually changes, rather than abruptly. Therefore, we introduce the edge probability envelope. The edge probability envelope is used to indicate the most likely position of the vessel edges and it consists of two Gaussian distributions, one around the expected left vessel edge position and one around the right vessel edge position. The envelope function is given by: where η 0 is the estimated location of the vessel center on L k , w k is the mean vessel width calculated over the last M iterations: the standard deviation of the Gaussian distributions is denoted with σ, and η 0 is used to align the envelope with the actual vessel profile. A robust value for η 0 is found by optimizing the cross-correlation of the envelope with the imaginary part of the actual profile: see Fig. 9c. The left and right edges are finally detected as the arguments corresponding to the minimum and  maximum points, η k lef t and η k right respectively, of the product of the envelope and the intensity profile: see Fig. 9d for an example. The new edge points can then be assigned by Step 2: Orientation θ k can be estimated by selecting the orientation that provides the highest orientation score response at both vessel edges. The orientation score response is a combination of the orienta-tion columns at the left and right edges (U f (u k , ·) and U f (v k , ·) resp.), and the optimal orientation is calculated as: Finally the new center point c k , which may not be equal toc k , is calculated as the point between the two edges, according to Eq. (24).

CTOS: Multi-scale vessel center-line tracking in orientation scores
In this section the scale-selective property of the Gabor wavelets is exploited in the design of a fast orientation score based method called CTOS. The potential presence of a central light reflex in a blood vessel makes the design of a simple and fast centerline tracking algorithm based on local minima tracking in the image nearly impossible. However, by using Gabor wavelets and appropriate scale selection, central light reflexes can be filtered out such that vessel center points can easily be detected as local minima on detection profiles (see Fig. 10). The CTOS algorithm follows the same geometric principle as the ETOS algorithm, for a mathematical underpinning see Appendix A, however CTOS is done on the real part of orientation scores in order to find the vessel center line, rather than the vessel edges. In this algorithm each iteration k consists of 3-steps: In step 1 the center point c k is detected (η-optimization), in step 2 the orientation θ k is detected (θ-optimization) and in step 3 the vessel scale a k is detected (a-optimization).
Step 1: Using the vessel center point c k−1 , orientation θ k−1 and scale a k−1 , which were detected during iteration k−1, phase 1 at iteration k is started by estimating the new center pointc k as given by Eq. (26). The new center point c k is selected from a set of candidate points p k (η) as given by Eq. (28). From the candidate points p k (η) a center point detection profile I η,Gabor k (η) is obtained by: with U Gabor f,a k−1 the orientation score generated by the Gabor wavelets at scale a k−1 . The new center point is detected as the coordinate belonging to the local minimum on the intensity profile nearest toc k : Step 2: The new vessel orientation is detected as the local maximum of the negative orientation response, nearest to the previous vessel orientation θ k−1 : Step 3: At the new center point c k and orientation θ k , scale a k is detected as the scale that gives the largest negative response: Compared to the ETOS algorithm, this algorithm is very fast since it only requires three basic (deterministic) detection steps. However the combination of scale and orientation detection makes the algorithm slightly less stable: orientation detection depends on the correct detection of scale and vice versa.

Validation
The algorithms were tested on the green channel of color fundus images. For each image the luminosity is normalized by disposing low frequency luminosity drifts. The low frequency drifts are detected by large scale Gaussian blurring of the image (typically σ = 32), and are subtracted from the original image.

Algorithm behavior at complex vessel junctions
A qualitative validation is done using a challenging set of 4 sub-images (see Fig. 11), which were taken from the high-resolution fundus images of the HRFI database [12]. This set of sub-images contains crossings, overlapping bifurcations with crossings, small vessels crossing large vessels, small vessels, curved vessels, parallel vessels, etc. In each sub-image we manually placed seed points at the start of each blood vessel and at each bifurcation. In total 27 seed points were marked. Each seed point contains initial vessel center position, left edge position, right edge position and orientation, denoted by c 0 , u 0 , v 0 and θ 0 . The initial scale for the CTOS algorithm is detected as the scale that provides the largest scale response at c 0 and θ 0 (see Eq. (39)).
The tracking experiments are conducted using the following set of tracking parameters: The step size is set to λ = 2 pixels; The width of the scan line is set to 2η max = 40 pixels; The number of orientations used to construct the orientation scores is set to N o = 36 and the standard deviation of the Gaussian distributions used in the edge probability envelope is set to σ = 3.
The ETOS algorithm was tested on both invertible orientation scores, which were constructed by cake wavelets, and non-invertible orientation scores, which were constructed by Gabor wavelets. The scale of the Gabor wavelets was chosen in such a way that the relevant vessel features were presented as well as possible in the orientation scores (e.g. a scale too large would only represent the very large blood vessels correctly, and a scale too low only the small vessels). We found that a = 3 2π 10 gave best results. For the CTOS algorithm we used a set of orientation scores constructed by Gabor wavelets at scales a = 3τ 2π with τ =5,10,15,20,25 and 30.  Results of the tracking experiments are shown in Fig. 11. From this figure we see that, at complex situations, the ETOS method (column 2 and 3) outperforms the CTOS method (column 4). Best results are obtained when ETOS is used with invertible orientation scores (column 2). The ETOS algorithm acting on the invertible orientation scores generated by the cake kernels only fails to correctly track blood vessel nr 5 from image 3. The algorithm gives excellent results for all other vessels and manages to track the blood vessels through all complex situations, even when the contrast of the vessel edges is very low.
The performance of the ETOS algorithm is slightly decreased when applied to non-invertible orientation scores based on Gabor filters. It now fails to track 3 vessels correctly. The scale selective property of the Gabor wavelets, resulting in non-invertible orientation scores, causes the ETOS algorithm to perform less accurately compared to the application to invertible orientation scores.
The CTOS algorithm, which relies on a multi-scale orientation score approach, has the lowest performance. It fails to correctly track 5 blood vessels. In some cases, incorrect scale selection causes small parallel blood vessel to be detected as one large blood vessel (vessel 2 with 4, and 3 with 6 in the first image). Other tracks failed as a result of incorrect orientation detection.
In conclusion we can state that ETOS outperforms CTOS and that it gives best results when applied on invertible orientation scores.

Validation of width measurements
In the previous section we showed that the ETOS algorithm is very well capable of tracking blood vessels through complex situations. In this section we quantitatively validate the reliability of the measured vessel widths that are provided by the ETOS algorithm. This is done by comparing the measured widths to ground truth width measurements provided by the RE-VIEW database [2]. The REVIEW database consists of 16 color fundus images, which can be divided into 4 subsets: 1) Kick point image set (KPIS), 2) Central light reflex image set (CLRIS), 3) Vascular disease image set (VDIS) and 4) the high resolution image set (HRIS). Each image set represents images of different quality and resolution, and all are provided with manual width measurements that are performed by three individual observers. A ground truth of the vessel widths is constructed by averaging the measurements of the observers. The HRIS set contains high resolution images (3584x2438) and were down-sampled by a factor of four before submission to the ETOS algorithm. For more information on the dataset we would like to refer to [2].
When testing our algorithm, we tracked each segment by initializing the algorithm using the first pair of manually marked edge points. The same parameters that are described in Section 3.2.1 were used. Fig. 13 shows a selection of the tracking results in comparison to ground truth vessel edge labeling.
In total 5066 vessel width measurements are available. The error between automated measurements and the ground truth measurements is defined as where w i is the estimated width as measured by the ETOS algorithm (recall Eq. (25)), and w GT i is the ground truth width of the ith profile. To be able to compare our method with others we follow the same validation procedure as described in [2,6,64], where the main focus is on the standard deviation of the errors. This is motivated by the idea that different implicit definitions of vessel widths may lead to consistent errors. If this bias however is consistent enough, the error could easily be accounted for by subtraction of a bias constant. A low standard deviation of the errors indicates that the error is consistent. Fig. 12: A scatter plot, plotting 5059 ground truth widths against the widths measured by our ETOS algorithm. The linear regression model y = 0.85 + 0.88x indicates an offset of less then a pixel, suggesting that ETOS slightly overestimates the vessel widths. The slope of 0.88 indicates a strong positive relation between the ground truth and measured widths. Table 1 shows the validation results of our ETOS algorithm, in comparison with methods by other authors that published their results using the same database [2,6,64]. The first four rows of the table show the results of the manual annotations (observer 1, 2 and 3) and the golden standard. The next four rows show results of four classic approaches to vessel width measurements: -Gregson: a rectangle is fitted to a vessel intensity profile, and the width is set such that the area under the rectangle and profile [36] are equal. -Half Height Full Width (HHFW): the standard halfheight method, which uses thresholds set half-way between the maximum and minimum intensities at either side of an estimated center point [11]. The last row shows the results we achieved using our ETOS algorithm.
The column labeled with % shows the success rate, it indicates how many width measurements could successfully be validated (for more detail see [2]). The success percentage is smaller then 100% whenever mea-surements failed to converge, e.g. when the distance between the ground truth and measured edge pair was too large. The column labeled with Mean indicates the mean vessel width of all the measured vessel profiles. The column labeled with σ χ indicates the standard deviation of the error (Eq. (40)), a lower σ χ is favorable since it indicates that the error is consistent.
From Table 1 it can be observed that ESP, Graph, ARIA and our ETOS algorithm all outperform the classic width measurement techniques. Also compared to the state of the art methods our algorithm scores very well. The ETOS algorithm performs remarkably well on the CLRIS dataset, which contains a large number of vessels with the central light reflex. For these images, the standard deviation of the errors is even lower than those of the observers. For other datasets, our method's performance is comparable to the state of the art. Fig. 12 shows a scatter plot of the ground truth widths against the widths measured by our ETOS algorithm, together with a linear regression model that was fit through these points. The points are very much centered around the line y = x, indicating a strong positive correlation. This is confirmed by the slope of the linear regression model y = 0.85 + 0.88x, which is near to 1. The low number of outliers in the scatter plot confirms the low standard deviation in errors, as demonstrated by Table 1. The offset of 0.85, together with slope 0.88, indicate that the ETOS algorithm has the tendency to slightly overestimate for vessels of size up to 7 pixels and underestimates for larger vessel sizes.
We conclude that our ETOS algorithm, which is highly capable of tracking blood vessels through all sorts of complex situations, also provides reliable width measurements.

Vasculature tracking
In this section we describe additions to the ETOS algorithm, so as to be able to construct models of the complete retinal vasculature. Our vasculature tracking algorithm consists of: 1. Optic disk detection. 2. Seed point detection in the optic disk region. 3. Correct Initialization of the ETOS algorith by robust initial edge detection. 4. Automatic termination based on a set of stopping criteria. 5. Junction detection, classification and numbering. 6. Junction resolving.
Each of these items is described in Sections 4.1.1 to 4.1.6 and the complete algorithm is validated in Section 4.2.

Optic disk detection
Since the blood vessels of interest all enter the eye through the optic disk, vasculature tracking is initiated in this region. The detection of the optic disk occurs in two phases. In the first phase a rough estimation of the optic disk position is made based on a method using variance filtering, as proposed by Sinthanayothin [54]. This method is based on the significant variance in pixel intensities that occur in the optic disk region.
Next, the estimated optic disk position is refined by edge focusing [8] on the optic disk boundary, followed by a weighted Hough transform [38] for circles. Here we assume that the optic disk is approximately circular. To avoid disturbance of blood vessels during the detection of the optic disk boundary, the vessels are first removed by applying a closing operator on the red channel of the image (Fig. 14b). Edge focussing is performed on a star-shaped set of profiles (Fig. 14c) and the detected edge positions are used as input for the weighted Hough transform for circles. The weight of each edge position is determined by the scale up to which the edge can be traced in Gaussian scale space (before it reaches a so-called toppoint [29,40]).
The optic disk radius R OD , found by the Hough transform for circles, describes the size of the optic disk. For adult human eyes, the average radius of the optic disk is known to be R ref OD = 0.92mm [57]. Vessels in the optic disk region typically have a caliber of w ref av = 0.15mm. While the resolution of retinal images may vary from camera to camera, the physical dimensions of the human eye are rather constant between individuals. In order to make our algorithm generally applicable to retinal images of different resolution, we will normalize distances used in our routines with where, based on real physical values, w av describes typical retinal vessel calibers in pixels.

Vessel likelihood map and seed point detection
For the detection of the initial seed points, a vessel likelihood map V : R 2 → R of the optic disk region is constructed using invertible orientation scores. Before images are subjected to our algorithms, they have their DC-component removed by means of high-pass filtering. As a result vessel pixels have negative values and background pixel values are assumed to be zero. In effect the real part of the orientation score at a certain vessel position and orientation is negative, the orientation score is approximately zero at background area's. This is also demonstrated in Fig. 7b, were the vessel in the score can be seen as a dark (negative value) blob on the plane. Based on these observations, we define a vessel likelihood map V (x) as: Seed points are detected as local maxima on circular intensity profiles centered around the optic disk ( Fig. 15a) with radii R = {R OD , 1.5R OD }, where R OD is the detected optic disk radius. A seed point is discarded whenever its value in the vessel likelihood map V is smaller then the average value of all points on the circle. For each remaining seed point c 0 , the initial orientation θ 0 is detected as the orientation that provides the highest modulus of the orientation scores: θ 0 = argmax |U f (c 0 , θ)|. An additional filtering step, in which the seed points are classified as either true or false positive, is described in Section 4.1.3.

Initial edge detection
The ETOS algorithm is initialized with a starting vessel center point c 0 , orientation θ 0 and edges u 0 and v 0 . Starting with an already detected initial center point c 0 and orientation θ 0 , intensity profile I η can be obtained from the orientation scores using Eq. (29). Candidate edges are detected as local optima on the imaginary part of I η . Beside the main vessel edges that we are interested in, it is very likely that multiple other candidate edges are detected as well (as a result of noise or a central light reflex). Therefore, we use an edge focussing approach to detect the dominant edges. Each combination of neighboring left and right edges will form potential vessel patches (see Fig. 14a). Note that a blood vessel with a central light reflex consists of two neighboring vessel patches. We start initial edge detection by detecting the main vessel patch of interest by scoring each edge pair (u,v), based on the initial center point estimate c 0 and initial orientation θ 0 , as follows: with The function ν θ (u, v) provides the so called vessel value and basically is the average value of the modulus of the orientation score at orientation θ, calculated from the left to the right edge. This value is high for elongated/vessel structures, and low for background structures in the score. The exponential in Eq. (43) penalizes the distance of the edge pair to the initialized center point. The edge pair with the highest score will be our main pair of interest.
Neighboring pairs are considered only if the distance between the nearest edges of the main and neighboring pair is smaller than the width of the main patch, Fig. 14a shows typical results. All edges from the main patch and its neighbors are now traced (in scale) in the Gaussian scale space of profile I η . The edges are traced up to the scale of the first appearing toppoint. The strongest edges at this scale are chosen as the true vessel edge points u 0 and v 0 (Fig. 14c).
An initialized seed point whose vessel value is lower than a threshold T ν is regarded as a false positive (Fig. 15b). The threshold T ν is defined as: where using all N sp initialized seed points, each with initial edge points (u i 0 , v i 0 ) and orientation (θ i 0 ), an average vessel value ν av is calculated. All true positive seed points are submitted to the ETOS algorithm to start expanding a model of the retinal vasculature.

Stopping criteria
For ETOS to stop tracking a single vessel, three stoppingcriteria are defined: 1. Tracking stops whenever a vessel being tracked leaves a prescribed region of interest. For this purpose a mask image is generated that covers the (typically circular) field of view in the fundus image. 2. Tracking stops whenever a blood vessel is already tracked. Each tracked segment is used to construct a pixel map, in which each pixel within the vessel edge contours is set to 1 and outside to 0. Whenever a point within the vessel being tracked lies within the pixel map for 4 w av /λ iterations in a row, the ETOS algorithm terminates. Recall w av defined in Eq. (41) and λ being the step size (Fig. 8). 3. Tracking stops whenever the vessel value ν θ (Eq. (44)) drops below threshold value T ν . Here we assumed that T ν , which is based on the average vessel value ν av calculated in the optic disk region, is a good indicator of the vessel values of the vessels elsewhere in the retina.

Junction detection, classification and numbering
In order to model the complete vasculature, starting by expanding a model from a set of initial seed points, the vasculature tracking algorithm should also be able to automatically detect junctions. Junction points are points where blood vessels bifurcate/branch or where two blood vessels cross. Either way, at a junction point multiple orientations may be observed. At a point on a straight vessel the modulus of the directional orientation column ideally contains two local maxima; one in the positive direction, say θ + , and one (with a π phase difference) in the negative direction, say θ − . A candidate vessel junction point at the left edge is detected whenever there exists another local maximum in between θ + and θ − . Similarly, a candidate junction point is detected at the right vessel edge whenever another local maximum is present between θ − and θ + . During vessel tracking, the orientation columns at the left and right edge are scanned for the presence of a junction point, and the location and orientation are stored (top image Fig. 17). The detected junction points are clustered on position by grouping all points whose distance to one another is smaller than w av . Within each cluster, the candidate junction points are clustered on orientation to prevent merging two proximate junction points. Clustering on orientation is done according to the number of local maxima in the histogram of orientations. Finally, clusters are merged to a single junction point by averaging the positions and by taking the most common orientation within the cluster. The found center points and orientations are then subjected to the edge initialization method described in  To classify between junctions and bifurcations we check for alignment of local orientations in SE(2) via the sub-Riemanian metric [34,9], see Eq. (51) in Appendix A.
Each detected seed point is assigned two ID numbers, a number that is unique for each vessel segment and the ID number of the vessel from which it originates. This way the relation between each vessel segment remains known, and the vessel segments can be organized in a hierarchical fashion. Vasculature tracking terminates whenever all detected bifurcations and crossings are evaluated.

Junction resolving
As mentioned in Section 4.1.4 vessel tracking is terminated whenever the algorithm is tracking a vessel that is already tracked. This criterion can be met at several situations, where for each situation appropriate actions need to be taken in order to maintain correct topological models of the vasculature. A detected point that suggests inappropriate modeling of the vasculature will be called an unresolved junction point. The detection of an unresolved junction point, together with the corresponding actions that are necessary to solve it will be called junction resolving. The appropriate actions necessary for junction resolving are based on the position of the junction point on the already established track, whether or not the two overlapping segments have the same with and in the case the The junction resolving steps are described using the following labeling: S old is the tracked segment of the established track before the junction point (Source), T old is the tracked segment of the same track after the junction point (Target), S new is the segment of the new track before the junction point, and T new is the segment of the new track after the junction point. The appropriate actions necessary for junction resolving are based on 1. the position of the junction point on the already established track (S new -T new ), 2. whether or not the two overlapping segments (T new and T old ) have the same width, 3. and the similarity between the source and target tracks.
The complete junction resolving scheme and all detail can be found in [7].

Validation
In Section 3.2.2 we demonstrated the reliability of the width measurements provided by the ETOS algorithm.
In the following section we validate the topological correctness of the complete vasculature models that are generated by our algorithm. The correctness of the models is validated by analyzing the junction points. The results discussed in this section are generated with the same parameters for the ETOS algorithm that are described in Section 3.2.1. A typical model generated by the vasculature tracking algorithm with these parameters is shown in Fig. 18.

Validation of topological correctness
For each vessel it is known from which parent vessel it originates and bifurcations can thus be directly extracted from the model. Crossings can easily be extracted by detecting overlapping vessel segments. Fig. 18b provides an overview of detected junction points for one   image of the HRFI-database, Fig. 19 shows detailed views of a bifurcations and a crossing. We evaluated the junction points of the first three images from each of the three HRFI-datasets (healthy, diabetes and glaucoma). The 9 generated vasculature models provided 495 junction points, of which 381 were bifurcation points and 114 were crossing points. The following types of errors for bifurcations were identified: the bifurcation was actually part of a crossing (E1), the vessel originating from the bifurcation did not represent a blood vessel according to the ground truth pixel map provided by the HRFI-database (E2). E2 errors indicate the presence of incorrect single vessel models, non-vessel elongated structures such as the optic disk border or pathological features such as aneurysms. Crossings are extracted by searching the vasculature model for overlapping vessel segments. A false positive crossing can thus only occur if false positive vessel segments exists within the models. The results are summarized in Table 2.
In total 290 out of 381 bifurcations and 109 out of 114 crossings were correctly detected, corresponding to precision rates of 76.54% and 96.03% respectively. Most of the incorrectly identified bifurcations were correct in the sense that they represent a true blood vessel, however they were actually part of a crossing. Only 5.25% of the bifurcations were incorrect in the sense that they did not represent a blood vessel. While the first kind of false positive bifurcations only affect the topological correctness of the model, the latter also pollutes the model with false positive vessel segments. The low percentage of E2 errors indicates that the generated vasculature models are very clean in the sense that almost all vessel segments actually represent true blood vessels.
Note that E1 errors are introduced by incorrect junction classification and that each miss-classified crossing introduces two E1 errors. Additional post-processing steps that try to solve for correct topology [3] may be used to filter out these errors.

Conclusion
In this paper we demonstrated that by representing image information in an invertible orientation score, one can exploit the disentanglement of crossing structures to track blood vessels through crossing points. We introduced a new algorithm that tracks vessel edges through the orientation score of an image (ETOS). The ETOS algorithm can generally be used with both invertible and non-invertible orientation scores, which were in this paper constructed with cake wavelets and Gabor wavelets respectively. We demonstrated that best results were obtained using invertible orientation scores. We also introduced a fast alternative method based on vessel centerline tracking through a multi-scale set of noninvertible orientation scores (CTOS). While the CTOS algorithm is very fast, the multi-scale approach makes the algorithm less stable at critical vessel points (crossings and parallel vessels) compared to ETOS.
The ETOS algorithm was used as a basis for our vasculature tracking algorithm, which we used to construct detailed hierarchical models of the retinal vasculature. Within this paper we validated the reliability of the width measurements provided by the models using ground truth data, and showed that our method performs excellent in comparison to other state of the art algorithms. Validation of the topology of the models showed that our method constructs clean topological models of the vasculare tree, i.e. they contain very few false positive vessels.
Tracking within orientation scores relies on a novel and basic geometrical principle (V-plane optimization) within a sub-Riemanian manifold within SE (2). In Appendix A we relate probabilistic approaches to optimal curves in SE (2) to this geometric principle. Here we include supporting examples of analytic and numeric computations on completion fields on SE (2). Such completion fields are obtained from collision probabilities on SE(2), while the orientation score provide a whole distribution of these particles. The next step in future work is to apply tracking by geometric control [19] in enhanced orientation scores, obtained by applying stochastic processes for [24,18] contour completion and enhancement directly on the score.
A A mathematical underpinning of optimization in the e η − e θ tangent planes V Before we provide the key-geometrical principle behind our ETOS algorithm, we make some comments on the moving frame of reference that lives in the domain of an orientation score.
The domain of an orientation requires a curved (non-Euclidean) geometry as can be seen in Fig. 1. To this end we identify the coupled space of positions and orientations R 2 S 1 with SE(2). This needs a bit of explanation. To see whether a local line element (x , θ ) ∈ R 2 × S 1 with position x ∈ R 2 and orientation θ is aligned with its neighbors one applies a rigid body motion g = (x, R θ ) ∈ SE(2), with counter-clockwise rotation and translation via Eq. (10). This puts a 1-to-1 correspondence between rigid body motions (x, y, R θ ) and elements from R 2 × S 1 as every local line element (x, y, θ) can be obtained by a local line element (0, 0, 0) positioned at the origin 0 = (0, 0) pointing in say x-direction: Via this identification we see that the curved domain of an orientation score is due to the non-commutative group structure of the rigid body motion (also known as roto-translation group or Euclidean motion group) SE(2) in the plane. Indeed, a concatenation of two rigid body motions is again a rigid body motion and this induces the group-product given by Eq. (11) It is a semi-direct product of rotations SO(2) and translations R 2 , as the rotation part affects the position part. As a result this group product is non-commutative, i.e. gg can be very different from g g. Therefore it is conceptually wrong to consider the domain of an orientation scores as the flat Cartesian space R 2 × S 1 with Euclidean metric and it is wrong to just take derivatives of orientation scores only in ∂ x direction or only in ∂ y direction. For formal underpinning of this statement see [17,Thm.21], for a short illustration see [31, Fig.2.6]. Instead one must differentiate along the moving frame of ref- This a priori moving frame of reference coincides with the socalled Lie algebra of left-invariant vector fields on the group SE(2). The curved geometry in the orientation score requires us to apply parallel transport along this moving frame of reference. For instance in Fig. 1(d) we see examples of autoparallel curves whose tangent vector is always constant to the moving frame of reference. The non-commutative structure of SE(2) (due to the fact that rotations and translations do not commute) is reflected by the commutators in the Lie-algebra As a result first moving in ξ direction and then moving in θ direction brings us to a different point in SE(2) ≡ R 2 ×S 1 than first moving in θ direction and then moving in ξ direction (as can be deduced in Fig. 1(b)).
To this end we recall Fig. 7 where one can observe that each tangent vector to the lifted curve s → (x(s), y(s), θ(s)) with θ(s) = arg(ẋ(s) + iẏ(s)) (e.g. the curve following the right edge of a blood vessel) in SE(2) is pointing orthogonal to a plane V (plotted in yellow). Within V we see that the maximum values of the absolute value of the imaginary part of the score is located at the origin of the yellow plane in the tangent space. (2) is called the lifted curve of a smooth planar curve iff for all s ∈ [0, ] we have θ(s) = arg(ẋ(s) + iẏ(s)).
Definition 2 If a curve γ is equal to the lift of its spatial projection (i.e. if its satisfies Eq. (46)) it is called horizontal.
Intuitively, such a sub-Riemannian manifold M equals the group SE(2) where one restricts oneself to horizontal curves with a constant relative penalty for bending and stretching determined by β > 0 which has physical dimension one over length.
For special cases of C we can show that our geometrical principle indeed produces optimal curves in SE (2), as we will show in the subsequent section.

A.2 Application of the geometrical principle to completion fields
In general the real part, imaginary part, or absolute value of an orientation score is a complicated function on SE(2). Therefore we will consider the case where C : SE(2) → R + is a so-called "completion field" [60,67,4,59,5]. This corresponds to collision probability densities of a source particle g 1 ∈ SE(2) and a sink particle g 2 ∈ SE(2).
There exist remarkable relations [21,18] between optimal curves (i.e. curves minimizing an optimal control problem on SE(2)) and solutions of Eq. (47) for special cases where C denotes a so-called completion distribution (or "completion field"). Given two sources at the origin g 1 = (0, 0, 0) and 6 In such a way that the perturbed curve is again horizontal.  (2)) and its planar θintegrated version. Bottom: A contour enhancement process the corresponding Green's function (in SE (2)) and its planar θ-integrated version. In the completion process one has randomness in θ with variance 2D 11 > 0 and non-random advection in ξ. In the enhancement process one has randomness both in θ-direction (with variance 2D 11 > 0) and in ξ-direction (with variance 2D 22 > 0). at g 2 = (x 1 , y 1 , θ 1 ), such completion fields are defined as products of resolvent Green's functions of stochastic processes for contour completion [45] and contour enhancement [24,34] in SE(2): where R λ (g) denotes the probability density of finding a random walker g in the underlying stochastic process [24,45] given that it started at g = (x, y, θ) = (0, 0, 0) regardless its memoryless traveling time T which is negatively exponentially distributed with expectation E(T ) = λ −1 . Furthermore, g → R λ, * (g) denotes the adjoint resolvent kernel (i.e. the resolvent that arises by taking the adjoint of the generator [18, ch:4.4]. For a concise overview on (Green's functions of ) contour completion and enhancement see [51]. Exact formulas for the resolvent Green's function R λ for contour enhancement can be found in [24], whereas exact formulae for resolvent kernels R λ for contour completion (direction process) can be found in [18]. In both cases there are representations involving 4 Mathieu functions. For a visual impression of exact Green's functions see Fig. 20 (where in dashed lines we have depicted level sets of the corresponding Heisenberg approximations that we will discuss and employ in the next subsection).
These completion fields relate to the well-known Brownian bridges (in probability theory) where the traveling time is integrated out. This relation is relevant, since it is known that such Brownian bridges concentrate on geodesics [62].
On the other hand, in his paper [45] Mumford showed that the modes of the direction process (also known as the contour completion process) coincide with elastica curves which are the solutions to the following optimal control problem inf x ∈ C 1 [0, ], > 0 (x(0),ẋ(0)) = g 0 , (x(l),ẋ( )) = g 1 0 κ 2 (s) + β 2 ds (52) Similar to the Onsager-Machlup approach to optimal paths [56] he obtains these modes by looking at the most probable/likely realization of discretized versions of the direction process. As this cannot be (efficiently) realized in practice, one needs a more tangible description of the mode. To this end, we will call solution curves of (47) the "modes". In case of the direction process and in case one uses the completion distribution (Eq. (50)) for the function C(g), the solution curves of (47) indeed seem to numerically coincide with the elastica. This experimentally underpins our conjecture in [18], where we have shown such a result does hold exactly for the corresponding Heisenberg approximations as we will explain next. 7 This means g −1 2 g 1 is within the range of the exponential map of the sub-Riemannian control problem, i.e. g 1 and g 2 can be connected via a stationary curve without cusps [10,9,19].

A.3 In the Heisenberg approximation of completion fields our approach produces B-splines
The Heisenberg group approximation (obtained by contraction [24]) of the Green's functions and induced completion field arises by replacing the moving frame of left-invariant vector fields {cos θ∂x + sin θ∂ y , − sin θ∂ x + cos θ∂ y , ∂ θ } on SE(2) by the moving frame of reference of left-invariant vector fields on the Heisenberg group {∂ x + θ∂ y , ∂ y , ∂ θ } and by replacing spatial arc-lengt parametrization via s by spatial coordinate x. Intuitively, such replacement boils down to replacing the space of positions and orientations by the space of positions and velocities.

A.4 Concluding remarks
By the results of the previous section the conjecture rises whether elastica curves coincide with Eq. (50) as such relation holds for their counterparts in the Heisenberg group. Numeric computations seem to provide a confirmation of this conjecture, see Fig. 22.
On the one hand, we expect with respect to the contour enhancement process that our approach produces the sub-Riemannian geodesics (based on the results in [62,23,21]), but this is a point for future investigation. On the other hand, the conjecture together with the result in [45] (and the result that B-splines solve Eq. (53)) would underpin our alternative applicable definition of modes as solution curves of Eq. (47).
In fact this means that optimization in each (η, θ)-plane V, as is done in our ETOS algorithm, produces the most probable curves (in the sense of [45,56]) in direction processes.
B The closure of the distributional orientation score transform is an L 2 -isometry on the whole L 2 (R 2 ) space.
Under a certain condition on ψ, we show that operator W is an isometry from H I (R 2 ) (with L 2 -norm) into L 2 (SE(2)). Therefore this operator is closable and its closure is an isometry. This bring us to the main result.
Concluding Remark By the result of the previous theorem, ψ ∈ H −I (R 2 ) with M D −1 β ψ = β −2 can be called proper distributional wavelets. When insisting on an L 2 -isometric mapping between image and score one has to fall back on these kind of wavelets. In case of cake wavelets (proper wavelets of class I [17, Ch. 4.6.1]), when → ∞ such wavelets typically become oriented δ-distributions.
In case of proper wavelets of class II [17, Ch. 4.6.2] (including the kernel proposed by Kalitzin [41]) such wavelets concentrate around and explode along the x-axis when N → ∞, [17,Fig. 4.11 and Ch. 7.3]. In both cases the limits do not exists in L 2 -sense, but they do exist both pointwise and in H −I -sense.
We conclude from the results in this section that the orientation score framework does not insist on images to be bandlimited, and remains valid regardless the sampling size/rate.