Global Minimum for a Finsler Elastica Minimal Path Approach

# Global Minimum for a Finsler Elastica Minimal Path Approach

Da Chen Da ChenLaurent D. Cohen University Paris Dauphine, PSL Research University
CNRS, UMR 7534, CEREMADE, 75016 Paris, France
Jean-Marie Mirebeau
Laboratoire de mathématiques d’Orsay, Université Paris-Sud, CNRS,
Université Paris-Saclay, 91405 Orsay, France
44email: jean-marie.mirebeau@math.u-psud.fr
Jean-Marie Mirebeau Da ChenLaurent D. Cohen University Paris Dauphine, PSL Research University
CNRS, UMR 7534, CEREMADE, 75016 Paris, France
Jean-Marie Mirebeau
Laboratoire de mathématiques d’Orsay, Université Paris-Sud, CNRS,
Université Paris-Saclay, 91405 Orsay, France
44email: jean-marie.mirebeau@math.u-psud.fr
Laurent D. Cohen Da ChenLaurent D. Cohen University Paris Dauphine, PSL Research University
CNRS, UMR 7534, CEREMADE, 75016 Paris, France
Jean-Marie Mirebeau
Laboratoire de mathématiques d’Orsay, Université Paris-Sud, CNRS,
Université Paris-Saclay, 91405 Orsay, France
44email: jean-marie.mirebeau@math.u-psud.fr
Received: 19-December-2015/ Accepted: 17-November-2016
###### Abstract

In this paper, we propose a novel curvature penalized minimal path model via an orientation-lifted Finsler metric and the Euler elastica curve. The original minimal path model computes the globally minimal geodesic by solving an Eikonal partial differential equation (PDE). Essentially, this first-order model is unable to penalize curvature which is related to the path rigidity property in the classical active contour models. To solve this problem, we present an Eikonal PDE-based Finsler elastica minimal path approach to address the curvature-penalized geodesic energy minimization problem. We were successful at adding the curvature penalization to the classical geodesic energy (Caselles et al, 1997; Cohen and Kimmel, 1997). The basic idea of this work is to interpret the Euler elastica bending energy via a novel Finsler elastica metric that embeds a curvature penalty. This metric is non-Riemannian, anisotropic and asymmetric, and is defined over an orientation-lifted space by adding to the image domain the orientation as an extra space dimension. Based on this orientation lifting, the proposed minimal path model can benefit from both the curvature and orientation of the paths. Thanks to the fast marching method, the global minimum of the curvature-penalized geodesic energy can be computed efficiently.

We introduce two anisotropic image data-driven speed functions that are computed by steerable filters. Based on these orientation-dependent speed functions, we can apply the proposed Finsler elastica minimal path model to the applications of closed contour detection, perceptual grouping and tubular structure extraction. Numerical experiments on both synthetic and real images show that these applications of the proposed model indeed obtain promising results.

###### Keywords:
Minimal Path Eikonal Equation Curvature Penalty Anisotropic Fast Marching Method Image Segmentation Tubular Structure Extraction
journal: International Journal of Computer Vision

## 1 Introduction

Snakes or active contour models have been studied considerably, and used for object segmentation and feature extraction for almost three decades, since the pioneering work of the snakes model proposed by Kass et al (1988). A snake is a parametrized curve (locally) that minimizes the energy:

 E(Γ)=∫10(w1∥Γ′(t)∥2+w2∥Γ′′(t)∥2+P(Γ(t)))dt

with appropriate boundary conditions at the endpoints and . and are the first- and second-order derivatives of the curve respectively. The positive constants and relate to the elasticity and rigidity of the curve and, hence, weight its internal force. This approach models contours as curves locally minimizing an objective energy functional that consists of an internal and an external force. The internal force terms depend on the first- and second-order derivatives of the curves (snakes), and, respectively, account for a prior of small length and of low curvature of the contours. The external force is derived from a potential function , which depends on image features such as the gradient magnitude, and it is designed to attract the active contours or snakes to the image features of interest such as object boundaries.

The drawbacks of the classical active contours or snakes model (Kass et al, 1988) are its sensitivity to initialization, the difficulty of handling topological changes, and the difficulty of minimizing the strongly non-convex path energy as discussed by Cohen and Kimmel (1997). Regarding initialization, the active contours model requires an initial guess that is close to the desired image features, and preferably enclosing them because energy minimization tends to shorten the snakes. The introduction of an expanding balloon force allows the model to be less demanding on the initial curve given inside the objective region (Cohen, 1991). The issue of topology changes leads, on the other hand, to the development of active contour methods, which represent object boundaries as the zero level set of the solution to a PDE (Osher and Sethian, 1988; Caselles et al, 1993; Malladi et al, 1994; Caselles et al, 1997; Yezzi et al, 1997).

The difficulty of minimizing the non-convex snakes energy (Kass et al, 1988) leads to important practical problems because the curve optimization procedure often becomes stuck at a local minimum of the energy functional, which makes the results sensitive to curve initialization and image noise. This limitation is still the case for the level set approach on geodesic active contours (Malladi et al, 1995; Caselles et al, 1997). To address this local minimum sensitivity issue, Cohen and Kimmel (1997) proposed an Eikonal PDE-based minimal path model to find the global minimum of the geodesic active contours energy (Caselles et al, 1997), in which the penalty associated to the second-order derivative of the curve was removed from the snakes energy. Thanks to this simplification, a fast, reliable and globally optimal numerical method allows to find the energy minimizing curve with prescribed endpoints, namely, the fast marching method (Sethian, 1999), which is based on the formalism of viscosity solutions to Eikonal PDE. These mathematical and algorithmic guarantees of Cohen and Kimmel’s minimal path model (Cohen and Kimmel, 1997) have important practical consequences, leading to various approaches for image analysis and medical imaging (Peyré et al, 2010; Cohen, 2001).

In the basic formulations of the minimal paths-based interactive image segmentation models (Appleton and Talbot, 2005; Appia and Yezzi, 2011; Mille et al, 2014), the common proposal is that the object boundaries can be delineated by a set of minimal paths constrained by user-provided points. In (Li and Yezzi, 2007; Benmansour and Cohen, 2011), vessels were extracted under the form of minimal paths over the radius-lifted space. Therefore, each extracted minimal path consists of both the centerline positions and the corresponding thickness values of a vessel.

To reduce the user intervention, Benmansour and Cohen (2009) proposed a growing minimal path model for object segmentation with unknown endpoints. This model can recursively detect keypoints, each of which can be used as a new source point for the fast marching front propagation. Thus this model requires only one user-provided point to start the keypoints detection procedure. Kaul et al (2012) applied the growing minimal path model to detect complicated curves with arbitrary topologies and developed criteria to stop the growth of the minimal paths. Rouchdy and Cohen (2013) proposed a geodesic voting model for vessel tree extraction by a voting score map that is constructed from a set of geodesics with a common source point.

Recently, improvements of the minimal path model have been devoted to extend the isotropic Riemannian metrics to the more general anisotropic Riemannian metrics by taking into account the orientation of the curves (Bougleux et al, 2008; Jbabdi et al, 2008; Benmansour and Cohen, 2011). Such orientation enhancement can solve some shortcuts problems suffered by the isotropic minimal path models (Cohen and Kimmel, 1997; Li and Yezzi, 2007). Kimmel and Sethian (2001) designed an orientation-lifted Riemannian metric for the application of path planning, providing an alternative way to take advantage of the orientation information. This isotropic orientation-lifted Riemannian metric (Kimmel and Sethian, 2001) was built over a higher dimensional domain by adding an extra orientation space to the image domain.

Bekkers et al (2015) considered a data-driven extension of the sub-Riemannian metric on SE(2), which shows that the sub-Riemannian structure outperforms the isotropic Riemannian structures on SE(2) in retinal vessel tracking. The numerical solver used by Bekkers et al (2015) is based on a PDE approach with an upwind discretization and iterative evolution scheme, which requires expensive computation time. To solve this problem, Sanguinetti et al (2015) used the fast marching method (Mirebeau, 2014a) as the Eikonal solver to track the sub-Riemannian geodesics. The sub-Riemannian geodesic model (Petitot, 2003) reintroduced curvature penalization to the geodesic energy, similar to the Euler elastica bending energy (Nitzberg and Mumford, 1990) considered in this paper, yet it differs in two ways: firstly, the Euler elastica bending energy involves the squared curvature, a stronger penalization than Petitot’s sub-Riemannian geodesic energy which is roughly linear in the curvature. Secondly, minimal geodesics for Petitot’s sub-Riemannian model occasionally feature cusps, which sometimes may not be desirable for the applications of interest. In contrast, Euler elastica curves can avoid such cusps.

Schoenemann et al (2012) proposed a model to address the problems of using curvature regularization for region-based image segmentation by a graph-based combinational optimization method. This curvature regularization model can find a solution that corresponds to the globally optimal segmentation which is initialization-independent, and has proven to obtain promising segmentation and inpainting results especially for objects with long and thin structures. Ulen et al (2015) proposed a curvature and torsion regularized shortest path model for tubular structure segmentation, where the curvature and the torsion were approximately computed by B-Splines. The solution of their energy functional including curvature and torsion penalization terms can be efficiently obtained by using line graphs and Dijkstra’s algorithm (Dijkstra, 1959).

Tai et al (2011) presented an efficient method to solve the minimization problem of the Euler elastica energy, and demonstrated that this fast method can be applied to image denosing, inpainting, and zooming. Other approaches of interest using the curvature penalization include the image segmentation models such as (El-Zehiry and Grady, 2010; Schoenemann et al, 2011; Zhu et al, 2013) and the image inpainting model introduced by Shen et al (2003).

### 1.1 Motivation

In contrast with the classical snakes energy (Kass et al, 1988), Eikonal PDE-based minimal path methods are first-order models, which do not penalize the second-order derivative of a curve, i.e., the curvature, and therefore do not enforce the smoothness of the geodesic, leading sometimes to undesired result as shown in Fig. 1, in which we would like to extract a boundary that is as smooth as possible between the two given points indicated by red and green dots. In Fig. 1a we show the edge saliency map. Figs. 1b and 1c are the minimal paths obtained by using the isotropic Riemannian metric (Cohen and Kimmel, 1997) and the anisotropic Riemannian metric (Bougleux et al, 2008), respectively, both of which are unable to find expected smooth boundaries and suffer from the shortcut problem due to the lack of curvature penalization in these metrics. In contrast, the minimal path model presented in this paper reintroduces the curvature, in the form of weighted Euler elastica curves as studied in (Nitzberg and Mumford, 1990; Mumford, 1994). Therefore, the geodesics extracted by the proposed metric can catch the smooth object boundary, as shown in Fig. 1d, with arrows indicating the corresponding tangents at the given positions.

### 1.2 Contributions

The contribution of this paper is three fold:

1. Firstly, we propose a novel globally minimized minimal path model, namely, the Finsler elastica minimal path model, with curvature penalization and Finsler metric. We establish the connection between the Euler elastica curves and the minimal paths with respect to a Finsler elastica metric. With an adequate numerical implementation, leveraging orientation lifting, asymmetric Finsler metrics and anisotropic fast marching method, the proposed model still allows to find the globally minimizing curves with prescribed points and tangents.

2. As a second contribution, we present the mathematical convergence analysis of the proposed Finsler elastica metrics, and of the associated Finsler elastica minimal paths. Furthermore, we discuss numerical options for geodesic distance and minimal paths computation, and settle for an adaptation of the fast marching method proposed by Mirebeau (2014b).

3. Finally, we provide two types of image data-driven speed functions that are computed by steerable filters. These speed functions are therefore orientation dependent, by which we apply the proposed Finsler elastica minimal path model to interactive closed contour detection, perceptual grouping and tubular structure extraction.

More contributions have been added regarding the original conference paper presented in (Chen et al, 2015), such as the applications of interactive closed contour detection for image segmentation and perceptual grouping, respectively.

### 1.3 Paper Outline

In Section 2 we briefly introduce the existing minimal path models, the concept of Finsler metric, and algorithms for distance computation and path extraction. The relationship between the Euler elastica bending energy and the Finsler elasica metric is analyzed in Section 3. In Section 4 we introduce two data-driven speed functions which are dependent of orientations. The applications of the Finsler elastica minimal path model are presented in Section 5. Experiments and Conclusion are presented in Sections 6 and 7, respectively.

## 2 Background on the Minimal Path Models

### 2.1 Cohen-Kimmel Minimal Path Model

Let () denote the image domain. The classical Cohen-Kimmel minimal path model (Cohen and Kimmel, 1997) was designed to find the global minimum of the geodesic energy which measures the length of Lipschitz paths as follows

 LI(γ)=∫10(w+P(γ(t)))∥γ′(t)∥dt, (1)

where denotes the canonical Euclidean norm. is a positive constant and is a potential function which usually depends on the image gradient magnitudes or intensities as suggested by Cohen and Kimmel (1997).

A geodesic , joining the source point to a point , is a path that globally minimizes the length (1):

 Cs,x=argminγ∈As,x{LI(γ)},

where is the collection of all Lipschitz paths with and .

The minimal action map associated to a source point , is the global minimum of the length (1) among all paths involved in the collection :

 (2)

The Cohen-Kimmel minimal path model considers an isotropic Riemannian metric, which is defined by

 FI(x,u):=√⟨u,MI(x)u⟩,∀x∈Ω,∀u∈Rn, (3)

where denotes the scalar product over and is a symmetric positive definite tensor field which is proportional to the identity matrix :

 MI(x)=(w+P(x))2Id,∀x∈Ω. (4)

The minimal action map defined in (2) satisfies the following Eikonal PDE (Cohen and Kimmel, 1997):

 {∥∇U(x)∥=w+P(x),∀x∈Ω∖{s},U(s)=0, (5)

where is the Euclidean gradient of with respect to the position in the domain.

A geodesic parameterized by its arc-length can be obtained through the gradient descent ordinary differential equation (ODE) till reaching the point :

 {^C′x,s(t)∝−∇U(^Cx,s(t)),^Cx,s(0)=x, (6)

where denotes positive collinearity. Then the geodesic is obtained by reversing the path with and .

### 2.2 Anisotropic Riemannian Metric Extension

Bougleux et al (2008) and Jbabdi et al (2008) extended the isotropic Riemannian metric (3) to the anisotropic Riemannian case invoking a symmetric positive definite tensor field

 MA(x)=n∑i=1Pi(x)vi(x)vTi(x),∀x∈Ω, (7)

where is the dimension of the domain and the vector fields are related to the image data. are the image data-driven functions associated to .

Based on the tensor field  (7), the anisotropic Riemannian metric can be expressed as:

 FA(x,u):=√⟨u,MA(x)u⟩,∀x∈Ω,∀u∈Rn. (8)

The minimal action map associated to the anisotropic Riemannian metric satisfies the Eikonal PDE:

 {∥∇U(x)∥M−1A(x)=1,∀x∈Ω∖{s},U(s)=0, (9)

where we denote that .

The geodesic is obtained by reversing the geodesic which is the solution to the following ODE

 {^C′x,s(t)∝−M−1A(^Cx,s(t))∇U(^Cx,s(t)),^Cx,s(0)=x. (10)

Benmansour and Cohen (2011) proposed an anisotropic radius-lifted Riemannian metic, as a special case of the metric , to address the problem of extracting both the centerline positions and the thickness values of a tubular structure simultaneously.

### 2.3 Isotropic orientation-lifted Riemannian Metric Extension

In order to take into account the local orientation in the image, it is possible to include orientation information in the energy minimization. For this purpose, the image domain can be extended to an orientation-lifted space by product with an abstract orientation space  (Kimmel and Sethian, 2001), i.e., and the problem is to find a geodesic in the new lifted space . Each point in the orientation-lifted path is thus a pair composed of a point in the image domain and an orientation , i.e., .

For any point and any vector , where denotes the orientation variation, the isotropic orientation-lifted Riemannian metric is expressed by

 FIL(¯x,¯u):=1ΦIL(¯x)√∥u∥2+ρ|ν|2, (11)

where is an orientation-dependent speed function relying on the image data and is a constant.

Following the form of (3), the metric defined in (11) can be written by using a tensor field

 (12)

where is a diagonal matrix. Note that can be extended to a positive scalar function that is dependent of the image data. For simplicity, we set to be a positive constant in this paper, as suggested in (Kimmel and Sethian, 2001; Péchaud et al, 2009).

The curve length of an orientation-lifted curve with respect to the isotropic orientation-lifted Riemannian metric can be expressed as

 LIL(γ) =∫10FIL(γ(t),γ′(t))dt =∫101ΦIL(γ(t))√⟨γ′(t),MIL(γ(t))γ′(t)⟩dt =∫101ΦIL(γ(t))√∥Γ′(t)∥2+ρ|θ′(t)|2dt, (13)

where , for all .

The idea of orientation lifting was applied to tubular structure extraction by Péchaud et al (2009) with additional radius lifting, typically accounting for the radius of a tubular structure in the processed image. Orientation lifting often improves the results from (Cohen and Kimmel, 1997; Li and Yezzi, 2007), but suffers from the fact that for the curve length of an orientation-lifted curve , nothing in (2.3) constrains the path direction to align with the orientation vector associated to , where . In other words, this isotropic orientation-lifted Riemannian metric cannot penalize the curvature of the geodesic, a point which is addressed in this paper.

### 2.4 General Minimal Path Model and Finsler Metric

The general minimal path problem is posed on a bounded domain equipped with a metric depending on positions and orientations . This metric defines at each point an asymmetric norm

 Fx(u):=F(x,u). (14)

These norms must be positive whenever , -homogeneous, and obey the triangular inequality. However, in general we allow them to be asymmetric: such that

 Fx(u)≠Fx(−u). (15)

One can measure the curve length of a regular curve with respect to the metric :

 LF(γ)=∫10F(γ(t),γ′(t))dt, (16)

The minimal action map is defined by:

 U(x):=inf{LF(γ);γ∈As,x}, (17)

which is the unique viscosity solution to an Eikonal PDE (Lions, 1982; Sethian and Vladimirsky, 2003):

 {F∗x(∇U(x))=1,∀x∈Ω∖{s},U(s)=0, (18)

where is the dual norm of defined for all by

 F∗x(u):=supv≠0⟨u,v⟩Fx(v). (19)

Based on the definition of the dual norm in (19), the corresponding optimal direction map is then obtained by

 Ψ(x,u):=argmaxv≠0⟨u,v⟩Fx(v),∀x∈Ω,∀u∈Rn. (20)

Again, the geodesic is obtained by reversing the geodesic with and , where is tracked through the following ODE involving the minimal action map and the optimal direction map

 ⎧⎨⎩^C′x,s(t)∝−Ψ(^Cx,s(t),∇U(^Cx,s(t))),^Cx,s(0)=x. (21)

Numerically, the ODE expressed in (21) is solved by using Heunâs or Runge-Kutta’s methods, or more robustly using the numerical method proposed by Mirebeau (2014a).

The metric considered in this paper combines a symmetric part, defined in terms of a symmetric positive definite tensor field , and an asymmetric part involving a vector field :

 F(x,u):=√⟨u,M(x)u⟩−⟨→ω(x),u⟩, (22)

for all and any vector . The fields and should obey the following smallness condition to ensure that the Finsler metric is positive:

 ⟨→ω(x),M−1(x)→ω(x)⟩<1,∀x∈Ω. (23)

Equation (22) defines an anisotropic Finsler metric in general. This is an anisotropic Riemannian metric if the vector field is identically zero, and an isotropic metric if in addition the tensor field is proportional to a diagonal matrix like in (12). Therefore, the general Eikonal PDE (18) and the geodesic back tracking ODE (21) reduce respectively to equations (5) and (6) for the isotropic case, and respectively to equations (9) and (10) for the anisotropic case.

### 2.5 Computation of the Minimal Action Map

In order to estimate the minimal action map , presented in (17) and (18), a discretization grid of the image domain is introduced - or of the extended domain in the case of an orientation-lifted metric. For each point , a small mesh of a neighbourhood of with vertices in is constructed. For example, can be the square formed by the four neighbours of in the classical fast marching method (Sethian, 1999) on a regular orthogonal grid used in (Cohen and Kimmel, 1997). In contrast with Sethian’s classical fast marching method which solves the discrete approximation of the Eikonal PDE itself, an approximation of the action map , with the initial source point , is calculated by solving the fixed point system (Tsitsiklis, 1995):

 {U(x)=ΛU(x),%forallx∈Z∖{s},U(s)=0, (24)

where the involved Hopf-Lax update operator is defined by:

 (25)

where denotes the piecewise linear interpolation operator on the mesh , and lies on the boundary of . The equality , replacing in (24) the Eikonal PDE: of (18), is a discretization of Bellman’s optimality principle, which is similar in spirit to the Tsitsiklis approach (Tsitsiklis, 1995). It reflects the fact that the minimal geodesic , from to , has to cross the mesh boundary at least once at some point ; therefore it is the concatenation of a geodesic from to , which length is approximated by piecewise linear interpolation, and a very short geodesic from to , approximated by a segment of geodesic curve length . The -dimensional fixed point system (24), with the number of grid points, can be solved by a single pass fast marching method as described in Algorithm 1 in Appendix.

The classical fast marching methods (Sethian, 1999; Tsitsiklis, 1995) using the square formed neighbourhood have difficulty to deal with the computation of geodesic distance maps with respect to anisotropic metrics, especially when the anisotropy gets large. An adaptive construction method of such stencils was introduced in (Mirebeau, 2014a) for anisotropic 3D Riemannian metric, and in (Mirebeau, 2014b) for anisotropic 2D Finsler metric, providing that the stencils or mesh at each point or satisfies some geometric acuteness property depending on the local metric . Such adaptive stencils based fast marching methods lead to breakthrough improvements in terms of computation time and accuracy for strongly anisotropic geodesic metrics. When the above mentioned geometric properties do not hold, the fast marching method is in principle not applicable, and slower multiple pass methods must be used instead, such as the Adaptive Gauss Siedel Iteration (AGSI) of Bornemann and Rasch (2006). The present paper involves a 3D Finsler metric (37), for which we constructed stencils by adapting the 2D Finsler metric construction method proposed by Mirebeau (2014b). Although these stencils lack the geometric acuteness condition, we found that the fast marching method still provided good approximations of the paths, while vastly improving computation performance. In Section 3.3, the complexity will be discussed more.

Note that whenever we mention fast marching method in the next sections, we mean the fast marching method with adaptive stencils proposed by Mirebeau (2014b).

## 3 Finsler Elastica Minimal Path Model

In this section, we present the core contribution of this paper: the orientation-lifted Finsler metric embedding curvature penalty term, defined over the orientation-lifted domain , where ) denotes the angle space with periodic boundary conditions and denotes the image domain.

### 3.1 Geodesic Energy Interpretation of the Euler Elastica Bending Energy via a Finsler Metric

The Euler elastica curves were introduced to the field of computer vision by Nitzberg and Mumford (1990) and Mumford (1994). They minimized the following bending energy:

 L(^Γ):=∫L01Φ0(^Γ(s))(1+ακ2(s))ds, (26)

where is a regular curve with non-vanishing velocity, is the curvature of curve , is the Euclidean curve length of , and is the arc-length. Parameter is a constant. is an image data-driven speed function.

For the sake of simplicity, we set , yielding the simplified Euler elastica bending energy

 ℓ(^Γ)=∫L0(1+ακ2(s))ds, (27)

where the general case will be studied in Section 3.4.

The goal of this section is to cast the Euler elastica bending energy  (27) in the form of curve length with respect to a relevant asymmetric Finsler metric. We firstly transform the elastica problem into finding a geodesic in an orientation-lifted space. For this purpose, we denote by

 vθ=(cosθ,sinθ) (28)

the unit vector corresponding to .

Let be a regular curve with non-vanishing velocity and be its canonical orientation lifting. For any , is defined as being such that:

 vζ(t):=Γ′(t)∥Γ′(t)∥. (29)

According to the definition of in (28), one has

 (Γ′(t)∥Γ′(t)∥)⊥ =(vζ(t))⊥ =(−sin(ζ(t)),cos(ζ(t))), (30)

where denotes the vector that is perpendicular to a vector . It is known that

 ddt(Γ′(t)∥Γ′(t)∥)=κ(t)∥Γ′(t)∥(Γ′(t)∥Γ′(t)∥)⊥, (31)

where is the curvature of path . Then we have the following equations:

 ddtvζ(t) =ddt(cos(ζ(t)),sin(ζ(t))) =ζ′(t)(−sin(ζ(t)),cos(ζ(t))) =ζ′(t)(Γ′(t)∥Γ′(t)∥)⊥.

Thus, using (29) and (31), we have

 (32)

which yields to

 ζ′(t)=κ(t)∥Γ′(t)∥,∀t∈[0,1]. (33)

Using equations (27), (29) and (33), one has

 ℓ(Γ) =∫L0(1+ακ2(s))ds =∫10(1+α|ζ′(t)|2∥Γ′(t)∥2)∥Γ′(t)∥dt =∫10(∥Γ′(t)∥+α|ζ′(t)|2∥Γ′(t)∥)dt, (34)

where the Euclidean arc-length is defined as

 ds=∥Γ′(t)∥dt.

By the definition of , for any we have and

 ℓ(Γ)=∫10F∞(γ(t),γ′(t))dt, (35)

where we define the Finsler metric on the orientation-lifted domain by

 F∞(¯x,¯u):=⎧⎨⎩∥u∥+α|ν|2∥u∥, if u∝vθ,+∞,otherwise. (36)

for any orientation-lifted point , any vector in the tangent space, and where denotes positive collinearity. Note that any other lifting of would by construction of (36) have infinite energy, i.e., .

### 3.2 λ Penalized Asymmetric Finsler Elastica Metric Fλ

The Finsler metric  (36) is too singular to compute the global minimum of  (27) by directly applying the numerical algorithm such as the fast marching method (Mirebeau, 2014b). Hence we introduce a family of orientation-lifted Finsler elastica metrics over the orientation-lifted domain , depending on a penalization parameter as follows:

 Fλ(¯x,¯u):=√λ2∥u∥2+2αλ|ν|2−(λ−1)⟨vθ,u⟩, (37)

for any orientation-lifted point and any vector , and where is the unit vector associated to which denotes the position of in the orientation space .

As , the penalized Finsler elastica metric can be expressed as:

 Fλ(¯x,¯u)= √λ2∥u∥2+2αλ|ν|2−(λ−1)⟨vθ,u⟩ = λ∥u∥ ⎷1+α2|ν|2λ∥u∥2−(λ−1)⟨vθ,u⟩ = λ∥u∥(1+α|ν|2λ∥u∥2+O(λ−2))−(λ−1)⟨vθ,u⟩ = ∥u∥+α|ν|2∥u∥+(λ−1)(∥u∥−⟨vθ,u⟩) +O(λ−1) (38)

The term will vanish if vector is positively proportional to . Therefore, one has for any and any

 Fλ(¯x,¯u)→F∞(¯x,¯u),as λ→∞.

The metric (37) has precisely the required form formulated in (22), with tensor field as:

 MF(¯x)=diag(λ2,λ2,2αλ), (39)

and vector field

 →ωF(¯x)=(λ−1)(vθ,0), (40)

for any . Note that the definiteness constraint (23) is satisfied:

 ⟨→ωF(¯x),M−1F(¯x)→ωF(¯x)⟩=(1−λ−1)2<1,∀¯x∈¯Ω.

The anisotropy ratio characterizes the distortion between different orientations induced by a metric. Letting , and , the anisotropy ratio of the Finsler elastica metric (37) can be defined by:

 μ(Fλ):=sup¯x∈¯Ω{max∥¯w∥=∥¯v∥=1{Fλ¯x(¯w)Fλ¯x(¯v)}}, (41)

where the norm . As an example, for the Finsler elastica metric  (37) with and , we can show that in (41) is obtained by choosing and , so that .

Moreover, one can define the physical anisotropy ratio of the Finsler elastica metric by replacing by and the variables and in (41), respectively. In this case, for any , the physical anisotropy ratio is equal to and only depends on .

A crucial object for studying and visualizing the geometry distortion induced by a metric is Tissot’s indicatrix defined as the collection of unit balls in the tangent space. For point and , we define the unit balls for metrics and respectively by

 B∞¯x:={¯u=(u,ν)∈R2×R;F∞(¯x,¯u)≤1}, (42)

and

 Bλ¯x:={¯u=(u,ν)∈R2×R;Fλ(¯x,¯u)≤1}. (43)

Then any vector can be characterized by

 u⊥=0,u∥>0,and u∥+α|ν|2u∥≤1, (44)

where we introduce and as follows:

 u∥:=⟨u,vθ⟩,u⊥:=⟨u,v⊥θ⟩.

Using (44), one has

 (u∥−12)2+α|ν|2≤14. (45)

Thus is a flat 2D ellipse embedded in the 3D tangent space, and containing the origin on its boundary. Particularly, when , turns to a flat 2D disk of radius as shown in Fig. 2a.

On the other hand, when , a short computation shows that any vector is characterized by a quadratic equation

 λ2u2⊥+aλ(u∥−bλ2)2+α|ν|2≤cλ4, (46)

where are all . Hence is an ellipsoid, for instance see Fig. 2b with , almost flat in the direction of due to the large factor , which converges to the flat disk in the Haussdorf distance as .

Tissot’s indicatrix is also the control set in the optimal control interpretation of the Eikonal PDE (18). The Haussdorf convergence of the control sets guarantees that the minimal action map and minimal paths for the metric converge towards those of as . Elements of proof of convergence can be found in Appendix B.

### 3.3 Numerical Implementations

Numerically, anisotropy is related to the problem stiffness, hence to its difficulty. In Table 1, we show the computation time and the average number of Hopf-Lax updates required for each grid point by the adaptive stencils based fast marching method (Mirebeau, 2014b) for and different values of on a grid. This experiment was performed with a C++ implementation running on a standard 2.7GHz Intel I7 laptop with 16Gb RAM.

We observe on Table 1 a logarithmic dependence of computation time and average number of the Hopf-Lax updates per grid point with respect to anisotropy. These observations agree with the complexity analysis of the fast marching method presented in (Mirebeau, 2014b), yielding the upper bound , depending poly-logarithmically on the anisotropy ratio  (41), and quasi-linearly on the number of discretization points in the orientation-lifted domain . In contrast, numerical methods such as (Sethian and Vladimirsky, 2003) displaying a polynomial complexity in the anisotropy ratio would be unworkable. The iterative AGSI method (Bornemann and Rasch, 2006), on the other hand, requires hundreds of evaluations of the Hopf-Lax operator (25) per grid point to converge for large anisotropies, which also leads to prohibitive computation time, thus impractical. For or , the average numbers of the Hopf-Lax updates per grid required by the AGSI method are approximately and , respectively, while the numbers of Hopf-Lax from the fast marching method are only and , respectively, as demonstrated in Table 1.

In Fig. 3a, we show different Finsler elastica minimal paths, computed by the fast marching method (Mirebeau, 2014b), with (see (37)) and different values of . The arrows indicate the initial source and end points tangents. The cyan point denotes the initial position and the blue point indicates the end position. In Figs. 3b and 3c, we show the Finsler elastica minimal paths for different values of , with and respectively. In this experiment, the angle resolution is and the image size is . When , the metric is constant over the domain and degenerates to the isotropic orientation-lifted metric (11), since the coefficient in front of the term in (37) vanishes. Hence the minimal geodesics are straight lines, see Fig. 3a, that do not align with the prescribed endpoints tangents. From Fig. 3, one can point out that as and increasing, curvature penalization forces the extracted paths to gradually align with the prescribed endpoints tangents and take the elastica shape.

### 3.4 Image Data-Driven Finsler Elastica Metric P

We use in Section 3.1 for the sake of simplicity. In the general case, the metric  (36) and its approximation  (37) should be respectively replaced by and . Furthermore, in order to take into account the orientation information, we use an orientation-dependent speed function to replace . In this case, the data-driven Finsler elastica metric can be defined by

 P(¯x,¯u):=1Φ(¯x)Fλ(¯x,¯u),∀¯x∈¯Ω,∀¯u∈R3, (47)

and minimizing Eq. (26) is approximated for large values of by minimizing

 L(Γ) =∫101Φ(γ(t))Fλ(γ(t),γ′(t))dt =∫10P(γ(t),γ′(t))dt,

where is the orientation-lifted curve of . The data-driven Finsler elastica metric is asymmetric in the sense that for most vectors , one has

 P(¯x,¯u)≠P(¯x,−¯u),∀¯x∈¯Ω. (48)

This asymmetric property can help to build a closed contour passing through a collection of orientation-lifted points as discussed in Section 5.1.

The minimal action map associated to the data-driven Finsler elastica metric and an initial source point is the unique viscosity solution to the Eikonal PDE (18(Lions, 1982; Sethian and Vladimirsky, 2003):

 {P∗¯x(∇W¯s(¯x))=1,∀¯x∈¯Ω∖{¯s},W¯s(¯s)=0, (49)

where is the dual norm of defined by (19). The geodesic distance value between and depends on both the curvature and the orientation information of the minimal paths. When is sufficiently large, the spatial and angular resolutions are sufficiently small, the fixed point system (24) is properly solved, and the minimal paths are properly extracted by (21).

## 4 Computation of Data-Driven Speed Functions by Steerable Filters

In this section, we introduce two types of anisotropic speed functions for boundary detection and tubular structure extraction, both of which are based on the steerable filters.

### 4.1 Steerable Edge Detector

Jacob and Unser (2004) proposed a new class of edge detection filters based on the computational framework and the steerable property. Letting be a 2D isotropic Gaussian kernel with variance and , the computational steerable filter with order can be expressed as

 FMθ(x)=M∑τ=1τ∑ξ=0Kτ,ξ(θ)∂(τ−ξ)∂x(τ−ξ)∂ξ∂yξGσ(x), (50)

where and are the orientation-dependent coefficients which can be computed in terms of some optimality criteria. Particularly, when , the steerable filter becomes the classical Canny detector (Canny, 1986). For higher order steerable filters, for example, or , the orientation-dependent responses of the filters will be more robust to noise. Therefore, we set for the relevant experiments. Regarding the details of the computation of , we refer to (Jacob and Unser, 2004).

A color image is regarded as a vector valued map , for each . A multi-orientation response function of a color image can be computed by the steerable filter (50)

 h(x,θ)=133∑i=1|Ii(x)∗FMθ(x)|,∀θ∈[0,2π). (51)

For a gray level image , we have the simple computation of :

 h(x,θ)=|I(x)∗FMθ(x)|,∀θ∈[0,2π). (52)

The multi-orientation response function is symmetric in the sense that for any orientation , one has

 h(x,θπ)=h(x,θπ+π),∀x∈Ω.

### 4.2 Multi-Orientation Optimally Oriented Flux Filter

Optimally oriented flux filter is used to extract the local geometry of the image (Law and Chung, 2008). The oriented flux of an image , of dimension , is defined by the amount of the image gradient projected along the orientation flowing out from a 2D circle at point with radius :

 f(x;r,v)=∫∂Cr(∇(Gσ∗I)(x+rn)⋅v)(v⋅n)ds, (53)

where is a Gaussian with variance . is the outward normal of the boundary and is the infinitesimal length on . According to the divergence theory, one has

 f(x;r,v)=vTQ(x,r)v,

for some symmetric matrix :

 Q(x,r)=(∂xxGσ∂xyGσ∂yxGσ∂yyGσ)∗\mathbbm1r∗I(x), (54)

where is an indicator function of the circle .

Let be the eigenvalues of symmetric matrix and suppose that the intensities inside the tubular structures are darker than the background regions. If a point is inside the tubular structure, one has and , where is the optimal scale map

 r∗(x):=argmaxr{1rλ1(x,r)},∀x∈Ω. (55)

The scale normalized factor in (55) provides the responses of the optimally oriented flux filter a scale invariant property as discussed in (Law and Chung, 2008).

As shown in (Benmansour and Cohen, 2011), the optimally oriented flux filter is steerable, which means that we can construct a multi-orientation response function for any by:

 g(x,θ)=max{uTθQ(x,r∗(x))uθ,0},∀x∈Ω. (56)

where is a unit vector associated to .

In addition, the vesselness map , which indicates the probability of a pixel belonging to a vessel, can be calculated by:

 (57)

The vesselness map will be used to compute the isotropic Riemannian metric in the following section.

### 4.3 Computation of Data-Driven Speed Functions

A requirement on the speed function used by the data-driven Finsler elastica metric defined in (47) is that it should depend on the position and the orientation. Therefore, based on the orientation-dependent response function , defined in Section 4.1, the speed function that is used for object boundary detection can be computed by

 Φ(x,θ)=1+η(h(x,mod(θ+π/2,2π))∥h∥∞)p, (58)

for any and any orientation , where is defined as .

Similarly, we define the speed function for tubular structure extraction by using the function (56)

 Φ(x,θ)=1+η(g(x,mod(θ+π/2,2π))∥g∥∞)p, (59)

where and are positive constants. In this paper, we use for all the experiments.

## 5 Closed Contour Detection and Tubular Structure Extraction

We use the following convention in the remaining part of this paper: if is a point in , then we use to denote the orientation-lifted point that has the same physical position with but the opposite direction, where .