Da Chen
email: chenda@ceremade.dauphine.fr

Jean-Marie Mirebeau
email: jean-marie.mirebeau@math.u-psud.fr

Laurent D Cohen
email: cohen@ceremade.dauphine.fr

Global Minimum for a Finsler Elastica Minimal Path Approach

Keywords: Minimal Path, Geodesic, Eikonal Equation, Curvature Penalty, Euler Elastica Curve, Finsler Metric, Anisotropic Fast Marching Method, Image Segmentation, Perceptual Grouping, Tubular Structure Extraction

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

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 [START_REF] Kass | Snakes: Active contour models[END_REF]. A snake is a parametrized curve (locally) that minimizes the energy:

E( ) = Z 1 0 w 1 k 0 (t)k 2 + w 2 k 00 (t)k 2 + P (t) dt
with appropriate boundary conditions at the endpoints t = 0 and t = 1. 0 and 00 are the first-and second-order derivatives of the curve respectively. The positive constants w 1 and w 2 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 E 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 P , 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 [START_REF] Kass | Snakes: Active contour models[END_REF] 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 [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]. 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 [START_REF] Cohen | On active contour models and balloons[END_REF]. 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 [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Caselles | A geometric model for active contours in image processing[END_REF][START_REF] Malladi | Evolutionary fronts for topology-independent shape modeling and recovery[END_REF][START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Yezzi | A geometric snake model for segmentation of medical imagery[END_REF].

The difficulty of minimizing the non-convex snakes energy [START_REF] Kass | Snakes: Active contour models[END_REF] 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 [START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF][START_REF] Caselles | Geodesic active contours[END_REF]. To address this local minimum sensitivity issue, [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] proposed an Eikonal PDE-based minimal path model to find the global minimum of the geodesic active contours energy [START_REF] Caselles | Geodesic active contours[END_REF], 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 [START_REF] Sethian | Fast marching methods[END_REF], 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 [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] have important practical consequences, leading to various approaches for image analysis and medical imaging [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF][START_REF] Cohen | Multiple contour finding and perceptual grouping using minimal paths[END_REF].

In the basic formulations of the minimal paths-based interactive image segmentation models [START_REF] Appleton | Globally optimal geodesic active contours[END_REF][START_REF] Appia | Active geodesics: Region-based active contour segmentation with a global edge-based constraint[END_REF][START_REF] Mille | Combination of piecewise-geodesic paths for interactive segmentation[END_REF], the common proposal is that the object boundaries can be delineated by a set of minimal paths constrained by user-provided points. In [START_REF] Li | Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centrelines[END_REF][START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], 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, [START_REF] Benmansour | Fast object segmentation by growing minimal paths from a single point on 2D or 3D images[END_REF] proposed a growing minimal path model for object segmentation with unknown endpoints. This model can re-cursively 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. [START_REF] Kaul | Detecting curves with unknown endpoints and arbitrary topology using minimal paths[END_REF] applied the growing minimal path model to detect complicated curves with arbitrary topologies and developed criteria to stop the growth of the minimal paths. [START_REF] Rouchdy | Geodesic voting for the automatic extraction of tree structures. Methods and applications[END_REF] 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 [START_REF] Bougleux | Anisotropic Geodesics for Perceptual Grouping and Domain Meshing[END_REF][START_REF] Jbabdi | Accurate anisotropic fast marching for diffusion-based geodesic tractography[END_REF][START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF]. Such orientation enhancement can solve some shortcuts problems suffered by the isotropic minimal path models (Cohen and Kimmel, 1997;[START_REF] Li | Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centrelines[END_REF]. [START_REF] Kimmel | Optimal algorithm for shape from shading and path planning[END_REF] 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 [START_REF] Kimmel | Optimal algorithm for shape from shading and path planning[END_REF] was built over a higher dimensional domain by adding an extra orientation space to the image domain. [START_REF] Bekkers | A PDE Approach to Data-driven Sub-Riemannian Geodesics in SE (2)[END_REF] 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 [START_REF] Bekkers | A PDE Approach to Data-driven Sub-Riemannian Geodesics in SE (2)[END_REF] is based on a PDE approach with an upwind discretization and iterative evolution scheme, which requires expensive computation time. To solve this problem, [START_REF] Sanguinetti | Sub-Riemannian Fast Marching in SE (2)[END_REF] used the fast marching method (Mirebeau, 2014a) as the Eikonal solver to track the sub-Riemannian geodesics. The sub-Riemannian geodesic model [START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF] reintroduced curvature penalization to the geodesic energy, similar to the Euler elastica bending energy [START_REF] Nitzberg | The 2.1-D sketch[END_REF] 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. [START_REF] Schoenemann | A linear framework for region-based image segmentation and inpainting involving curvature penalization[END_REF] proposed a model to address the problems of using curvature regularization for regionbased 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 inpaint- ing results especially for objects with long and thin structures. [START_REF] Ulen | Shortest Paths with Higher-Order Regularization[END_REF] 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. [START_REF] Tai | A fast algorithm for Euler's elastica model using augmented Lagrangian method[END_REF] introduced an efficient method to solve the minimization problem of the Euler elastica energy with applications to image denosing, inpainting, and zooming. Other approaches of interest using the curvature penalization include the image segmentation models such as [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF][START_REF] Schoenemann | The Elastic Ratio: Introducing Curvature Into Ratio-Based Image Segmentation[END_REF][START_REF] Zhu | Image segmentation using Euler's elastica as the regularization[END_REF] and the image inpainting model introduced by [START_REF] Shen | Euler's elastica and curvature-based inpainting[END_REF].

Motivation

In contrast with the classical snakes energy [START_REF] Kass | Snakes: Active contour models[END_REF], 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 and1c are the minimal paths obtained by using the isotropic Riemannian metric [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] and the anisotropic Riemannian metric [START_REF] Bougleux | Anisotropic Geodesics for Perceptual Grouping and Domain Meshing[END_REF], 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 [START_REF] Nitzberg | The 2.1-D sketch[END_REF][START_REF] Mumford | Elastica and computer vision[END_REF]. 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.

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 the corresponding 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 closed contour detection, perceptual grouping and tubular structure extraction. Closed contour detection is performed in an interactive manner, where the contour is concatenated by a set of piecewise smooth geodesics. It connects a set of user-provided orientation lifted points. With a procedure similar to the closed contour detection method, we apply our model to perceptual grouping based on the criteria of connectivity and smoothness. Moreover, we also provide a method for simplifying the tubular structure extraction application.

More contributions have been added regarding the original conference paper presented in [START_REF] Chen | Global Minimum for Curvature Penalized Minimal Path Method[END_REF], such as the applications of interactive closed contour detection for image segmentation and perceptual grouping, respectively.

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

Cohen-Kimmel Minimal Path Model

Let ⌦ ⇢ R 2 or R 3 denote the image domain. The classical Cohen-Kimmel minimal path model [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] was designed to find the global minimum of the geodesic energy which measures the length of Lipschitz paths : [0, 1] ! ⌦ as follows

L I ( ) = Z 1 0 ⇣ w + P (t) ⌘ k 0 (t)k dt, (1)
where k • k denotes the canonical Euclidean norm. w is a positive constant and P : ⌦ ! R + is a potential function which usually depends on the image gradient magnitude or gray levels as suggested by [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF].

The minimal action map U : ⌦ ! R + associated to a source point s, is the global minimum of the length (1) among all paths joining s to any point x 2 ⌦:

U (x) := inf{L I ( ); 2 A s,x }, 8x 2 ⌦, (2)
where A s,x is the collection of all Lipschitz paths. The Cohen-Kimmel minimal path model considers an orientation-independent isotropic Riemannian metric. For any vector u 2 R 2 or R 3 , this metric is defined by

F I (x, u) := p hu, M I (x) ui, 8x 2 ⌦, (3)
where h•, •i denotes the scalar product over R 2 or R 3 and M I is a symmetric positive definite tensor field which is proportional to the identity matrix

I d : M I (x) = w + P (x) 2 I d , 8x 2 ⌦. (4)
The minimal action map U defined in (2) satisfies the following Eikonal PDE [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]:

( krU (x)k = w + P (x), 8x 2 ⌦\{s}, U (s) = 0, ( 5
)
where rU is the Euclidean gradient of U with respect to the position in the domain. The geodesic C s,x , linking the source point s to x, is obtained through the gradient descent ordinary differential equation (ODE), solved backwards from t = 1 to t = 0:

( C 0 s,x (t) / rU(C s,x (t)), C s,x (1) = x, (6)
where / denotes positive collinearity.

Anisotropic Riemannian Metric Extension

Bougleux et al ( 2008) and [START_REF] Jbabdi | Accurate anisotropic fast marching for diffusion-based geodesic tractography[END_REF] extended the isotropic Riemannian metric (3) to the anisotropic Riemannian case invoking a symmetric positive definite tensor field

M A (x) = n X i=1 P i (x) v i (x) v T i (x), 8 x 2 ⌦, ( 7
)
where n is the dimension of the domain ⌦ and the vector fields v i are related to the image data. P i : ⌦ ! R + are the image data-driven functions associated to v i . Based on the tensor field M A (7), the anisotropic Riemannian metric F A can be expressed as:

F A (x, u) := p hu, M A (x) ui. (8)
The minimal action map U associated to the anisotropic Riemannian metric F A satisfies the Eikonal PDE:

( krU (x)k M 1 A (x) = 1, 8 x 2 ⌦\{s}, U (s) = 0, (9)
where we denote that kuk M = p hu, M ui. The geodesic C s,x can be tracked by solving the ODE, again solved backwards from t = 1 to t = 0:

( C 0 s,x (t) / M 1 A (C s,x (t))rU (C s,x (t)), C s,x (1) = x.
(10) [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] proposed an anisotropic radius lifted Riemannian metic, as a special case of F A (8), to deal with the problem of extracting both the centerline positions and the corresponding radius values of a tubular structure simultaneously.

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 ⌦ ⇢ R 2 can be extended to the orientation lifted space ⌦ by product with an abstract orientation space S 1 [START_REF] Kimmel | Optimal algorithm for shape from shading and path planning[END_REF], i.e., ⌦ = ⌦ ⇥ S 1 ⇢ R 3 and the problem is to find a minimal path in the new lifted space ⌦. Each point x in the orientation lifted path is thus a pair composed of a point x in the image domain ⌦ and an orientation ✓, i.e., x = (x, ✓).

For any point x and any vector ū = (u, ⌫) 2 R 2 ⇥ R, where ⌫ denotes the orientation variation, the isotropic orientation lifted Riemannian metric F IL is expressed by

F IL (x, ū) := 1 IL (x) p kuk 2 + ⇢|⌫| 2 , ( 11
)
where IL : ⌦ ! R + is an orientation-dependent speed function relying on the image data and ⇢ > 0 is a constant. Following the form of (3), the metric F IL defined in (11) can be written by using a tensor field

M IL M IL (x) = 1 2 IL (x) D ⇢ , (12)
where D ⇢ = diag(1, 1, ⇢) 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 [START_REF] Kimmel | Optimal algorithm for shape from shading and path planning[END_REF][START_REF] Péchaud | Extraction of tubular structures over an orientation domain[END_REF].

The curve length of an orientation lifted curve := ( , ✓) : [0, 1] ! ⌦ with respect to the isotropic orientation lifted Riemannian metric F IL can be expressed as

L IL ( ) = Z 1 0 F IL ( (t), 0 (t)) dt = Z 1 0 1 IL ( (t)) p h 0 (t), M IL ( (t)) 0 (t)i dt = Z 1 0 1 IL ( (t)) p k 0 (t)k 2 + ⇢|✓ 0 (t)| 2 dt, (13)
where 0 (t) = ( 0 (t), ✓ 0 (t)) 2 R 3 , for all t 2 [0, 1]. The idea of orientation lifting was applied to tubular structure extraction by [START_REF] Péchaud | Extraction of tubular structures over an orientation domain[END_REF] with additional radius lifting, typically accounting for the radius r of the tubular structure in the processed image. Orientation lifting often improves the results from [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF][START_REF] Li | Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centrelines[END_REF], but suffers from the fact that for the geodesic curve length of an orientation lifted curve = ( , ✓) : [0, 1] ! ⌦, nothing in (13) constrains the path direction 0 (t) to align with the orientation vector associated to orientation ✓(t), where t 2 [0, 1]. In other words, this isotropic orientation lifted Riemannian metric F IL cannot penalize the curvature of the geodesic, a point which is addressed in this paper.

General Minimal Path Model and Finsler Metric

The general minimal path problem is posed on a bounded domain ⌦ ⇢ R n equipped with a metric F(x, u) depending on positions x 2 ⌦ and orientations u 2 R n . This metric F defines at each point x a norm

F x (u) := F(x, u). (14)
These norms must be positive F x (u) > 0 whenever u 6 = 0, 1-homogeneous, and obey the triangular inequality. However, in general we allow them to be asymmetric :

F x (u) 6 = F x ( u), (15)
One can measure the curve length of a regular curve with respect to the metric F:

L F ( ) = Z 1 0 F (t), 0 (t) dt, (16)
The minimal action map U (x) is defined by:

U (x) := inf{L F ( ); 2 A s,x }, (17)
which is the unique viscosity solution to an Eikonal PDE [START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF][START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms[END_REF]:

( F ⇤ x rU (x) = 1, 8x 2 ⌦\{s}, U (s) = 0, (18)
where F ⇤ x is the dual norm of F

x defined for all u 2 R n by

F ⇤ x (u) := sup v6 =0 hu, vi F x (v) . (19)
Based on the definition of the dual norm in (19), the corresponding optimal direction map is then obtained by

(x, u) := arg max v6 =0 hu, vi F x (v) , 8 x 2 ⌦, 8u 2 R n . ( 20
)
The geodesic C s,x , joining the source point s to any point x 2 ⌦, is the solution to the following ODE involving the minimal action map U and the optimal direction map

( C 0 s,x (t) / ⇣ C s,x (t), rU C s,x (t) ⌘ , C s,x (1) = x. (21)
Numerically, the ODE expressed in ( 21) is solved backwards from t = 1 to t = 0, using Heuns or Runge-Kutta's methods, or more robustly using the numerical method proposed by Mirebeau (2014a).

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

F(x, u) := p hu, M(x) ui h!(x), ui, (22)
for all x 2 ⌦ and any vector u 2 R n . The asymmetric part should obey the following smallness condition to ensure that the Finsler metric F is positive:

h!(x), M 1 (x) !(x)i < 1, 8x 2 ⌦. (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 M is proportional to the identity matrix or diagonal matrix D ⇢ in [START_REF] Cohen | Multiple contour finding and perceptual grouping using minimal paths[END_REF]. Therefore, the general Eikonal PDE (18) and the geodesic back tracking equation (21) reduce respectively to equations ( 5) and ( 6) for the isotropic case, and respectively to equations ( 9) and (10) for the anisotropic case.

Computation of the Minimal Action Map

In order to estimate the minimal action map U , presented in ( 17) and ( 18), a discretization grid Z of the image domain ⌦ is introduced -or of the extended domain ⌦ in the case of an orientation lifted metric. For each point x 2 Z, a small mesh S(x) of a neighbourhood of x with vertices in Z is constructed. For example, S(x) can be the square formed by the four neighbours of x in the classical fast marching method [START_REF] Sethian | Fast marching methods[END_REF] on a regular orthogonal grid used in [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]. 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 U , with the initial source point s, is calculated by solving the fixed point system [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF]:

( U (x) = ⇤U (x), for all x 2 Z\{s}, U (s) = 0, (24)
where the involved Hopf-Lax update operator is defined by:

⇤U (x) := min y2@S(x) n F(x, x y) + I S(x) U (y) o , (25)
where I S(x) denotes the piecewise linear interpolation operator on the mesh S(x), and y lies on the boundary of S(x). The equality U (x) = ⇤U (x), replacing in (24) the Eikonal PDE: F

x (rU (x)) = 1 of ( 18), is a discretization of Bellman's optimality principle, which is similar in spirit to the Tsitsiklis approach [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF]. It reflects the fact that the minimal geodesic C s,x , from s to x, has to cross the mesh boundary @S(x) at least once at some point y; therefore it is the concatenation of a geodesic C s,y from s to y, which length is approximated by piecewise linear interpolation, and a very short geodesic C y,x from y to x, approximated by a segment of geodesic curve length F(x, x y). The N -dimensional fixed point system (24), with N = #Z 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 [START_REF] Sethian | Fast marching methods[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] using the square formed neighbourhood S 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 S 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 S(x) at each point x 2 ⌦ or ⌦ satisfies some geometric acuteness property depending on the local metric F(x, •). 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 [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF]. 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).

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 ) denotes the angle space with periodic boundary conditions and ⌦ ⇢ R 2 denotes the image domain.

⌦ = ⌦ ⇥ S 1 ⇢ R 3 , where S 1 = [0, 2⇡

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 [START_REF] Nitzberg | The 2.1-D sketch[END_REF] and [START_REF] Mumford | Elastica and computer vision[END_REF]. They minimized the following bending energy:

L( ˆ ) := Z L 0 1 0 ˆ (s) 1 + ↵ 2 (s) ds, (26)
where ˆ : [0, L] ! ⌦ is a regular curve with non-vanishing velocity,  is the curvature of curve ˆ , L is the Euclidean curve length of ˆ , and s is the arc-length. Parameter ↵ > 0 is a constant. 0 is an image data-driven speed function.

For the sake of simplicity, we set 0 ⌘ 1, yielding the simplified Euler elastica bending energy 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

`( ˆ ) = Z L 0 1 + ↵ 2 (s) ds, (27)
v ✓ = (cos ✓, sin ✓) (28)
the unit vector corresponding to ✓ 2 S 1 . Let : [0, 1] ! ⌦ be a regular curve with non-vanishing velocity and = ( , ✓) : [0, 1] ! ⌦ be its canonical orientation lifting. For any t 2 [0, 1], ✓(t) is defined as being such that:

v ✓(t) := 0 (t) k 0 (t)k . ( 29
)
According to the definition of v ✓ in ( 28), one has

✓ 0 (t) k 0 (t)k ◆ ? = (v ✓(t) ) ? = ( sin ✓(t), cos ✓(t)), (30)
where u ? denotes the vector that is perpendicular to a vector

u. It is known that d dt ✓ 0 (t) k 0 (t)k ◆ = (t)k 0 (t)k ✓ 0 (t) k 0 (t)k ◆ ? , ( 31
)
where  is the curvature of path . Then we have the following equations:

d dt v ✓(t) = d dt (cos ✓(t), sin ✓(t)) = ✓ 0 (t) sin ✓(t), cos ✓(t) = ✓ 0 (t) ✓ 0 (t) k 0 (t)k ◆ ? .
Thus, using ( 29) and ( 31), we have

✓ 0 (t) ✓ 0 (t) k 0 (t)k ◆ ? = (t)k 0 (t)k ✓ 0 (t) k 0 (t)k ◆ ? , (32)
which yields to

✓ 0 (t) = (t)k 0 (t)k, 8t 2 [0, 1]. (33)
Using equations ( 27), ( 29) and ( 33), one has

`( ) = Z L 0 1 + ↵ 2 (s) ds = Z 1 0 ✓ 1 + ↵ |✓ 0 (t)| 2 k 0 (t)k 2 ◆ k 0 (t)k dt = Z 1 0 ✓ k 0 (t)k + ↵ |✓ 0 (t)| 2 k 0 (t)k ◆ dt, (34)
where the Euclidean arc-length is defined as

ds = k 0 (t)kdt. By the definition of , for any t 2 [0, 1] we have 0 (t) = ( 0 (t), ✓ 0 (t)) and `( ) = Z 1 0 F 1 (t), 0 (t) dt, (35)
where we define the Finsler metric F 1 on the orientation lifted domain ⌦ by

F 1 (x, ū) := ( kuk + ↵ |⌫| 2 kuk , if u / v ✓ , +1, otherwise. ( 36
)
for any orientation lifted point x = (x, ✓) 2 ⌦, any vector ū = (u, ⌫) 2 R 2 ⇥ R in the tangent space, and where / denotes positive collinearity. Note that any other lifting

˜ (t) = ( (t), ✓(t)) of (t)
would by construction of (36) have infinite energy, i.e., `(˜ ) = 1.

Penalized Asymmetric Finsler Elastica Metric F

The Finsler metric F 1 (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 1 as follows:

F (x, ū) := p 2 kuk 2 + 2↵ |⌫| 2
( 1)hv ✓ , ui, (37) for any orientation lifted point x = (x, ✓) 2 ⌦ and any vector ū = (u, ⌫) 2 R 2 ⇥ R, and where v ✓ = (cos ✓, sin ✓) is the unit vector associated to ✓ which denotes the position of x in the orientation space S 1 .

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

F (x, ū) = p 2 kuk 2 + 2↵ |⌫| 2 ( 1)hv ✓ , ui = kuk s 1 + ↵ 2|⌫| 2 kuk 2 ( 1)hv ✓ , ui = kuk(1 + ↵|⌫| 2 kuk 2 + O( 2 )) ( 1)hv ✓ , ui =kuk + ↵|⌫| 2 kuk + ( 1)(kuk hv ✓ , ui) + O( 1 ) (38)
The term kuk hv ✓ , ui will vanish if vector u is positively proportional to v ✓ . Therefore, one has for any x and any u F (x, ū) ! F 1 (x, ū), as ! 1. The metric F (37) has precisely the required form formulated in ( 22), with tensor field M := M F as:

M F (x) = diag( 2 , 2 , 2↵ ), (39)
and vector field

! := ! F ! F (x) = ( 1)(v ✓ , 0), (40)
for any x = (x, ✓) 2 ⌦. Note that the definiteness constraint ( 23) is satisfied:

h! F (x), M 1 F (x) ! F (x)i = (1 1 ) 2 < 1, 8 x 2 ⌦.
The anisotropy ratio characterizes the distortion between different orientations induced by a metric. Letting x = (x, ✓),

w = (w, ⌫ w ) 2 R 2 ⇥ R and v = (v, ⌫ v ) 2 R 2 ⇥ R
, the anisotropy ratio µ(F ) of the Finsler elastica metric F (37) can be defined by:

µ(F ) := sup x2 ⌦ ⇢ max k wk=kvk=1 ⇢ F x ( w) F x (v) , (41)
where the norm F x (•) = F (x, •). As an example, for the Finsler elastica metric F (37) with 2 and ↵ = 1, we can show that µ(F ) in ( 41) is obtained by choosing w = ( v ✓ , 0) and v = (v ✓ , 0), so that µ(F ) = 2 1. Moreover, one can define the physical anisotropy ratio of the Finsler elastica metric F by replacing by ws = (w, 0) and vs = (v, 0) the variables w and v in (41), respectively. In this case, for any ↵, the physical anisotropy ratio is equal to 2 1 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 x = (x, ✓) 2 ⌦ and 2 [1, 1), we define the unit balls for metrics F 1 and F respectively by

B 1 x := {ū = (u, ⌫) 2 R 2 ⇥ R; F 1 (x, ū)  1}, (42)
and

B x := {ū = (u, ⌫) 2 R 2 ⇥ R; F (x, ū)  1}. ( 43
)
Then any vector ū = (u, ⌫) 2 B 1 x can be characterized by

u ? = 0, u k > 0, and u k + ↵ |⌫| 2 u k  1, (44)
where we introduce u k and u ? as follows:

u k := hu, v ✓ i, u ? := hu, v ? ✓ i. Using (44), one has ✓ u k 1 2 ◆ 2 + ↵ |⌫| 2  1 4 . ( 45
)
Thus B 1

x is a flat 2D ellipse embedded in the 3D tangent space, and containing the origin on its boundary. Particularly, when ↵ = 1, B 1

x turns to a flat 2D disk of radius 1/2 as shown in Fig. 2a.

On the other hand, when < 1, a short computation shows that any vector ū = (u, ⌫) 2 B x is characterized by a quadratic equation

2 u 2 ? + a ✓ u k b 2 ◆ 2 + ↵ |⌫| 2  c 4 , ( 46
)
where a , b , c are all 1 + O(1/ ). Hence B x is an ellipsoid, for instance see Fig. 2b with ↵ = 1, almost flat in the direction of v ? ✓ due to the large factor /2, which converges to the flat disk B 1

x in the Haussdorf distance as ! 1. 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 F converge towards those of F 1 as ! 1. Elements of proof of convergence can be found in Appendix B.

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 ↵ = 500 and different values of on a 300 2 ⇥ 108 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 O(N (ln µ) 3 +N ln N ), depending poly-logarithmically on the anisotropy ratio µ (41), and quasi-linearly on the number N of discretization points in the orientation lifted domain ⌦. In contrast, numerical methods such as [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms[END_REF] displaying a polynomial complexity O(µ 2 N ln N ) in the anisotropy ratio would be unworkable. The iterative AGSI method [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], 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 = 30 or 100, the average numbers of the Hopf-Lax updates per grid required by the AGSI method are approximately 86 and 182, respectively, while the numbers of Hopf-Lax from the fast marching method are only 6.49 and 7.27, 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 ↵ = 500 (see ( 37 blue point indicates the end position. In Figs. 3b and3c, we show the Finsler elastica minimal paths for different values of ↵, with = 100 and = 300 respectively. In this experiment, the angle resolution is ✓ s = 2⇡/108 and the image size is 300 ⇥ 300. When = 1, the metric F is constant over the domain ⌦ and degenerates to the isotropic orientation lifted metric F I (11), since the coefficient in front of the term hv ✓ , ui 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.

Image Data-Driven Finsler Elastica Metric P

We use 0 ⌘ 1 in Section 3.1 for the sake of simplicity. In the general case, the metric F 1 (36) and its approximation F (37) should be respectively replaced by 1 0 F 1 and 1 0 F . Furthermore, in order to take into account the orientation information, we use an orientation dependent speed function : ⌦ ! R + to replace 0 . In this case, the datadriven Finsler elastica metric can be defined by

P(x, ū) := 1 (x) F (x, ū), 8x 2 ⌦, 8ū 2 R 3 , (47)
and equation ( 26) becomes

L( ) = Z 1 0 1 ( (t)) F ( (t), 0 (t))dt = Z 1 0 P( (t), 0 (t))dt,
where is the orientation lifted curve of . The data-driven Finsler elastica metric P is asymmetric in the sense that for most vectors ū 6 = 0, one has

P(x, ū) 6 = P(x, ū), 8 x 2 ⌦. ( 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 W s associated to the data-driven Finsler elastica metric P and an initial source point s is the unique viscosity solution to the Eikonal PDE (18) [START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF][START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations: Theory and algorithms[END_REF]:

( P ⇤ x rW s(x) = 1, 8x 2 ⌦\{s}, W s(s) = 0, (49)
where P ⇤ x is the dual norm of P x(•) := P(x, •) defined by ( 19). The geodesic distance value W s(x) between x and s 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).

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.

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 G be a 2D isotropic Gaussian kernel with variance and x = (x, y), the computational steerable filter F M ✓ with order M (Jacob and Unser, 2004) can be expressed as:

F M ✓ (x) = M X ⌧ =1 ⌧ X ⇠=0 K ⌧,⇠ (✓) @ (⌧ ⇠) @x (⌧ ⇠) @ ⇠ @y ⇠ G (x), (50)
where ✓ 2 [0, 2⇡) and K ⌧,⇠ are the orientation-dependent coefficients which can be computed in terms of some optimality criteria. Particularly, when M = 1, the steerable filter F 1 ✓ becomes the classical Canny detector [START_REF] Canny | A computational approach to edge detection[END_REF].

For higher order steerable filters, for example, M = 3 or M = 5, the orientation-dependent responses of the filters will be more robust to noise. Therefore, we set M = 5 for the relevant experiments. Regarding the details of the computation of K ⌧,⇠ , we refer to [START_REF] Jacob | Design of steerable filters for feature detection using canny-like criteria[END_REF]. A color image is regarded as a vector valued map I :

⌦ ! R 3 , I(x) = [I 1 (x), I 2 (x), I 3 (x)] for each x 2 ⌦. A multi-orientation response function h : ⌦ ! R + of a color image I can be computed by the steerable filter F M ✓ (50) h(x, ✓) = 1 3 3 X i=1 |I i (x) ⇤ F M ✓ (x)|, 8 ✓ 2 [0, 2⇡). ( 51
)
For a gray level image I : ⌦ ! R, we have the simple computation of h:

h(x, ✓) = |I(x) ⇤ F M ✓ (x)|, 8 ✓ 2 [0, 2⇡). ( 52
)
The multi-orientation response function h is symmetric in the sense that for any orientation ✓ ⇡ 2 [0, ⇡), one has

h(x, ✓ ⇡ ) = h(x, ✓ ⇡ + ⇡), 8x 2 ⌦.

Multi-Orientation Optimally Oriented Flux Filter

Optimally oriented flux filter is used to extract the local geometry of the image. The oriented flux [START_REF] Law | Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux[END_REF] of an image I : ⌦ ! R + , of dimension 2, is defined by the amount of the image gradient projected along the orientation v flowing out from a 2D circle at point x = (x, y) 2 ⌦ with radius r:

f (x; r, v) = Z @Cr (r(G ⇤ I)(x + rn) • v)(v • n) ds, ( 53
)
where G is a Gaussian with variance and n is the outward unit normal vector along @C r . ds is the infinitesimal length on the boundary of C r . According to the divergence theory, one has

f (x; r, v) = v T Q(x, r) v, for some symmetric matrix Q(x, r): Q(x, r) = ✓ @ xx G @ xy G @ yx G @ yy G ◆ ⇤ r ⇤ I(x), ( 54
)
where r is an indicator function of the circle C r . Let 1 (•)

2 (•) be the eigenvalues of symmetric matrix Q(•) and suppose that the intensities inside the tubular structures are darker than the background regions. If a point

x is inside the tubular structure, one has 1 (x, r ⇤ (x)) 0 and 2 (x, r ⇤ (x)) ⇡ 0, where r ⇤ is the optimal scale map defined by:

r ⇤ (x) = arg max r ⇢ 1 r 1 (x, r) , 8x 2 ⌦. ( 55
)
The scale normalized factor 1/r in (55) provides the responses of the optimally oriented flux filter a scale invariant property [START_REF] Law | Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux[END_REF].

As shown in [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], the optimally oriented flux filter is steerable, which means that we can construct the multi-orientation response function g :

⌦ ! R + for any ✓ 2 [0, 2⇡) by: g(x, ✓) = max u T ✓ Q(x, r ⇤ (x)) u ✓ , 0 , 8 x 2 ⌦. ( 56
)
where u ✓ is a unit vector associated to ✓. In addition, the vesselness map V n : ⌦ ! R, which indicates the probability of a pixel x belonging to a vessel, can be calculated by:

V n (x) = max ⇢ max r ⇢ 1 r 1 (x, r) , 0 . (57)
The vesselness map V n will be used to compute the isotropic Riemannian metric in the following section.

Computation of Data-Driven Speed Functions

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

(x, ✓) = 1 + ⌘ ✓ h(x, ✓) khk 1 ◆ p , (58)
for any x 2 ⌦ and any orientation ✓ 2 [0, 2⇡).

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

(x, ✓) = 1 + ⌘ ✓ g(x, ✓) kgk 1 ◆ p , (59)
where ⌘ and p are positive constants. In this paper, we use p = 2 for all the experiments.

Closed Contour Detection and Tubular Structure Extraction

We use the following convention in the remaining part of this paper: if p = (p, ✓) is a point in ⌦, then we use p † = (p, mod (✓ + ⇡), 2⇡) to denote the orientation lifted point that has the same physical position p with p but the opposite direction, where ✓ 2 [0, 2⇡).

Closed Contour Detection as a Set of Piecewise Finsler Elastica Minimal Paths

In this section, we propose an interactive closed contour detection method based on the Finsler elastica minimal paths constrained by a set H of m user-provided physical points H := {x i 2 ⌦, i = 1, 2, ..., m; m 2}, all of which are on the target object boundary. As discussed in Section 3, the Finsler elastica metric P is defined on the orientation lifted space ⌦. Thus, we build an orientation lifted collection D of H by

D := n xi = (x i , ✓ i ), x † i = x i , mod(✓ i + ⇡, 2⇡) ; i = 1, 2, ..., m, x i 2 H, and ✓ i 2 [0, 2⇡) o , (60)
where the directions ✓ i are manually specified in this paper. Corresponding to each physical point x i 2 H, there exist two orientation lifted points: xi , x † i 2 D, which have opposite tangents. We show an example of these orientation lifted points in Fig. 4a, where the physical positions and the corresponding tangents are denoted by dots and arrows, respectively.

The basic idea of the proposed closed contour detection method is to identify a set of ordered vertices from the collection D, and to join these detected vertices by a set of piecewise Finsler elastica minimal paths.

We start the closed contour detection procedure by selecting a physical position from H, say x 1 . The corresponding orientation lifted points of x 1 are denoted by x1 , x † 1 2 D. Once x 1 is specified, we remove both x1 and x † 1 from D. As shown in Fig. 4a, x1 andx † 1 are denoted by a red dot and two arrows with opposite directions.

Let ā⇤ 2 D be the closest orientation lifted point to x1 in terms of the geodesic distance W x1 (49) with respect to the data-driven Finsler elastica metric P, i.e.

ā⇤ = arg min

z2D W x1 (z). (61)
Similarly to ā⇤ , the closest orientation lifted point c⇤ 2 D of x † 1 can be detected. By the detected points ā⇤ and c⇤ , the first pair of successive vertices (q 1 , q2 ) is determined simultaneously using the following criterion:

(q 1 , q2 ) = ( (x 1 , ā⇤ ), if W x1 (ā ⇤ ) < W x † 1 (c ⇤ ), (x † 1 , c⇤ ), otherwise. (62)
In Fig. 4b, we illustrate the vertices q1 and q2 by the red and green dots with arrows, respectively. If the minimal action map W x1 (resp. W x †

1

) is computed via the fast marching method (Mirebeau, 2014b), ā⇤ (resp. c⇤ ) is the first vertex reached by the fast marching front, which is monotonically advancing. Once the first pair of successive vertices (q 1 , q2 ) is found, the geodesic C q1, q2 (red curve in Fig. 4b) can be tracked by ( 21), and both q2 , q † 2 are removed from D. If the number of physical points is 2, i.e. m = 2, the closed contour detection procedure can be terminated. The geodesic C q2, q1 is tracked by taking q2 as the initial source point and q1 as the end point.

(a) q 1 q 2 (b) q 1 q 2 q 3 (c) q 1 q 2 q 3 q 4 (d) q 1 q 2 q 3 q 4 (e) (f)
If m > 2, the subsequent vertex qi with i 3 is identified from the remaining points of D by searching for the nearest neighbour of the vertex qi 1 in terms of geodesic distance, i.e.

qi = arg min

z2D W qi 1 (z). ( 63
)
After the detection of the vertex qi , we remove both qi and q † i from D. Again the geodesic C qi 1 ,qi can be tracked, as denoted by the green curve in Fig. 4c.

The procedure of detecting the nearest neighbor from the set of remaining orientation lifted points is recursively carried out according to the criterion (63) until m ordered vertices have been identified. Then the geodesic C qm, q1 , which is denoted by the cyan curve in Fig. 4e, is computed by simply allowing qm to be the initial source point and q1 to be the end point. The final closed contour, denoted by C, is defined as the concatenation of all of the detected Finsler elastica minimal paths as demonstrated in Fig. 4f.

In summary, the proposed closed contour detection procedure aims to seeking a set of m pairs of successive ori-entation lifted points from D:

= m 1 [ i=1 (q i , qi+1 ) [ (q m , q1 ) , (64)
and a closed contour C contains a set of Finsler elastica minimal paths, joining all the pairs of vertices in . This method simply matches orientation lifted points by pairs, joining a vertex to the remaining nearest neighbor with respect to the curvature-penalized geodesic distance, so as to form a closed contour located at the expected object boundaries. Note importantly, that the obtained piecewise geodesic contour is smooth (C 1 differentiable) since the initial source and end orientation lifted points of the consecutive geodesics have both matching positions q i and orientations ✓ i . In fact, we find a closed contour passing through all the orientation lifted points in a greedy manner. Instead of trying out all possible combinations of Finsler elastica minimal paths, we use a greedy searching strategy that is performed with a low complexity. The problem we solve here is similar to the NPhard traveling salesman problem, where the cities are represented by the orientation lifted points qi 2 D.

Perceptual Grouping via the Curvature-Penalized Geodesic Distance

Perceptual grouping is relevant to the task of curve reconstruction and completion [START_REF] Cohen | Multiple contour finding and perceptual grouping using minimal paths[END_REF]. The geodesic distance based perceptual grouping model was firstly introduced by Cohen (2001) using the concept of saddle point. The basic idea of this model is to identify each pair of points which has to be linked by a minimal path from a set of key points. Later on, [START_REF] Bougleux | Anisotropic Geodesics for Perceptual Grouping and Domain Meshing[END_REF] improved this grouping model by using path orientations and structure tensors. However, neither of the mentioned grouping methods considered the curvature penalization.

In this section, we address the perceptual grouping problem of finding n closed contours, each of which is formed by a set of piecewise Finsler elastica minimal paths. Each closed contour C

i passes through all the ordered orientation lifted points involved in D i ✓ D, where D is defined in (60) and i = 1, 2, 3, ..., n.

We initialize the proposed perceptual grouping procedure by specifying a physical position x 1 , where the orientation lifted points of x 1 , denoted by x1 and x † 1 , are involved in D. Both x1 and x † 1 are removed from D after detecting the respective nearest vertices that correspond to x1 and x⇤ 1 by ( 61). As a consequence, the first two vertices q1 , q2 are identified using the criterion of ( 62), and the geodesic C q1, q2 is recovered by solving the ODE in (21). Once the vertices q1 and q2 are detected, we add q1 , q2 to D 1 , remove q2 , q † 2 from D and compensate q1 to D.

Similar to the closed contour detection method, the next vertex qi with i 3 is found based on the criterion of ( 63) and the detected vertex qi 1 . Following the detection of vertex qi , we add qi to D 1 , remove qi , q † i from D and track the geodesic C qi 1 ,qi that joins the vertices qi 1 and qi . This perceptual grouping procedure is carried out by recursively searching for new vertices. Once the vertex q1 is detected again according to the criterion (63), we stop the construction of D 1 after removing q1 from D, and back track the geodesic from q1 . The desired closed contour C 1 can be obtained by concatenating all the detected Finsler elastica minimal paths with initial source points and end points in D 1 .

We start to build the collection D 2 by choosing a new physical point as initialization. This initial physical point is obtained from the remaining orientation lifted points of D. Similar to the procedure of constructing D 1 , we build the collection D 2 from the remaining orientation lifted points of D. The procedure of building the collections D i can be terminated when n such collections have been identified or when the collection D is empty. One can note that the constructed collections D i follow

D i \ D j = ?, 8i 6 = j.
In contrast to the closed contour detection method described in Section 5.1, we do not enforce all of the orientation lifted points in D to be used in the perceptual grouping procedure.

Tubular Structure Extraction

In this section, we apply the proposed Finsler elastica minimal path model to the tubular structure extraction, where the centerlines of the tubular structures are represented by the Finsler elastica minimal paths.

The minimal paths with the proposed data-driven Finsler elastica metric depend on the tangents of both the initial source point and the end point. To simplify the initialization procedure, we firstly compute the optimal orientation map, denoted by ⇥ : ⌦ ! [0, ⇡), which minimizes the multiorientation response function g in ( 56):

⇥(x) = arg min ✓2[0,⇡) {g(x, ✓)}, 8x 2 ⌦. ( 65
)
Once the optimal orientation map ⇥ is obtained, for the initial position s 2 ⌦, one can obtain two orientation lifted points s = (s, ⇥(s)) and s † . Additionally, for any end position p i 2 ⌦ (i = 1, 2, • • • , n), the corresponding orientation lifted end points are defined by pi = (p i , ⇥(p i )) and p † i . For each set of orientation lifted end points {p i , p † i }, we can extract four possible geodesics, each of which joins an initial source point in {s, s † } to an end point in {p i , p † i }. The goal in this section is to search for a geodesic C ⇤ i with minimal geodesic curve length associated to the metric P, among all the four possible geodesics.

Let us denote the initial source point and the end point of the geodesic C ⇤ i by ā⇤ and c⇤ i , respectively. If the geodesic curve length is estimated by the fast marching method (Mirebeau, 2014b), this procedure can be simplified as follows: starting the fast marching front propagation from both of the initial source points s and s † , the orientation lifted point

c⇤ i 2 {p i , p † i }
is the first point that is reached by the front. The desired geodesic C ⇤ i can be determined by solving the ODE (21). As a result, a set {C ⇤ i ; 1  i  n} of all the desired geodesics can be extracted from the same minimal action map generated by a single fast marching propagation.

In these applications, the geodesic distance maps associated to the Finsler elastica metric are computed in a manner of early abort, i.e., once all the specified orientation lifted endpoints are reached by the fast marching front, the geodesic distance computation will be terminated. This early abort trick can greatly reduce the computation time. It is similar to the partial front propagation described in [START_REF] Deschamps | Fast extraction of minimal paths in 3D images and applications to virtual endoscopy[END_REF] with a simple extension to multiple points.

Experimental Results

We show the advantages of using curvature penalization for minimal paths extraction in the following experiments involving a study of the proposed metric itself, and comparative results against the isotropic Riemannian (IR) metric, the anisotropic Riemannian (AR) metric and the isotropic orientation lifted Riemannian (IOLR) metric in the applications of closed contour detection and tubular structure extraction. Note that in this section, whenever we mention Finsler elastica metric, we mean the data-driven Finsler elastica metic.

Riemannian Metrics Construction

We adopt the color image gradient introduced by Di Zenzo (1986) to construct the IR and AR metrics for closed boundaries detection of objects in color images.

Considering a color image I = (I 1 , I 2 , I 3 ) and a Gaussian kernel G with variance , the respective Gaussiansmoothed x-and y-derivatives of I are defined by

I x := @ x G ⇤ I and I y := @ y G ⇤ I, ( 66
)
where I x (•) and I y (•) should be understood as 1 ⇥ 3 vectors. Following [START_REF] Zenzo | A note on the Gradient of a Multi-Image[END_REF], a tensor E(x) can be defined by

E(x) = 0 @ kI x (x)k 2 hI x (x), I y (x)i hI x (x), I y (x)i kI y (x)k 2 1 A , 8 x 2 ⌦.
It is known that the tensor E(x) can be decomposed in terms of its eigenvalues and eigenvectors for all x 2 ⌦ by

E(x) = ' 1 (x) g 1 (x) g T 1 (x) + ' 2 (x) g 2 (x) g T 2 (x),
where ' 1 (x), ' 2 (x) are the eigenvalues of E(x) and g 1 (x), g 2 (x) are the associated eigenvectors. Without loss of generality, we assume that ' 1 (x)  ' 2 (x), 8 x 2 ⌦. In this case, ' 2 denotes the color gradient magnitude and g 2 denotes the unit color gradient vector field.

Based on the scalar field ' 2 , the IR metric F I can be constructed for all x 2 ⌦ by

F I (x, u) = ⇣ 1 + 2 ' p 2 (x) ⌘ 1 kuk, ( 67
)
where 1 , 2 and p are positive constants. In the following relevant experiments, we set 1 = 1 and p = 2.

The tensor field M A for the AR metric F A (8) can be computed by

M A (x) = exp(⌧ ' 2 (x)) g 1 (x) g T 1 (x) + exp(⌧ ' 1 (x)) g 2 (x) g T 2 (x), 8x 2 ⌦, ( 68
)
where ⌧ is a negative constant.

In the tubular structure extraction experiments, for the construction of the IR metric F I (67), we simply replace the scalar field ' 2 by the vesselness map V n defined in (57).

Moreover, regarding the construction of the AR metric, we make use of a radius lifted tensor field as introduced by Benmansour and [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], instead of using the tensor field M A defined in (68). In this case, the AR metric is regarded as the anisotropic radius lifted Riemannian (ARLR) metric defined over the radius lifted domain. In the following related experiments, we use the optimally oriented flux filter [START_REF] Law | Three Dimensional Curvilinear Structure Detection Using Optimally Oriented Flux[END_REF] to compute the ARLR metric as suggested in [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF]. For further details on the construction of the ARLR metric, we refer the reader to [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF].

The speed function IL for the IOLR metric F IL (11) should be dependent of the orientations. Simply, we compute the speed function IL by

IL (x, ✓) = (x, ✓), 8 ✓ 2 [0, ⇡), 8 x 2 ⌦, (69)
where is the orientation-dependent speed function defined in ( 58) and ( 59). The parameter ⇢ of the IOLR metric F IL penalizes the variations of the orientation ✓, and is set as

⇢ = ↵,
where ↵ is the parameter for the curvature term in the bending energy L (26).

Parameters Setting

Curvature penalization in the proposed Finsler elastica metric relies on two parameters, ↵ and (37). The choice of is dictated by algorithmic compromises. Indeed, minimal paths with respect to the Finsler elastica metric P converge to the elastica curves in the limit ! 1, hence a large value of is desirable. However, large values of yield metrics with strong anisotropy ratio µ(P). As a result, the numerical algorithm used in this paper, adapted from Mirebeau (2014b), uses larger discretization stencils, which increases its numerical cost and reduces its locality. For instance, = 30 (resp. 100 or 300) leads to stencils with a radius of 4 pixels (resp. 8 or 13). We typically use = 100.

On the other hand, the parameter ↵ is used to weight the curvature penalty in the Finsler elastica metric P. In the course of the fast marching method, a large value of ↵ makes the front to propagate slowly along the orientation dimension, implying that the obtained geodesics tend to be smooth, i.e., with low curvature. When ↵ is very small, the extracted geodesics mainly depend on the image data-driven speed functions defined in Section 4.3. Therefore, the choice of ↵ should depend on the desired image features. Basically, we make use of the following heuristics. There is a natural candidate ↵ ⇤ for the parameter ↵, dictated by the physical units of the parameters, namely ↵ ⇤ = (R ⇤ / ⇤ ) 2 , where R ⇤ is the smallest radius of curvature of the image features to be extracted, measured in pixels, and ⇤ is the typical value of the speed function around these features.

The angular resolution is set as ✓ s = ⇡/36 for both the IOLR metric and the proposed Finsler elastica metric. The parameter ⌘ for image data-driven speed function is set for each tested image individually. The parameter 2 that is used in the IR metric F I (67) is set as 2 = 2⌘ for all the comparative experiments. Unless otherwise specified, we set the anisotropy ratio values for the AR metric and the ARLR metric to be 20.

Smoothness and Asymmetry of the Finsler Elastica Minimal Paths

The Finsler elastica metric invoking orientation lifting and curvature penalty benefits from the smooth and asymmetric properties of the minimal paths. We demonstrate these properties in a synthetic image as shown in Fig. 5, where two ellipse-like shapes cross each other. The red dots and green dots are the initial source and end positions respectively. The arrows indicate the tangents at the corresponding positions. One can see that for the fixed initial source and end positions, changing the corresponding tangents will give different minimal paths. Moreover, the minimal paths with the same initial source and end positions could form a complete ellipse.

In Fig. 6, we design a spiral that has high anisotropy. The initial source position and end position are placed at the ends of the spiral. In the top row we add high noise to the spiral while in the bottom row we blur the spiral. In columns 1-4, we show the minimal paths obtained from the IR metric, the ARLR metric, the IOLR metric and the Finsler elastica metric, respectively. One can see that by using the mentioned Riemannian metrics, the shortcuts occur as shown in columns 1 to 3. In the top row, the minimal path (shown in column 3) obtained by using the IOLR metric is improved compared to the respective paths from the IR metric and the AR metric. However, a segment of the spiral is also missed due to the shortcuts problem. In contrast, the minimal paths shown in column 4 extracted by the Finsler elastica metric can completely avoid the problem of shortcuts thanks to the curvature penalization embedded in the proposed metric. In this experiment, we use an anisotropy ratio value of 100 for the ARLR metric. For the Finsler elastica metric, we set ↵ = 500 to ensure the geodesics to be smooth enough. s q 1 q 2 q 3 q 4 q 5 q 6 (a) In Fig. 7a, we specify six orientation lifted candidates qi , i = 1, 2..., 6, which are denoted by green dots with arrows, and an initial source point s (red dot with arrow). Among all of these candidates, we aim to find the closest orientation lifted candidate to the initial source point s, in terms of the geodesic distance associated to the metric P defined in (47). In Fig. 7b, it is shown that the closest orientation lifted point to s is the candidate q6 , even though the geodesic (red curve), joining the orientation lifted points s and q6 , passes through the vicinity of the physical position of q1 . Moreover, one can claim that the Euclidean distance value between the physical positions of s and q6 is the largest among all of the values between the physical positions of s and each orientation lifted candidate qi . This experiment demonstrates the asymmetric and smooth properties of the proposed Finsler elastica minimal path model.

In Fig. 8, we show the minimal path extraction results on three natural images, where each pair of the prescribed initial source positions and end positions is very close to each other in terms of the Euclidean distance. For each image, we expect to detect a long boundary between the two given orientation lifted points. It can be observed that the extracted geodesics associated to the Finsler elastica metric are able to catch the desired boundaries. In Fig. 8, the images shown in columns 1 and 2 are from the Berkeley Segmenta-tion Dataset [START_REF] Arbelaez | Contour Detection and Hierarchical Image Segmentation[END_REF] and the image in column 3 is from the Weizmann dataset [START_REF] Alpert | Image segmentation by probabilistic bottom-up aggregation and cue integration[END_REF].

6.4 Closed Contour Detection and Image Segmentation Fig. 9 shows the closed contour detection results with three prescribed physical positions using different metrics, where each position are assigned two opposite orientations1 . In this experiment, we firstly build the collection (64) by the proposed contour detection procedure using the Finsler elastica metric as described in Section 5.1, where the detection results are shown in column 5. Columns 2-4 show the closed contour detection results using the IR metric, the AR metric, and the IOLR metric, respectively. The minimal paths shown in columns 2 to 4 are obtained by simply linking each pair of vertices by the respective metrics involved in . The red, yellow, and green dots are the physical positions of the vertices q1 , q2 , q3 , respectively. The arrows shown in column 5 indicate the tangents of the geodesics at the corresponding positions. We assign each geodesic the same color as its initial source position. In these images, most parts of the desired boundaries appear to be weak edges which can be observed from the edge saliency maps from column 1. The detected contours associated to the Finsler elastica metric succeed at catching the desired boundaries due to the curvature penalization and asymmetric property. In contrast, the three Riemannian metrics without curvature penalization fail to extract the expected boundaries. The images used in this experiment are from the Weizmann dataset.

In Fig. 10, we show the closed contour detection results obtained by the proposed method with only two given physical positions and the corresponding orientations. One can see that the proposed method can indeed reduce the user intervention at least for objects with smooth boundaries.

For the Finsler elastica metric P defined in (47), the curvature penalization depends on the parameter ↵ ( is fixed to 100). In Fig. 11, we show the closed contour detection results by varying ↵ to demonstrate the influence of the curvature term in our approach. In column 1, we show the closed contour detection results with suitable values of ↵, say ↵ 0 . In columns 2 and 3, the closed contour detection results using ↵ 0 /10 and 5↵ 0 are demonstrated. One can see that it could lead to shortcuts by using small values of ↵ in rows 1-3 of column 2. In contrast, with a larger ↵, the detected closed contour can catch the optimal boundaries of the objects, which supports the effect of using curvature penalization. The edge saliency maps for each image in this experiments can be found from the first column of Fig. 9.

Perceptual Grouping

The perceptual grouping result on a synthetic noisy image is shown in Fig. 12. In Fig. 12a, we demonstrate the original image, which consists of a set of edges. The red and blue dots with arrows are the orientation lifted points provided by the user as initializations, where the red dot is the selected initial physical position. 65). In Fig. 20, the corresponding tangents for the physical positions are provided manually because high noise could lead to failure of the optimal orientation detection using (65). We use the extraction strategy described in Section 5.3 for the Finsler elastica metric.

In Fig. 15, the retinal vessels are extracted by the IR metric, the ARLR metric, the IOLR metric and the Finsler elastica metric as shown in columns 1 to 4, respectively. In columns 1 to 3, the minimal paths suffer from the short branches combination problem. In other words, these paths pass through some unexpected vessel segments. In contrast, the minimal paths obtained by the proposed metric can determine correct combinations of vessel branches. Similar extraction results are observed in Fig. 16. Again, the short branches combination occurs in columns 1 to 3, where the minimal paths in these columns are extracted by the IR metric, the ARLR metric and the IOLR metric, respectively. Instead, the Finsler elastica minimal paths can follow the correct vessel segments as shown in column 4.

In Fig. 17, we present the extraction results of the retinal artery centerlines in three patches of retinal images2 . The centerline of a retinal artery usually appears as a smooth curve. In column 1, we show the retinal artery-vein ground truth maps, where the red and blue regions indicate the arteries and veins, respectively. Note that the small vessels have been removed from the ground truth maps. Columns 2 to 5 show the centerlines extraction results by using the IR metric, the ARLR metric, the IOLR metric and the Finsler elastica metric, respectively. One can see that the minimal paths demonstrated in columns 2 to 4 pass through the wrong vessels due to the low gray-level contrast of the retinal arteries. The proposed model can obtain the correct artery centerlines as shown in column 5, thanks to the curvature penalization.

In Fig. 18, we demonstrate the retinal vein extraction results (blue curves) in the same patches which are used in rows 1-2 of Fig. 17. The extracted minimal paths (blue curves) are shown in columns 2 to 5 by using the IR metric, the ARLR metric, the IOLR metric and the Finsler elastica metric, respectively. From this figure, we can see that all of the extracted minimal paths can successfully follow the retinal veins as indicated by red regions in the artery-vein ground truth maps shown in column 1.

In Fig. 19, the vessel extraction results on a patch of a retinal image are demonstrated, where the vessels are blurred by other tissues. The minimal paths using the three Riemannian metrics fail to extract the desired vessel as shown in columns 1 to 3. In contrast, the Finsler elastica minimal paths can successfully delineate the targeted vessel as shown in column 4, thanks to the curvature penalization.

In Fig. 20, we show the road segmentation results on an aerial image by the proposed Finsler elastica metric. The road images are blurred by Gaussian noise with different variances. One can claim that our method is able to obtain smooth and accurate minimal paths on noisy images.

Conclusion

The core contribution of this paper lies at the introduction of the curvature penalization to the Eikonal PDE-based minimal path framework, which is accomplished by establishing the connection between the Euler elastica bending energy and the geodesic energy via a family of orientation lifted Finsler elastica metrics. Solving the Eikonal PDE with respect to the Finsler elastica metric, the minimal path model thus can determine globally minimizing curves that blend the benefits of the orientation lifting and the curvature penalization. Given a set of the user-provided orientation lifted points, we have demonstrated the ability of the proposed model for interactive image segmentation, perceptual grouping and tubular structure extraction. The experimental re- 25) but not frozen. The Accepted points are these points that have been updated and frozen. This method can solve the fixed point system (24) in a single pass manner. Compute U new (y) using ( 25). 5:

Algorithm 1 Fast Marching Method

Set U (y) min{U new (y), U (y)}. 6: end for Appendix B: Elements of Proof of Finsler Elastica Minimal Paths Convergence Let X ⇢ R d be a compact domain and = be the collection of compact convex sets of R d . We will later specialize to d = 3 for the application to the Euler elastica curves. The set = is metric space, equipped with the Haussdorff distance which is defined as follows.

Definition 1 The Euclidean distance map from a set

A ✓ R d is d(A, x) := inf y2A kx yk. ( 70
)
The Haussdorff distance between sets

A 1 , A 2 j R d is H(A 1 , A 2 ) := sup x2R d |d(A 1 , x) d(A 2 , x)|. ( 71
)
Definition 2 Let 2 C 0 ([0, 1], X) be a path and B 2 C 0 (X, =) be a collection of controls on X. The path is said B-admissible iff it is locally Lipschitz and 0 (t) 2 B( (t)) for a.e.

) (73
where C is the continuous cutoff function defined by

C(t) = 8 > < > : 1 if t  1 2 t if 1  t  2 0 if t 2.
Clearly B 7 ! diam(B) and (B, ✏) 7 ! ⌅(B, ✏) are continuous functions of B 2 C 0 (X, =) and ✏ 2 R + . In addition ⌅(B, ✏) is increasing w.r.t. ✏. Here and below, if A 1 , A 2 are metric spaces, and A 1 is compact, then C 0 (A 1 , A 2 ) is equipped with the topology of uniform convergence. This applies in particular to the space of paths C 0 ([0, 1], X) and of controls C 0 (X, =).

Lemma 1 If is B-admissible, then its Lipschitz constant is at most diam(B). A necessary and sufficient condition for to be B-admissible is: for all 0  p  q  1, d ✓ B( (p)), (q) (p) q p ◆  ⌅(B, (q p) diam(B)),

Proof Assume that is B-admissible. Then for any 0  p  q  1 one has k (p) (q)k  Z q p k 0 (%)kd%  |p q| diam(B),

hence is diam(B)-Lipschitz as announced. Denoting w % = p + (q p)%, for all % 2 [0, 1], one obtains (q) (p) q p = Z 1 0 0 (w % )d%. Conversely, assume that and B obey (74). Then (q) (p) q p  diam(B) + ⌅(B, diam(B)), for any 0  p  q  1. Thus is a Lipschitz path as announced, and therefore it is almost everywhere differentiable. If p 2 [0, 1] is a point of differentiability, then letting q ! p we obtain d(B( (p)), 0 (p))  ⌅(B, 0) = 0, which is the announced admissibility property 0 (p) 2 B( (p)). u t

The characterization (74) is written in terms of continuous functions of the path and control set B, hence it is a closed condition, which implies the following two corollaries. We denote T B(x) := {T u; u 2 B(x), 8x 2 X}.

Corollary 1 The set {( , B) 2 C 0 ([0, 1], X) ⇥ C 0 (X, =); is B-admissible} is closed.

Corollary 2 Let x n , y n , B n and T n be converging sequences in X, X, C 0 (X, =) and R + , with limits x 1 , y 1 , B 1 , and T 1 respectively. Let n 2 C 0 ([0, 1], X) be a (T n + 1/n)B n -admissible path with endpoints x n and y n . Then the sequence of paths ( n ) is equip-continuous, and the limit 1 of any converging subsequence is a T 1 B 1 -admissible path 2 C 0 ([0, 1], X) with endpoints x 1 and y 1 .

Proof Note that the map (T, B) 7 ! T B is continuous on R + ⇥ C 0 (X, =), hence the controls Bn := (T n + 1/n)B n converge to B1 := T 1 B 1 . Defining E := sup{diam( Bn ); n > 0} which is finite by continuity of diam(•) (72) and convergence of Bn , we find that the paths ( n ) n>0 are simultaneously E-Lipschitz, hence that a subsequence uniformly converges to some path 1 . The B1admissibility of 1 then follows from Corollary 1. u t

We next introduce the minimum-time optimal control problems. The minimum of ( 80) is attained by Corollary 2, which also immediately implies the Corollary 3.

Definition 4 For all x, y 2 X, and B 2 C 0 (X, =), we let T B (x, y) := min{T > 0; 9 2 C 0 ([0, 1], X), (0) = x, (1) = y, is T B-admissible}. (80)

Corollary 3 The map (x, y, B) 7 ! T B (x, y) is lower semi-continuous on X ⇥ X ⇥ C 0 (X, =). In other words, whenever (x n , y n , B n ) ! (x 1 , y 1 , B 1 ) as n ! 1 one has

lim inf T B n (x n , y n ) T B 1 (x 1 , y 1 ). ( 81
)
Proof For each n > 0 let T n = T B n (x n , y n ) and let T 1 = lim inf T n as n ! 1. Up to extracting a subsequence, we can assume that T 1 = lim T n as n ! 1. Denoting by n a path as in Corollary 2 for all n 0, we find that there is a converging subsequence which limit 1 is T 1 B 1 admissible and obeys 1 (0) = x 1 and 1 (1) = y 1 . Thus T 1 T B 1 (x 1 , y 1 ) as announced.

u t

Definition 5 Let B 1 , B 2 2 C 0 (X, =). These collections of controls are said included B 1 ✓ B 2 iff B 1 (x) ✓ B 2 (x) for all x 2 X.

The property B 1 ✓ B 2 clearly implies, for all x, y 2 X T B 1 (x, y) T B 2 (x, y).

(82)

Corollary 4 Assume that one has a converging sequence of controls B n ! B 1 obeying the inclusions B n ◆ B 1 for all n 0. Then lim T B n (x, y) = T B 1 (x, y) (83) for all x, y 2 X. Let T n := T B n (x, y) for all n 2 N [ {1}, and let n be an arbitrary (T n + 1/n)B n -admissible path from x to y.

If there exists a unique T 1 B 1 -admissible path 1 from x to y, then n ! 1 as n ! 1.

Proof The identity (83) follows from inequalities ( 81) and (82). By Corollary 2 the sequence of paths n is equi-continuous, and any converging subsequence tends to a T ⇤ B 1 admissible path ⇤ from x to y, with T ⇤ := lim T n . Since T ⇤ = T 1 and by uniqueness we have ⇤ = 1 , hence n ! 1 as announced.

Application to Finsler elastica geodesics convergence problem.

Consider an orientation lifted Finsler metric F : X⇥R 3 ! R + where X := ⌦ ⇢ R 3 . For any x 2 X, let B(x) := {ū 2 R 3 ; F (x, ū)  1} be the unit ball of F . B(x) is compact and convex, due to the positivity, continuity and convexity of the metric F and the map B : X ! = is continuous. Furthermore, using the homogeneity of the metric, one obtains for all x, ȳ 2 X:

L ⇤ F (x, ȳ) = T B (x, ȳ),
where L ⇤ (x, ȳ) is the minimal curve length between x and ȳ with respect to metric F . In the case of Finsler elastica problem, one has B 1 (x) := B 1

x and B (x) := B x , 8x 2 X, where B 1 x and B x are defined in equations ( 42) and ( 43) respectively. The Finsler elastica metrics F on X pointwisely tend to the metric F 1 as ! 1. Fortunately, the associated control sets B (x) ! B 1 (x) uniformly in C 0 (⌦, =), as can be seen from ( 46). Hence one has lim inf L ⇤ F (x, ȳ) = lim inf T B (x, ȳ) T B 1 (x, ȳ) = L ⇤ F 1 (x, ȳ),

as ! 1 for all x, ȳ 2 X. To show that equality holds, it suffices to prove that sequence B obeys B ◆ B 1 , equivalently to prove that F (x, ū)  F 1 (x, ū) for all x 2 X and any vector ū 2 R 3 . Indeed, let x = (x, ✓) 2 X, and ū = (u, ⌫) 2 R 2 ⇥ R:

F (x, ū) = p 2 kuk 2 + 2 |⌫| 2 ( 1)hu, v ✓ i = kuk 1 + s 1 + 2|⌫| 2 kuk 2 ! + kuk ( 1)hu, v ✓ i. = kuk + |⌫| 2 kuk 0 @ 2 1 + q 1 + 2|⌫| 2 kuk 2 1 A + ( 1)(kuk hu, v ✓ i). (84)
 F 1 (x, ū).

The last inequality holds because the denominator 1 +

q 1 + 2|⌫| 2 kuk 2
in (84) is greater than 2, and kuk hu, v ✓ i for any vector u and any orientation ✓. By Corollary 2, minimal paths C with endpoints x and ȳ for geodesic distance L ⇤ F (x, ȳ) converge as ! 1 to a minimal path C 1 for L ⇤ F 1 (x, ȳ). We finally point out that L ⇤ F 1 (x, ȳ) < 1 for all x, ȳ in the interior of X, provided this interior is connected, due a classical controllability result for the Euler elastica problem.

Fig. 1

 1 Fig. 1 (a) Edge saliency map. (b) Minimal path with isotropic Riemannian metric. (c) Minimal path with anisotropic Riemannian metric. (d) Minimal path with the proposed Finsler elastica metric. The red curves are the extracted minimal paths with the initial source positions and end positions indicated by the red and green dots respectively. The arrows in (d) indicate the tangents.

Fig. 2

 2 Fig. 2 Visualization for the metrics F 1 and F with ↵ = 1 by Tissot's indicatrix. (a) Tissot's indicatrix for the metric F 1 (37) with ↵ = 1 are flat 2D disks embedded in 3D space, aligned with the direction v ✓ (several directions ✓ are shown). (b) Tissot's indicatrix for the Finsler elastica metrics F are ellipsoids, which flatten and approximate a limit disk as ! 1.

Fig. 3

 3 Fig. 3 Approximating Euler elastica curves by Finsler elastica minimal paths. (a) Finsler elastica minimal paths with ↵ = 500 and different values of . (b) and (c) Finsler elastica minimal paths with = 100 and = 300 respectively, and different values of ↵.

  )) and different values of . The arrows indicate the initial source and end points tangents. The cyan point denotes the initial position and the

Fig. 4

 4 Fig. 4 Steps for the closed contour detection procedure. (a) Original image and all of the vertices in D denoted by dots and arrows. (b) The first pair of successive vertices (q 1 , q2 ) is detected. (c) The third vertice q3 is detected. (d) The final vertice q4 is detected and the closed contour detection procedure is terminated. (e) The minimal path joining q4 and q1 is tracked (cyan curve). (f) The final closed contour is obtained.

Fig. 5

 5 Fig. 5 Flexible Finsler elastica minimal paths extraction on ellipse-like curves. The red and green dots denote the source and end positions, respectively. Arrows indicate the tangents.

Fig. 6

 6 Fig. 6 Minimal paths extraction results on Spirals. Columns 1-4 Minimal paths extracted by the IR metric, the ARLR metric, the IOLR metric and the Finsler elastica metric, respectively. The red and green dots denote the source and end positions, respectively. The arrows indicate the tangents.

Fig. 7

 7 Fig. 7 Finding the closest orientation lifted candidate to the source point s in terms of the geodesic distance associated to the Finsler elastica metric. (a) The orientation lifted candidates qi , i = 1, 2, • • • , 6, and source point s are demonstrated. (b) The closest candidate q6 is detected. Red curve indicates the geodesic linking s to q6 .

Fig. 8

 8 Fig. 8 Finsler elastica minimal paths extraction results. The red and green dots indicate the source and end positions, respectively. The arrows indicate the corresponding tangents.

Fig. 9

 9 Fig. 9 Comparative closed contour detection results obtained by using different metrics. Column 1 Edge saliency map. Columns 2-5 The closed contour detection results from the IR metric, the AR metric, the IOLR metric and the Finsler elastica metric, respectively.

Fig. 10

 10 Fig. 10 Closed contour detection results by using only two given physical positions and the corresponding orientations.

Fig. 11

 11 Fig. 11 Closed contour detection results with different values of ↵. Column 1 Results by the suitable values of ↵ and . Columns 2 Results by small values of ↵. Column 3 Results by large values of ↵.

  Fig. 12b shows the perceptual grouping results of the proposed method. The identified orientation lifted points in the set D 1 are denoted by red dots with arrows. The red curves that links the ordered orientation lifted vertices in D 1 indicate the expected closed curves.

Fig. 12 (

 12 Fig. 12 (a) Initialization. The red and blue dots indicate the physical positions, where the red dot is the initial position. (b) Perceptual grouping result. The arrows denote the tangents for each physical position.

Fig. 13

 13 Fig. 13 illustrates the capacity of the proposed method to handle the grouping problem with spurious points. Different initializations are shown in Figs. 13a and 13c, where the red dots are the selected initial physical positions. Figs. 13b and 13d are the grouping results. The red curves indicate the detected closed curves. The proposed perceptual grouping method can detect multiple closed curves by specifying the number of expected curves. In Fig. 14, three curves are detected by the proposed method. Column 1 shows different initializations, where the red dots are the selected initial positions. Columns 2 to 4 demonstrate the intermediate grouping results that correspond to different initializations shown in column 1. The final perceptual grouping results are shown in column 5. One can claim that our algorithm indeed has the ability to detect curves intersecting with one another. The vertices shown in columns 2 to 4, denoted by red dots with arrows, make up the corresponding collections D 1 to D 4 , respectively.

Fig. 13

 13 Fig. 13 Perceptual grouping results. (a) Initialization 1. Red dot is the selected initial position. (b) Perceptual grouping result for initialization 1. (c) Initialization 2. (d) Perceptual grouping result for initialization 2.

Fig. 14

 14 Fig. 14 Perceptual grouping results by the proposed method where three groups are identified. Column 1 Initializations. The red dots denote the selected initial positions. Columns 2-4 Intermediate grouping results for the corresponding initializations. Column 5 Final grouping results.

Fig. 15

 15 Fig.15Comparative blood vessel extraction results on retinal images. Columns 1-4 The extracted minimal paths using the IR metric, the ARLR metric, the IOLR metric, and the proposed Finsler elastica metric, respectively.

Fig. 16

 16 Fig.16Comparative blood vessel extraction results on fluoroscopy images. Columns 1-4 The extracted minimal paths using the IR metric, the ARLR metric, the IOLR metric, and the proposed Finsler elastica metric, respectively.

Fig. 17

 17 Fig. 17Comparative blood vessels extraction results on retinal images. Column 1 The retinal artery-vein vessels ground truth maps. Columns 2-5 The extracted minimal paths by the IR metric, the ARLR metric, the IOLR metric, and the Finsler elastica metric, respectively.

  Fig. 17Comparative blood vessels extraction results on retinal images. Column 1 The retinal artery-vein vessels ground truth maps. Columns 2-5 The extracted minimal paths by the IR metric, the ARLR metric, the IOLR metric, and the Finsler elastica metric, respectively.

Fig. 18

 18 Fig. 18 Retinal veins extraction results from different metrics. Column 1 The retinal artery-vein vessels ground truth maps. Columns 2-5 The extracted minimal paths (blue curves) by the IR metric, the ARLR metric, the IOLR metric, and the Finsler elastica metric, respectively.

Fig. 19

 19 Fig. 19 Comparative blood vessel extraction results on a blurred retinal image. Columns 1-4 The extracted minimal paths by the IR metric, the ARLR metric, the IOLR metric, and the Finsler elastica metric, respectively.

Fig. 20

 20 Fig. 20 Roads extraction results by the proposed Finsler elastica metric in an aerial image blurred by Gaussian noise.

  Appendix A: Fast Marching Method Basically, the fast marching method introduces a boolean map b : Z ! {Trial, Accepted} and tags each grid point x 2 Z either Trial or Accepted. The Trial points are defined as grid points which have been updated at least once by Hopf-Lax operator (

  Input: Metric F , stencil map S, initial source points collection W. Output: Minimal action map U .Initialization:For each point x 2 Z, set U (x) +1 and b(x) Trial. For each point y 2 W, set U (y) 0. While: at least one grid point is tagged as Trial 1: Find x min , the Trial point which minimizes U . 2: b(x min )Accepted. 3: for All y such that x min 2 S(y) and b(y) 6 = Accepted do 4:

  t 2 [0, 1]. Definition 3 A collection of controls on X is a map B 2 C 0 (X, =). Its diameter diam(B) and modulus 3 of continuity ⌅(B, ✏), are defined by ⌅(B, 0) = 0 and for ✏ > 0 diam(B) := sup{kuk; u 2 B(x), 8x 2 X}, (72) ⌅(B, ✏) := sup x,y2X H(B(x), B(y))C(kx yk/✏),

  Hence by Jensen's inequality and the convexity of d(B( (p)), •), which follows the convexity of B( (p)), of the announced characterization. The inequality (78) follows from the admissibility property 0 (w % ) 2 B( (w % )) and the definition of the Haussdorff distance (71). The inequality (79) follows from the above established Lipschitz regularity of and the definition of the modulus of continuity (73).

Table 1

 1 Computation time (in seconds) and average number of Hopf-Lax updates required for each grid point by the fast marching method with ↵ = 500 and different values of on a 300 2 ⇥ 108 grid.

		1	10	20	30	100	200	1000
	time	13.9s 25.3s 27.3s 27.7s 31.7s 33.9s 36.8s
	number	3	5.49	6.06	6.49	7.27	7.82	8.12

For the IR metric and the AR metric, only the physical positions of these orientation lifted vertices are used.

Many thanks to Dr. Jiong Zhang to provide us these images.

We actually use a slight variant of the classical modulus of continuity because the latter one, obtained with the hard cutoff function C(t) = 1 if t  1, 0 otherwise, lacks continuity in general.

Acknowledgements The authors would like to thank all the anonymous reviewers for their detailed remarks that helped us improve the presentation of this paper. This work was partially supported by ANR grant NS-LBR, ANR-13-JS01-0003-01.