extfieldsofstudy:
Mathematics
Computer Science

text:
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.


Using context and phonetic features in models of etymological sound change

This paper presents a novel method for aligning etymological data, which models context-sensitive rules governing sound change, and utilizes phonetic features of the sounds. The goal is, for a given corpus of cognate sets, to find the best alignment at the sound level. We introduce an imputation procedure to compare the goodness of the resulting models, as well as the goodness of the data sets. We present evaluations to demonstrate that the new model yields improvements in performance, compared to previously reported models.

Introduction
This paper introduces a context-sensitive model for alignment and analysis of etymological data. Given a raw collection of etymological data (the corpus)-we first aim to find the "best" alignment at the sound or symbol level. We take the corpus (or possibly several different corpora) for a language family as given; different data sets are typically conflicting, which creates the need to determine which is more correct. Etymological data sets are found in digital etymological databases, such as ones we use for the Uralic language family. A database is typically organized into cognate sets; all elements within a cognate set are posited (by the database creators) to be derived from a common origin, which is a word-form in the ancestral proto-language.
Etymology encompasses several problems, including: discovery of sets of cognatesgenetically related words; determination of genetic relations among groups of languages, based on linguistic data; discovering regular sound correspondences across languages in a given lan-guage family; and reconstruction of forms in the proto-languages.
Computational methods can provide valuable tools for the etymological community. The methods can be judged by how well they model certain aspects of etymology, and by whether the automatic analysis produces results that match theories established by manual analysis.
In this work, we allow all the data-and only the data-to determine what rules underly it, rather than relying on external (and possibly biased) rules that try to explain the data. This approach will provide a means of measuring the quality of the etymological data sets in terms of their internal consistency-a dataset that is more consistent should receive a higher score. We seek methods that analyze the data automatically, in an unsupervised fashion, to determine whether a complete description of the correspondences can be discovered automatically, directly from raw etymological data-cognate sets within the language family. Another way to state the question is: what alignment rules are "inherently encoded" in the given corpus itself.
At present, our aim is to analyze given etymological datasets, rather than to construct new ones from scratch. Because our main goal is to develop methods that are as objective as possible, the models make no a priori assumptions or "universal" principles-e.g., no preference to align vowel with vowels, or a symbol with itself. The models are not aware of the identity of a symbol across languages, and do not try to preserve identity, of symbols, or even of features-rather they try to find maximally regular correspondences.
In Section 2 we describe the data used in our experiments, and review approaches to etymological alignment over the last decade. We formalize the problem of alignment in Section 3, give the Uralic tree technical details of our models in Section 4. We present results and discussion in Sections 5 and 6.

Data and Related Work
We use two large Uralic etymological resources. The StarLing database of Uralic, (Starostin, 2005), based on (Rédei, 1988(Rédei, 1991, contains over 2500 cognate sets. Suomen Sanojen Alkuperä (SSA), "The Origin of Finnish Words", a Finnish etymological dictionary, (Itkonen and Kulonen, 2000), has over 5000 cognate sets, (about half of which are only in languages from the Balto-Finnic branch, closest to Finnish). Most importantly, for our models, SSA gives "dictionary" word-forms, which may contain extraneous morphological material, whereas StarLing data is mostly stemmed.
One traditional arrangement of the Uralic languages 1 is shown in Figure 1. We model etymological processes using these Uralic datasets.
The methods in (Kondrak, 2002) learn regular one-to-one sound correspondences between pairs of related languages in the data. The methods in (Kondrak, 2003;Wettig et al., 2011) find more complex (one-to-many) correspondences. These models operate on one language pair at a time; also, they do not model the context of the sound changes, while most etymological changes are conditioned on context. The MCMC-based model proposed in (Bouchard-Côté et al., 2007) explicitly aims to model the context of changes, and op-1 Adapted from Encyclopedia Britannica and (Anttila, 1989) erates on more than a pair of languages. 2 We should note that our models at present operate at the phonetic level only, they leave semantic judgements of the database creators unquestioned. While other work, e.g. (Kondrak, 2004), has attempted to approach semantics by computational means as well, our model uses the given cognate set as the fundamental unit. In our work, we do not attempt the problem of discovering cognates, addressed, e.g., in, (Bouchard-Côté et al., 2007;Kondrak, 2004;Kessler, 2001). We begin instead with a set of etymological data (or more than one set) for a language family as given. We focus on the principle of recurrent sound correspondence, as in much of the literature, including (Kondrak, 2002;Kondrak, 2003), and others.
As we develop our alignment models at the sound or symbol level, in the process of evaluation of these models, we also arrive at modeling the relationships among groups of languages within the family. Construction of phylogenies is studied extensively, e.g., by (Nakhleh et al., 2005;Ringe et al., 2002;Barbançon et al., 2009). This work differs from ours in that it operates on manually pre-selected sets of characters, which capture divergent features of languages within the family, whereas we operate on the raw, complete data.
There is extensive work on alignment in the machine-translation (MT) community, and it has been observed that methods from MT alignment may be projected onto alignment in etymology. The intuition is that translation sentences in MT correspond to cognate words in etymology, while words in MT correspond to sounds in etymology. The notion of regularity of sound change in etymology, which is what our models try to capture, is loosely similar to contextually conditioned correspondence of translation words across languages. For example, (Kondrak, 2002) employs MT alignment from (Melamed, 1997;Melamed, 2000); one might employ the IBM models for MT alignment, (Brown et al., 1993), or the HMM model, (Vogel et al., 1996). Of the MT-related models, (Bodrumlu et al., 2009) is similar to ours in that it is based on MDL (the Minimum Description Length Principle, introduced below).

Aligning Pairs of Words
We begin with pairwise alignment: aligning pairs of words, from two related languages in our corpus of cognates. For each word pair, the task of alignment means finding exactly which symbols correspond. Some symbols may align with "themselves" (i.e., with similar or identical sounds), while others may have undergone changes during the time when the two related languages have been evolving separately. The simplest form of such alignment at the symbol level is a pair (σ : τ ) ∈ Σ × T , a single symbol σ from the source alphabet Σ with a symbol τ from the target alphabet T . We denote the sizes of the alphabets by |Σ| and |T |.
To model insertions and deletions, we augment both alphabets with a special empty symboldenoted by a dot-and write the augmented alphabets as Σ . and T . . We can then align word pairs such as vuosi-al (meaning "year" in Finnish and Xanty) , for example as any of:

Context Model with Phonetic Features
The context-aware alignment method we present here is built upon baseline models published previously, (Wettig et al., 2011), where we presented several models that do not use phonetic features or context. Similarly to the earlier ones, the current method is based on the Minimum Description Length (MDL) Principle, (Grünwald, 2007).
We begin with a raw set of (observed) datathe not-yet-aligned word pairs. We would like to find an alignment for the data-which we will call the complete data-complete with alignments, that make the most sense globally, in terms of embodying regular correspondences. We are after the regularity, and the more regularity we can find, the "better" our alignment will be (its goodness will be defined formally later). MDL tells us that the more regularity we can find in the data, the fewer bits we will need to encode it (or compress it). More regularity means lower entropy in the distribution that describes the data, and lower entropy allows us to construct a more economical code. That is, if we have no knowledge about any regularly of correspondence between symbols, the joint distribution over all possible pairs of symbols will be very flat (high entropy). If we know that certain symbol pairs align frequently, the joint distribution will have spikes, and lower entropy. In (Wettig et al., 2011) we showed how starting with a random alignment a good joint distribution can be learned using MDL. However the "rules" those baseline models were able to learn were very rudimentary, since they could not use any information in the context, and we know that many regular correspondences are conditioned by context.
We now introduce models that leverage information from the context to try to reduce the uncertainty in the distributions further, lowering the coding cost. To do that, we will code sounds in terms of their phonetic features: rather than coding the symbols (sounds) as atomic, we code them as vectors of phonetic features. Rather than aligning symbol pairs, we align the corresponding features of the symbols. While coding each feature, the model can make use of features of other sounds in its context (environment), through a special decision tree built for that feature.

Features
We will code each symbol, to be aligned in the complete data, as a feature vector. First we code the Type feature, with values: K (consonant), V (vowel), dot, and word boundary, which we denote as #. Consonants and vowels have their own sets of features, with 2-8 values per feature:

Vowel articulation V Vertical
high-low H Horizontal front-back R Rounding -, + L Length 1-5

Contexts
While coding any symbol, the model will be allowed to query a fixed, finite set of candidate contexts. A context is a triplet (L, P, F ), where L is the level-either source or target,-and P is one of the positions that the model may queryrelative to the position currently being coded; for example, we may allow positions as in Fig. 2. F is one of the possible features found at that position. Therefore, we will have about 2 levels * 8 positions * 2-6 features ≈ 80 candidate contexts that can be queried by the model, as explained below.
I itself, -P previous position -S previous non-dot symbol -K previous consonant -V previous vowel +S previous or self non-dot symbol +K previous or self consonant +V previous or self vowel

The Two-Part Code
We code the complete (i.e., aligned) data using a two-part code, following the MDL Principle. We first code which particular model instance we select from our class of models, and then code the data, given the defined model. Our model class is defined as: a set of decision trees (forest), with one tree to predict each feature on each level. The model instance will define the particular structures for each of the trees.
The forest consists of 18 decision trees, one for each feature on the source and the target level: the type feature, 4 vowel and 4 consonant features, times 2 levels. Each node in such tree will either be a leaf, or will be split by querying one of the candidate contexts defined above. The cost of coding the structure of the tree is one bit for every node-to encode whether this node was split (is an internal node) or is a leaf-plus ≈ log 80 times the number of internal nodes-to encode which particular context was chosen to split that node. We will explain how the best context to split on is chosen in Sec. 4.6.
Each feature and level define a tree, e.g., the "voiced" (X) feature of the source symbols corresponds to the source-X tree. A node N in this tree holds a distribution over the values of X of only those symbol instances in the complete data that have reached in N by following the context queries, starting from the root. The tree structure tells us precisely which path to followcompletely determined by the context. For example, when coding a symbol α based on another symbol found in the context of α-at some level (say, target), some position (say, -K), and one of its features (say, M)-the next edge down the tree is determined by that feature's value; and so on, down to a leaf. For an example of an actual decision tree learned by the model, see Fig. 5.
To compute the code length of the complete data, we only need to take into account the distributions at the leaves. We could choose from a variety of coding methods; the crucial point is that the chosen code will assign a particular numberthe cost-to every possible alignment of the data. This code-length, or cost, will then serve as the objective function-i.e., it will be the value that the algorithm will try to optimize. Each reduction in cost will correspond directly to reduction in the entropy of the probability distribution of the symbols, which in turn corresponds to more certainty (i.e., regularity) in the correspondences among the symbols, and to improvement in the alignment. This is the link to our goal, and the reason for introducing code lengths-it gives us a single number that describes the quality of an alignment.
We use Normalized Maximum Likelihood (NML), (Rissanen, 1996) as our coding scheme. We choose NML because it has certain optimality properties. Using NML, we code the distribution at each leaf node separately, and summing the costs of all leaves gives the total cost of the aligned data-the value of our objective function.
Suppose n instances end up in a leaf node N , of the λ-level tree, for feature F having k values (e.g., consonants satisfying N 's context constraints in the source-X tree, with k = 2 values: − and +), and the values are distributed so that n i instances have value i (with i ∈ {1, . . . , k}). Then this requires an NML code-length of Here i n i n n i is the maximum likelihood of the multinomial data at node N , and the term C(n, k) = is a normalizing constant to make P N M L a probability distribution.
In the MDL literature, e.g., (Grünwald, 2007), the term − log C(n, k) is called the stochastic complexity or the (minimax) regret of the model, (in this case, the multinomial model). The NML distribution provides the unique solution to the minimax problem posed in (Shtarkov, 1987), whereΘ(x n ) = arg max Θ P(x n ) are the maximum likelihood parameters for the data x n . Thus, P N M L minimizes the worst-case regret, i.e., the number of excess bits in the code as compared to the best model in the model class, with hind-sight. For details on the computation of this code length see (Kontkanen and Myllymäki, 2007). Learning the model from the observed data now means aligning the word pairs and building the decision trees in such a way as to minimize the two-part code length: the sum of the model's code length-to encode the structure of the trees,and the data's code length-to encode the aligned word pairs, using these trees.

Summary of the Algorithm
The full learning algorithm runs as follows: We start with an initial random alignment for each pair of words in the corpus, i.e., for each word pair choose some random path through the matrix depicted in Figure 3.
From then on we alternate between two steps: A. re-build the decision trees for all features on source and target levels, and B. re-align all word pairs in the corpus. Both of these operations monotonically decrease the two-part cost function and thus compress the data.
We continue until we reach convergence.

Re-alignment Procedure
To align source word σ consisting of symbols σ = [σ 1 ...σ n ], σ ∈ Σ * with target word τ = [τ 1 ...τ m ] we use dynamic programming. The tree structures are considered fixed, as are the alignments of all word pairs, except the one currently being aligned-which is subtracted from the counts stored at the leaf nodes. We now fill the matrix V , left-to-right, top-tobottom. Every possible alignment of σ and τ cor- Figure 3: Dynamic programming matrix V, to search for the most probable alignment responds to exactly one path through this matrix: starting with cost equal to 0 in the top-left cell, moving only downward or rightward, and terminating in the bottom-right cell. In this Viterbi-like matrix, every cell corresponds to a partially completed alignment: reaching cell (i, j) means having read off i symbols of the source word and j symbols of the target. Each cell V (i, j)-marked X in the Figure-stores the cost of the most probable path so far: the most probable way to have scanned σ through symbol σ i and τ through τ j : Each term V (·, ·) has been computed earlier by the dynamic programming; the term L(·)-the cost of aligning the two symbols, inserting or deleting-is determined by the change in data code length it induces to add this event to the corresponding leaf in all the feature trees it concerns.
In particular, the cost of the most probable complete alignment of the two words will be stored in the bottom-right cell, V (n, m), marked .

Building Decision Trees
Given a complete alignment of the data, we need to build a decision tree, for each feature on both levels, yielding the lowest two-part cost. The term "decision tree" is meant in a probabilistic sense here: instead of a single value, at each node we store a distribution of the corresponding feature values, over all instances that reach this node. The distribution at a leaf is then used to code an instance when it reaches the leaf in question. We code the features in some fixed, pre-set order, and source level before target level.
We now describe in detail the process of building the tree for feature X, for the source level, (we will need do the same for all other features, on both levels, as well). We build this tree as follows. First, we collect all instances of consonants on the source level, and gather the the counts for feature X; and build an initial count vector; suppose it is: value of X: + -1001 1002 This vector is stored at the root of the tree; the cost of this node is computed using NML, eq. 1.
Next, we try to split this node, by finding such a context that if we query the values of the feature in that context, it will help us reduce the entropy in this count vector. We check in turn all possible candidate contexts, (L, P, F ), and choose the best one. Each candidate refers to some symbol found on the source (σ) or the target (τ ) level, at some relative position P , and to one of that symbol's features F . We will condition the split on the possible values of F . For each candidate, we try to split on its feature's values, and collect the resulting alignment counts.
Suppose one such candidate is (σ, -V, H), i.e., (source-level, previous vowel, Horizontal feature), and suppose that the H-feature has two values: front/back. The vector at the root node (recall, this tree is for the X-feature) would then split into two vectors, e.g.: value of X: This would likely be a very good split, since it reduces the entropy of the distribution in each row almost to zero. The criterion that guides the choice of the best candidate to use for splitting a node is the sum of the code lengths of the resulting split vectors, and the code length is proportional to the entropy.
We go through all candidates exhaustively, and greedily choose the one that yields the greatest reduction in entropy, and drop in cost. We proceed recursively down the tree, trying to split nodes, and stop when the total tree cost stops decreasing.
This completes the tree for feature X on level σ. We build trees for all features and levels similarly, from the current alignment of the complete data.
We augment the set of possible values at every node with two additional special branches: =, meaning the symbol at the queried position is of the wrong type and does not have the queried feature, and #, meaning the query ran past the beginning of the word.

Evaluation and Results
One way to evaluate the presented models would require a gold-standard aligned corpus; the models produce alignments which could be compared to the gold-standard alignments, and we could measure performance quantitatively, e.g., in terms of accuracy. However, building a gold-standard aligned corpus for the Uralic data proved to be extremely difficult. In fact, it quickly becomes clear that this problem is at least as difficult as building a full reconstruction for all internal nodes in the family tree (and probably harder), since it requires full knowledge of all sound correspondences within the family. It is also compounded by the problem that the word-forms in the corpus may contain morphological material that is etymologically unrelated: some databases give "dictionary" forms, which contain extraneous affixes, and thereby obscure which parts of a given word form stand in etymological relationship with other members in the cognates set, and which do not. We therefore introduce other methods to evaluate the models.
Compression: In figure 4, we compare the context model, and use as baselines the standard data compressors, Gzip and Bzip, as well as the more basic models presented in (Wettig et al., 2011), (labeled "1x1 and "2x2"). We test the compression of up to 3200 Finnish-Estonian word pairs, from SSA. Gzip and Bzip compress data by finding regularities in it (i.e., frequent substrings). The comparison with Gzip is a "sanity check": we would like to confirm whether our models find more regularity in the data than would an off-the-shelf data compressor, that has no knowledge that the words in the data are etymologically related. Of course, our models know that they should align pairs of consecutive lines. This test shows that learning about the "vertical" correspondences achieves much better compression rates-allows the models to extract greater regularity from the data. Rules of correspondence: One our main goals is to model rules of correspondence among languages. We can evaluate the models based on how good they are at discovering rules. (Wettig et al., 2011) showed that aligning multiple symbols captures some of the context and thereby finds more complex rules than their 1-1 alignment model. However, certain alignments, such as t∼t/d, p∼p/b, and k∼k/g between Finnish and Estonian, cannot be explained by the multiple-symbol model. This is due to the rule of voicing of word-medial plosives in Estonian. This rule could be expressed in terms of Two-level Morphology, (Koskenniemi, 1983) as: a voiceless plosive in Finnish, may correspond to voiced in Estonian, if not word-initial. 3 The context model finds this rule, shown in Fig. 5. This tree codes the Target-level (i.e., Estonian) Voiced consonant feature. In each node, the counts of corresponding feature values are shown in brackets. In the root node-prior to knowing anything about the environment-there is almost complete uncertainty (i.e., high entropy) about the value of Voiced feature of an Estonian consonant: 821 voiceless to 801 voiced in our data. Redder nodes indicate higher entropy, bluer nodes-lower entropy. The query in the root node tells us to check the context Finnish Itself Voiced for the most informative clue about whether the current Estonian consonant is voiced or not. Tracing the options down left to right from the root, we obtain the rules. The leftmost branch says, if the Finnish is voiced (⊕), then the Estonian is almost certainly voiced as well-615 voiced to 2 voiceless in this case. If the Finnish is voiceless (Finnish Itself Voiced = ), it says voicing may occur, but only in the red nodes-i.e., only if preceded by a voiced consonant on Estonian level (the branch marked by ⊕, 56 cases), or-if previous position is not a consonant (the = branch indicates that the candidate's query does not apply: i.e., the sound found in that position is not a consonant)it can be voiced only if the corresponding Finnish is a plosive (P, 78 cases). The blue nodes in this branch say that otherwise, the Estonian consonant almost certainly remains voiceless.
The context models discover numerous complex rules for different language pairs. For example, they learn a rule that initial Finnish k "changes" (corresponds) to h in Hungarian, if it is followed by a back vowel; the correspondence between Komi trills and Udmurt sibilants; etc.
Imputation: We introduce a novel test of the quality of the models, by using them to impute unseen data, as follows. For a given model, and a language pair (L 1 , L 2 )-e.g., (Finnish, Estonian)-hold out one word pair, and train the model on the remaining data. Then show the model the hidden Finnish word and let it guess the corresponding Estonian. Imputation can be done for all models with a simple dynamic programming algorithm, similar to the Viterbi-like search used during training. Formally, given the hidden Finnish string, the imputation procedure selects from all possible Estonian strings the most probable Estonian string, given the model. We then compute an edit distance between the imputed sting and the true withheld Estonian word (e.g., using the Levenshtein distance). We repeat this procedure for all word pairs in the (L 1 , L 2 ) data set, sum the edit distances and normalize by the total size of the (true) L 2 data-this yields the Normalized Edit Distance N ED(L 2 |L 1 , M ) between L 1 and L 2 , under model M .
Imputation is a more intuitive measure of the model's quality than code length, with a clear practical interpretation. NED is also the ultimate test of the model's quality.
If model M imputes better than M -i.e., N ED(L 2 |L 1 , M ) < N ED(L 2 |L 1 , M )-then it is difficult to argue that M could be in any sense "worse" than Mit has learned more about the regularities between L 1 and L 2 , and it knows more about L 2 given L 1 . The context model, which has much lower cost than the baseline, almost always has lower NED. This also yields an important insight: it is an encouraging indication that optimizing the code length is a good approach-the algorithm does not optimize NED directly, and yet the cost correlates strongly with NED, which is a simple and intuitive measure of the model's quality.

Discussion
We have presented a novel feature-based contextaware MDL model, and a comparison of its performance against prior models for the task of alignment of etymological data. We have evaluated the models by examining the the rules of correspondence that they discovers, by comparing compression cost, imputation power and language distances induced by the imputation. The models take only the etymological data set as input, and require no further linguistic assumptions. In this regard, they is as objective as possible, given the data. The data set itself, of course, may be highly subjective and questionable.
The objectivity of models given the data now opens new possibilities for comparing entire data sets. For example, we can begin to compare the Finnish and Estonian datasets in SSA vs. Star-Ling, although the data sets have quite different characteristics, e.g., different size-3200 vs. 800 word pairs, respectively-and the comparison is done impartially, relying solely on the data provided. Another direct consequence of the presented methods is that they enable us to quantify uncertainty of entries in the corpus of etymological data. For example, for a given entry x in language L 1 , we can compute exactly the probability that x would be imputed by any of the models, trained on all the remaining data from L 1 plus any other set of languages in the family. This can be applied equally to any entry, in particular to entries marked dubious by the database creators.
We can use this method to approach the question of comparison of "competing" etymological datasets. The cost of an optimal alignment obtained over a given data set serves as a measure of its internal consistency.
We are currently working to combine the context model with 3-and higher-dimensional models, and to extend these models to perform diachronic imputation, i.e., reconstruction of protoforms. We also intend to test the models on databases of other language families.