Computing Geodesic Paths Encoding a Curvature Prior for Curvilinear Structure Tracking

The celebrated Euler-Mumford elastica model (Mumford, 1994; Nitzberg, 1990) characterizes smooth curves which minimize a functional featuring a length term and a bending energy term, in other words a curvature penalization. Let $\kappa:[0,1]\to\mathbb{R}$ be the curvature of a smooth curve $\gamma:[0,1]\to\Omega$ with non-vanishing velocity, where $\Omega\subset\mathbb{R}^{2}$ is a bounded domain. Given a data-driven cost function $\tilde{\psi}:\Omega\times\mathbb{R}^{2}\to]0,\infty[$, the energy functional considered in the classical Euler-Mumford elastica model reads

where $\gamma^{\prime}$ is the first-order derivative of the curve $\gamma$, and $\|\cdot\|$ denotes the standard Euclidean norm over the space $\mathbb{R}^{2}$. The constant parameter $\beta>0$ controls the relative importance of the bending energy term $\kappa(u)^{2}$, which penalizes curvature, w.r.t. the penalization of path length. By construction, Eq. 1 is invariant under smooth increasing re-parametrizations of the curve $\gamma$, and the second argument of $\tilde{\psi}$ is the unit tangent vector to the path. In the special case where $\tilde{\psi}$ is constant, the curves $\gamma$ for which Eq. 1 is extremal correspond to the rest positions of elastic rods of suitable length (with their endpoints clamped), however this physical interpretation is lost for non-constant $\tilde{\psi}$.

A widely used strategy for addressing the classical Euler-Mumford elastica optimal curve problem is to lift a planar smooth curve $\gamma$ into the higher dimensional configuration space $\mathbb{M}:=\Omega\times\mathbb{S}^{1}$ of positions and orientations, see (Chen et al, 2017; Duits et al, 2018; Chambolle and Pock, 2019), where $\mathbb{S}^{1}:=\mathbb{R}/2\pi\mathbb{Z}$ denotes the angular domain. In other words $\mathbb{S}^{1}=[0,2\pi[$ with periodic boundary conditions, and we can parametrize the unit circle via the unit vector $\dot{\mathrm{n}}(\theta):=(\cos\theta,\sin\theta)$ for any angle $\theta\in\mathbb{S}^{1}$. Furthermore, we denote by $\mathbb{E}:=\mathbb{R}^{2}\times\mathbb{R}$ the tangent vector space to $\mathbb{M}$, and by $\mathbb{E}^{*}$ the dual vector space of $\mathbb{E}$. Recall that $\gamma:[0,1]\to\Omega$ is a smooth curve with non-vanishing velocity. Following the literature (Chen et al, 2017; Mirebeau, 2018; Duits et al, 2018), an orientation-lifted curve $\bm{\rho}=(\gamma,\eta):[0,1]\to\mathbb{M}$, whose physical projection is the given planar curve $\gamma$, can be constructed by introducing an angular function $\eta:[0,1]\to\mathbb{S}^{1}$ such that

This leads to an alternative representation of the curvature

Inserting this expression of the curvature $\kappa$ into the energy defined in Eq. 1 yields an equivalent formula, measured along an orientation-lifted curve $\bm{\rho}=(\gamma,\eta)$ whose components $\gamma$ and $\eta$ obey Eq. 2:

where $\psi:\mathbb{M}\to]0,\infty[$ is a positive cost function defined by

For the task of curvilinear structure tracking, as considered in this work, the value of $\psi(x,\theta)$ should be relatively low in case the physical position $x$ is inside a curvilinear structure and the direction $\dot{\mathrm{n}}(\theta)$ matches its local geometry at the position $x$, i.e. $\dot{\mathrm{n}}(\theta)$ is approximately tangent to its centerline at $x$. Fig. 2 shows an example for the cost function $\psi$ derived from an image that consists of curvilinear structures. The construction and the numerical computation of $\psi$ is presented in Eq. 25 below, and relies on a tool known as the orientation score.

The energy $\mathcal{L}(\bm{\rho})$ formulated in Eq. 4 can be expressed in terms of a degenerate geodesic metric (Chen et al, 2017; Mirebeau, 2018), denoted by $\mathcal{F}:\mathbb{M}\times\mathbb{E}\to[0,+\infty]$ and referred to as the elastica metric. For any point $\mathbf{x}=(x,\theta)$ and any non-zero vector $\dot{\mathbf{x}}=(\dot{x},\dot{\theta})$ this elastica metric is defined as

Note that this metric is not symmetric, since $\mathcal{F}(\mathbf{x},\dot{\mathbf{x}})\neq\mathcal{F}(\mathbf{x},-\dot{\mathbf{x}})$ in general in view of the non-linear constraint $\dot{x}=\dot{\mathrm{n}}(\theta)\|\dot{x}\|$, and thus the corresponding distance defined below (see Eq. 11) is likewise non-symmetric. Such asymmetric objects are often referred to as quasi-metrics and quasi-distances, but the prefix “quasi-” is dropped in this paper for the sake of readability.

Several mathematical objects can be attached to such a metric. The control sets $\mathcal{B}(\mathbf{x})\subset\mathbb{E}$, for all points $\mathbf{x}\in\mathbb{M}$, allow visualizing the geometry, see Fig. 3, and are defined as

The Hamiltonian $\mathcal{H}:\mathbb{M}\times\mathbb{E}^{*}\to[0,\infty[$ is involved in the numerical computation of optimal geodesic paths through the HJB PDE, as formulated in next subsection, and is defined for any co-vector $\hat{\mathbf{x}}\in\mathbb{E}^{*}$ as

In the special case of the elastica metric, one can use Eq. 5 to obtain the explicit expressions of the control sets and Hamiltonian, see (Mirebeau, 2018) for the details of this computation. Specifically, for any point $\mathbf{x}=(x,\theta)\in\mathbb{M}$ one has

For any co-vector $\hat{\mathbf{x}}=(\hat{x},\hat{\theta})\in\mathbb{E}^{*}$ one obtains

The following equivalent expression of the Hamiltonian, in integral form, is established in Mirebeau (2018)

where $\langle\cdot,\cdot\rangle_{+}:=\max\{0,\langle\cdot,\cdot\rangle\}$ is the positive part of the Euclidean scalar product $\langle\cdot,\cdot\rangle$. Here and below we denote the vector ${\dot{\mathbf{n}}}(\theta,a):=(\dot{\mathrm{n}}(\theta),a)\in\mathbb{E}$ for any angle $\theta\in\mathbb{S}^{1}$ and any scalar $a\in\mathbb{R}$, and observe that in particular $\langle\hat{x},\dot{\mathrm{n}}(\theta)\rangle=\langle\hat{\mathbf{x}},{\dot{\mathbf{n}}}(\theta,0)\rangle$.

Consider a fixed source point $\mathbf{s}\in\mathbb{M}$ and an arbitrary target point $\mathbf{x}\in\mathbb{M}$, thus defining initial and final planar endpoints and tangent path directions. The classical Euler-Mumford elastica geodesic model (Chen et al, 2017; Mirebeau, 2018) aims to track a minimal geodesic path $\mathcal{G}_{\mathbf{s},\mathbf{x}}:[0,1]\to\mathbb{M}$ which links from the source point $\mathbf{s}$ to the target point $\mathbf{x}$ (i.e. $\mathcal{G}_{\mathbf{s},\mathbf{x}}(0)=\mathbf{s}$, $\mathcal{G}_{\mathbf{s},\mathbf{x}}(1)=\mathbf{x}$) and minimizes the energy $\mathcal{L}$ formulated in Eq. 4. For that purpose, we first consider the set $\operatorname{Lip}([0,1],\mathbb{M})$ of all Lipschitz continuous curves $\bm{\rho}:[0,1]\to\mathbb{M}$ as the search space of geodesic paths. The minimization of the energy $\mathcal{L}$ yields a geodesic distance map $\mathcal{D}_{\mathbf{s}}:\mathbb{M}\to[0,\infty)$, defined as

The geodesic distance map $\mathcal{D}_{\mathbf{s}}$ is a potentially discontinuous viscosity solution to the first-order static HJB PDE (Bardi and Capuzzo-Dolcetta, 2008)

where $d\mathcal{D}_{\mathbf{s}}$ denotes the differential of the geodesic distance map $\mathcal{D}_{\mathbf{s}}$. Outflow boundary conditions are applied on $\partial\mathbb{M}$.

A reversely parametrized minimal geodesic path $\tilde{\mathcal{G}}$, starting from a given target point $\mathbf{x}$ and going back to the source point $\mathbf{s}$, can be obtained by solving the following gradient descent ordinary differential equation (ODE), known as the geodesic backtracking ODE

with initial condition $\tilde{\mathcal{G}}(0)=\mathbf{x}$, and where $\partial_{2}\mathcal{H}:=\partial\mathcal{H}(\mathbf{x},\hat{\mathbf{x}})/\partial\hat{\mathbf{x}}$. By construction, this path is well-defined up to the minimal arrival time $T:=\mathcal{D}_{\mathbf{s}}(\mathbf{x})$, and $\tilde{\mathcal{G}}(T)=\mathbf{s}$ is the source point. The reparametrization $\mathcal{G}_{\mathbf{s},\mathbf{x}}(u):=\tilde{\mathcal{G}}(T(1-u))$ yields a minimal geodesic path $\mathcal{G}_{\mathbf{s},\mathbf{x}}:[0,1]\to\mathbb{M}$ such that $\mathcal{G}_{\mathbf{s},\mathbf{x}}(0)=\mathbf{s}$ and $\mathcal{G}_{\mathbf{s},\mathbf{x}}(1)=\mathbf{x}$.

We begin by formulating a new curvature penalization term, which involves a scalar-valued function $\omega:\mathbb{M}\to\mathbb{R}$. For a smooth curve $\gamma$ with curvature $\kappa$, we consider

where the angle $\eta(u)\in\mathbb{S}^{1}$ reflects the tangent path direction for all $u\in[0,1]$, see Eq. 2. In this model, the function $\omega$ is referred to as the curvature prior map. For the purposes of curvilinear structure tracking, as considered in this paper the value $\omega(x,\theta)$ may be defined as the curvature of a reference curve passing nearby the position $x$, with a tangent orientation close to $\theta$, and obtained using an independent skeletonization procedure, see Eq. 24.

We consider the variant of the Euler-Mumford elastic energy involving the curvature prior map $\omega$ and defined as:

For an orientation-lifted curve $\bm{\rho}=(\gamma,\eta)$ satisfying Eq. 2, the corresponding quantity is defined as

The energy $\mathfrak{L}$ leads to a new metric $\mathfrak{F}:\mathbb{M}\times\mathbb{E}\to[0,\infty]$ (Mirebeau and Portegies, 2019)

where $\mathbf{x}=(x,\theta)\in\mathbb{M}$, and where $\dot{\mathbf{x}}=(\dot{x},\dot{\theta})\in\mathbb{E}$, with the convention $\mathfrak{F}(\mathbf{x},\mathbf{0})=0$. The energy $\mathfrak{L}(\bm{\rho})$ can therefore be reformulated as

Path globally minimizing this energy can be obtained, similarly to the classical case and following the general principles of optimal control, by applying the geodesic backtracking ODE to the viscosity solution of a HJB PDE, both involving a suitably modified Hamiltonian described in the next section.

In order to describe our variant of the classical elastica model with curvature prior, we briefly adopt a more general point of view. A (sub-Finslerian) metric on the domain $\mathbb{M}$ is a mapping $\mathcal{F}:\mathbb{M}\times\mathbb{E}\to[0,\infty]$, obeying $\mathcal{F}(\mathbf{x},\lambda\dot{\mathbf{x}})=\lambda\mathcal{F}(\mathbf{x},\dot{\mathbf{x}})$ for any $\mathbf{x}\in\mathbb{M}$, $\dot{\mathbf{x}}\in\mathbb{E}$, $\lambda>0$, satisfying $\mathcal{F}(\mathbf{x},\mathbf{0})=0$, and such that the control sets $\mathcal{B}(\mathbf{x})$ defined by Eq. 6 are non-empty, convex, compact, and depend continuously on the base point $\mathbf{x}\in\mathbb{M}$ w.r.t. the Hausdorff distance (Chen et al, 2017). Some control sets obeying these properties are illustrated in the right column of Fig. 3, where the corresponding Hamiltonian $\mathcal{H}$ is defined by Eq. (7).

The integral form of the Hamiltonian $\mathfrak{H}$, formulated in Eq. 19 can be approximated using the $L$-point Fejer quadrature rule for any positive integer $L$, following (Mirebeau, 2018). This yields

where we denoted by $\mu_{l}>0$ the weights of the Fejer quadrature rule at the angles $\varphi_{l}=(2l-1-L)\pi/(2L)\in[-\pi/2,\pi/2]$, for any $1\leq l\leq L$. This approximation plays a crucial role in the HFM method, as presented in next section.

Centerlines serve as an effective descriptor for curvilinear structures, by which geometric features such as their tangent direction and curvature properties can be estimated. A broad variety of approaches have been introduced to handle the task of extracting centerlines from curvilinear structure networks which usually include complex junctions, e.g. (Kaul et al, 2012; Cetin et al, 2012; Türetken et al, 2016). Among those methods, the skeletonization schemes, adopted in this work, provide a powerful way for producing centerlines as the skeletons of the pre-segmented curvilinear structures (Kimmel et al, 1995; Siddiqi et al, 2002; Lam et al, 1992).

Let $\Xi:\Omega\to\{0,1\}$ represent the binary segmentation of curvilinear structures contained in an image, where $\Xi(x)=1$ implies that the position $x$ belongs to a curvilinear structure, while $\Xi(x)=0$ indicates a background position $x$. We compute the skeleton or medial axis of the elongated shapes $\{x\in\Omega\mid\Xi(x)=1\}$, and remove all the junction points of these skeleton structures, yielding $J$ disjoint skeleton segments, where $J$ is a positive integer. After applying a curve smoothing procedure to regularize these skeleton segments, we obtain a family of smooth curves $\rho_{j}:[0,1]\to\Omega$, indexed by $1\leq j\leq J$, as the parametrization of the regularized skeleton segments. In what follows we refer to these curves $\rho_{j}$ as candidate centerlines.

In practice, we rely on morphological filters (Lam et al, 1992), which benefit from a low computational complexity, to generate skeletons from pre-segmentation data. Fig. 4 illustrates an example for constructing the candidate centerlines $\{\rho_{j}\}_{1\leq j\leq J}$ using a synthetic image. Fig. 4b visualizes the corresponding binary segmentation map $\Xi$. In Fig. 4c, the candidate centerlines are illustrated using different colors, and arrows indicate their directions. The curvature of these candidate centerlines is depicted in Fig. 4d.

Note that a curvilinear structure of interest usually corresponds to the concatenation of several candidate centerlines, together with the crucial curves connecting pairs of adjacent centerlines. For this reason, the skeletonization process is not an appropriate solution to the curvilinear structure tracking problem considered, but rather an efficient descriptor for extracting the geometric features of parts of the curvilinear structures.

We propose to estimate the curvature prior map $\omega$ using the computed candidate centerlines $\rho_{j}$, whose curvature and unit normal are respectively denoted by $\kappa_{j}:[0,1]\to]-\infty,\infty[$ and $\mathcal{N}_{j}:[0,1]\to\mathbb{R}^{2}$, for $1\leq j\leq J$.

Let $0<U<\min_{j}\{1/\|\kappa_{j}\|_{\infty}\}$ be a bounded constant. Consider the mapping $\Psi_{j}:[0,1]\times\,[-U,U]\,\to\mathbb{R}^{2}$, depending on the curve parameter and the deviation from $\rho_{j}$, defined as

for any $u\in[0,1]$ and $\lambda\in[-U,U]$. The bounding constraint on the parameter $U$ implies that the mappings $\Psi_{j}$ are injective, and each of them parametrizes a tubular neighborhood $\mathfrak{T}_{j}\subset\Omega$ of width $U$ surrounding the candidate centerline $\rho_{j}$

However, such a tubular neighborhood $\mathfrak{T}_{i}$ may intersect with a distinct one $\mathfrak{T}_{j}$, i.e. $\mathfrak{T}_{i}\cap\mathfrak{T}_{j}\neq\emptyset$, when their centerlines $\rho_{i}$ and $\rho_{j}$ are close to each other. In order to tackle this problem, we partition the domain $\Omega$ into a series of disjoint regions $\mathcal{V}_{j}\subset\Omega$ for $1\leq j\leq J$

where $d_{\rho_{j}}(x)$ denotes the Euclidean distance between $x$ and the candidate centerline $\rho_{j}$. These regions are also known as the Voronoi cells. Accordingly, for each $\rho_{j}$ a modified tubular neighborhood $\mathcal{T}_{j}\subset\Omega$ can be constructed by

which satisfies $\mathcal{T}_{j}\cap\mathcal{T}_{i}=\emptyset$ for any $i\neq j$. In practice, the construction of the Voronoi cells $\mathcal{V}_{j}$ can be efficiently implemented by propagating an unsigned Euclidean distance map simultaneously from all the candidate centerlines Saye and Sethian (2011, 2012).

We consider two scalar-valued functions $\phi:\cup_{j}\mathcal{T}_{j}\to\mathbb{R}$ and $\vartheta:\cup_{j}\mathcal{T}_{j}\to\mathbb{S}^{1}$ using all the candidate centerlines $\rho_{j}$, $1\leq j\leq J$, such that for any physical position $x=\Psi_{j}(u,\lambda)\in\mathcal{T}_{j}$, one respectively has

The function $\phi$ extends the curvature of the candidate centerlines $\rho_{j}$ to their respective neighborhood regions $\mathcal{T}_{j}$ of width $U$. Furthermore, for a physical position $x\in\mathcal{T}_{j}$, the angles $\vartheta(x)$ and $\vartheta(x)+\pi$ correspond to the two possible orientations that the elongated structure should have at the position $x$. With these definitions, the curvature prior map $\omega$ can be generated in terms of the estimated functions $\phi$ and $\vartheta$. More precisely, we define the curvature prior map $\omega$ for any point $\mathbf{x}=(x,\theta)\in\mathbb{M}$ as follows

where $\operatorname{sign}(a)$ stands for the sign of a scalar value $a\in\mathbb{R}$. In the first line of Eq. 24, for any physical position $x\in\cup_{j}\mathcal{T}_{j}$ we set $\omega(x,\theta)=\omega(x,\vartheta(x))=\phi(x)$ if the angular position $\theta$ is closer to $\vartheta(x)$ in the sense of the Euclidean scalar product of the vectors $\dot{\mathrm{n}}(\theta)$ and $\dot{\mathrm{n}}(\vartheta(x))$, and $\omega(x,\theta)=\omega(x,\vartheta(x)+\pi)=-\phi(x)$ if $\theta$ is closer to $\vartheta(x)+\pi$.

The multi-orientation data-driven cost function $\psi(x,\theta)$ should have low values when both (i) the physical position $x$ is close to the centerline of a curvilinear structure, and (ii) the direction $\dot{\mathrm{n}}(\theta)$ closely aligns with the tangent of this centerline at the physical position $x$. For that purpose and following (Chen et al, 2017; Duits et al, 2018; Bekkers et al, 2015a), we build the cost function $\psi$ for any point $\mathbf{x}=(x,\theta)\in\mathbb{M}$ as follows

where $\alpha>0$ is a positive weighting parameter, and where $g:\mathbb{M}\to[0,\infty[$ is an orientation score map usually taken as an efficient tool in image analysis (Franken and Duits, 2009; Hummel and Zucker, 1983; Bekkers et al, 2014; Parent and Zucker, 1989). This map carries the curvilinear appearance and anisotropy features, and obeys $g(x,\theta)=g(x,\theta+\pi)$, $\forall(x,\theta)\in\mathbb{M}$. In this work, the steerable optimally oriented flux filter (Law and Chung, 2008) is chosen as the curvilinear feature extractor in order to compute the orientation score $g$. We refer to (Chen et al, 2017) for more details on the computation of the orientation score map $g$ from image data involving curvilinear structures.

The curvilinear structure tracking method can be implemented in an interactive manner, requiring a source point $\mathbf{s}=(s,\theta_{s})\in\mathbb{M}$ and a target point $\mathbf{y}=(y,\theta_{y})\in\mathbb{M}$ as initialization. In our model, we assume that the physical positions $s$ and $y$ are provided by the user. Furthermore, the corresponding angles $\theta_{s}$ and $\theta_{y}$ can be either given by the user in complex scenarios, or be estimated from the cost function $\psi$ as

In addition to the points $\mathbf{s}$ and $\mathbf{y}$ as described above, we generate another pair of source and target points $\tilde{\mathbf{s}}=\big{(}s,\theta_{s}+\pi\big{)}$ and $\tilde{\mathbf{y}}=\big{(}y,\theta_{y}+\pi\big{)}$. In other words, $\mathbf{s}$ and $\tilde{\mathbf{s}}$ (resp. $\mathbf{y}$ and $\tilde{\mathbf{y}}$) share the same physical position but different angular positions. In this case, the goal is to find a minimal geodesic path from the source point set $\mathfrak{s}=\{\mathbf{s},\tilde{\mathbf{s}}\}$ to the target point set $\mathfrak{y}=\{\mathbf{y},\tilde{\mathbf{y}}\}$, as considered in (Chen et al, 2017). For that purpose, a geodesic distance map $\mathcal{D}_{\mathfrak{s}}$ is estimated by solving the HJB PDE

for all $\mathbf{x}\in\mathbb{M}\setminus\mathfrak{s}$, with the boundary condition $\mathcal{D}_{\mathfrak{s}}(\mathbf{x})=0$ for any point $\mathbf{x}\in\mathfrak{s}$, which is a straightforward generalization of the single source case described in Eq. 28. The geodesic of interest is backtracked, as described in Eq. 13, from the point of the target set $\mathfrak{y}$ whose geodesic distance value is smallest, and by the nature of the viscosity solution the backtracking procedure ends at the point of the seed set $\mathfrak{s}$ which is closest w.r.t. the metric. For computational efficiency, the procedure for estimating the respective geodesic distances can be terminated immediately once either point of the target set $\mathfrak{y}$ is reached by the fast marching front, and this optimization is valid since the fast marching front advances monotonically. Note that this geodesic path tracking procedure allows the user not to distinguish between the source and target points, which is convenient in practice.

In Fig. 7, we compare results obtained with the classical elastica model, and with the curvature prior elastica model, on five synthetic images. A common feature for the curvilinear structures in those images is that each of them involves strongly bent segments. The first column of Fig. 7 shows the initialization information for each synthetic image, where the blue and red dots together with the arrows stand for the source and target points respectively used to define the sets $\mathfrak{s}$ and $\mathfrak{y}$. Columns $2$ and $3$ illustrate the physical projections of the computed geodesic paths from the classical elastica model and the curvature prior elastica model, respectively. More specifically, in column $2$ of rows $1$ to $3$, we can observe that the physical projections of the geodesic paths from the classical elastica model are trapped into a shortcut problem, i.e. the projected paths travel along non-target regions. The reason is that the “correct” path features strongly curved segments whose squared curvature accumulate and lead to a high value of the classical elastica energy, hence it is not selected by this model. In contrast, in column $3$ of the same rows, the curvature prior elastica model indeed leads to accurate extraction of the target structures, even through the two endpoints of the structures in rows $2$ and $3$ are very close to each other. In rows $4$ and $5$, each of the curvilinear structures is interrupted by an almost straight segment with slightly stronger appearance features, i.e. with higher values of the orientation score along this segment, yielding challenging crossing structures. In these tests, the physical projections of the geodesic paths from the classical elastica model pass through the interrupting segments, yielding a short branches combination problem (chen2018minimal), i.e. the obtained geodesic paths travel different segments not belonging to the targets. In contrast, the curvature prior elastica model produces satisfactory results which are able to accurately delineate the target curvilinear structures, as shown in column $3$.

Fundus vessel tracking is a task of fundamental importance in retinal imaging and the relevant disease diagnosis. In Fig. 8, we compare the curvature prior elastica model to the classical elastica model in $3$ patches of retinal fundus images, chosen from the IOSTAR dataset (Zhang et al, 2016). In this experiment, the goal is to extract artery blood vessels, given two endpoints of each target vessel, as shown in the left column of Fig. 8. In a typical retinal fundus image, the major artery vessels usually exhibit weaker appearance features than those of the nearby vein vessels, and often cross over vein vessels yielding junction structures. An additional difficulty is that the target artery vessels also consist of segments with strong tortuosity. For the classical Euler-Mumford elastica model, one needs to assign a low importance weight to the curvature term, in order to avoid the underlying shortcut problem at those segments. However, this in turn encourages the computed paths to pass through the stronger neighboring vein vessels (whose cost $\psi$ is lower), thus increasing the risk of the short branches combination problem, see the center column of Fig. 8. In contrast, the minimization of the energy (16) defining the curvature prior elastica model favors optimal geodesic paths whose curvature is close to $\omega$, allowing to handle scenarios where challenging structures such as crossings and strongly bent segments appear, see the right column of Fig. 8.

Another practical application for geodesic models is road tracking from aerial images, as demonstrated in Fig. 9. In this experiment, the centerlines of the objective roads have high Euclidean length, and simultaneously feature segments with low turning radii, as shown in Fig. 9a. We can see that the classical elastica model generates tracking results which suffer from a shortcut problem, see the center column of Fig. 9. In contrast, accurate results derived from the curvature prior elastica model are observed in the right column of Fig. 9, thanks to the use of the curvature prior term which alleviates difficulties associated to the complex scenarios.

We perform a quantitative evaluation of the performance and stability of the classical elastica model and of the curvature prior elastica model on synthetic images, when one varies the parameters $\beta$ and $\alpha$ involved in these models, see Eqs. (16) and (25). These parameters determine the balance between the image data and the curvature regularization. The test data are generated by incorporating different levels of additive Gaussian noise to the original binary images of those shown in Fig. 7. More specifically, we add zero-mean Gaussian noise with $16$ levels of normalized variance values between $0$ and $0.15$ to each clean image, producing $16\times 5=80$ test images. Both geodesic models are configured by choosing parameters $\beta\in\{4.5,\,6\}$ and $\alpha\in\{3,4,5\}$, thus yielding $6$ runs per test image. In each test, the accuracy of the curvilinear structure tracking is measured using the Jaccard score between two bounded regions $A,B\subset\mathbb{R}^{2}$, defined as

where $|A|$ represents the area of the region $A$. In our case, the regions $A$ and $B$ are defined as tubular neighborhoods with a fixed radius of $6$ grid points along the ground truth centerline and the physical projection of the computed geodesic path, respectively. Table 1 shows the quantitative results comparing the classical elastica model and the curvature prior elastica model, in terms of the Jaccard score between the ground truth centerlines and the physical projections of geodesic paths, see (27). In each row of Table 1, the average (Ave.) and standard deviation (Std.) values of $\operatorname{JS}$, which are derived from $480$ runs w.r.t. different values of the parameters $\alpha$ and $\beta$. Overall, the results from the classical Euler-Mumford elastica model in general exhibit low average values of $\operatorname{JS}$, and large standard deviation, reflecting the fact that various shortcut and short branches combination problems occur in most of the tests. In contrast, we observe much higher average values of $\operatorname{JS}$ using the curvature prior elastica model, and a smaller standard deviation, illustrating its ability to accurately and robustly extract curvilinear structures in the presence of strongly bent segments and junctions.