Jean-Marie Mirebeau

Jorg Portegies

Hamiltonian Fast Marching

We introduce a generalized Fast-Marching algorithm, able to compute paths globally minimizing a measure of length, defined w.r.t. a variety of metrics in dimension two to five. Our method applies in particular to arbitrary Riemannian metrics, and implements features such as second order accuracy, sensitivity analysis, and various stopping criteria. We also address the singular metrics associated with several non-holonomic control models, related with curvature penalization, such as the Reeds-Shepp's car with or without reverse gear, the Euler-Mumford elastica curves, and the Dubins car. Applications to image processing and to motion planning are demonstrated.

Introduction

In this paper we provide the numerical details involved in a generalized fast marching solver, able to compute the distance from a point within a domain, and the related minimal paths. Path length is measured w.r.t. a given metric which may take a variety of forms and involve data-driven parameters. For instance the metric may be Isotropic (locally proportional to the Euclidean norm), or Riemannian. Combining a dimension lifting trick with the choice of an adequate singular metric, we are also able to penalize path curvature, according to the Reeds-Shepp [START_REF] Reeds | Optimal paths for a car that goes both forwards and backwards[END_REF][START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF], Euler-Mumford [START_REF] Mumford | Elastica and computer vision[END_REF], or Dubins [START_REF] Dubins | On Curves of Minimal Length with a Constraint on Average Curvature, and with Prescribed Initial and Terminal Positions and Tangents[END_REF] models. Our software is referred to as the Hamiltonian Fast Marching (HFM) library 1 , because it relies on an original and specific representation of Hamiltonian of the addressed problems (13), and uses the fast marching method (a generalization of Dijkstra's algorithm) to solve eikonal-type equations in a single pass.

The search for optimal paths appears in many applications, such as motion planning and image analysis, see the numerical section of this paper and e.g. [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF]. The numerical approach chosen in this paper is to first compute a distance map, which is the viscosity solution to a Partial Differential Equation (PDE) of eikonal type [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF], and then extract the minimal paths, which obey an Ordinary Differential Equation (ODE) defined in terms of the distance map. In the case of isotropic metrics, the first task is essentially a solved problem [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], however it becomes challenging when considering (strongly) anisotropic metrics [KS98, SV01, BR06, AM11, [START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF]. For that purpose, we rely on an original Eulerian and causal discretization of the eikonal equation [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF][START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF]. Two robust geodesic backtracking methods are provided for the second task. The HFM software is made available as an opensource C++ library with interfaces to Python R , Matlab R , and Mathematica R . It is fast (albeit single-threaded by nature), built for extensibility (for introducing additional classes of metrics, stopping criteria, interfaces with other programming languages, etc), features uncommon but useful methods such as reverse mode sensitivity analysis, and comes with a thorough series of introductory notebooks.

Curve optimization via eikonal PDEs

In order to make our contributions clear and to relate them to the literature, we need to formally state the addressed mathematical problem. Let E = R d be the ambient space, and let Ω ⊂ E be a domain. Our numerical method requires Ω to be box shaped, possibly equipped with periodic boundary conditions, and with inner obstacles. The geometry of the domain is described using a quasi-metric F : Ω × E → [0, ∞], which specific form is discussed later. This object defines a measure of path length L F : Γ → [0, ∞], where Γ := Lip([0, 1], Ω) is the set of locally Lipschitz paths, and a related quasi-distance d F : Ω × Ω → [0, ∞], as follows. For any path γ ∈ Γ, and any points p, q ∈ Ω L F (γ) := Depending on the choice of quasi-metric F the quasi-distance d F may or may not satisfy d F (p, q) = d F (q, p) (symmetry), d F (p, q) < ∞ (global controllability), or d F (p, q) → 0 as p → q (local controllability). For readability we choose to drop in the following the prefix "quasi-", used conventionally to emphasize that a metric or a distance may lack symmetry. Note that we reserve the word metric for functions F measuring the length of tangent vectors, and the word distance for functions such as d F measuring the minimal path length between two points.

Our objective is to compute the distance map u : Ω →] -∞, ∞] from the domain boundary, which is defined as follows ∀p ∈ Ω, u(p) = inf q∈∂Ω u(q) + d F (q, p), ∀p ∈ ∂Ω, u(p) = σ(p),

where σ : Ω →] -∞, ∞] is a given initial delay. The distance from e.g. a single point of interest q * ∈ ∂Ω is obtained by setting σ(q * ) = 0 and σ ≡ ∞ on ∂Ω \ {q * }. Under suitable assumptions, the function u is known to be the unique viscosity solution to the (generalized) eikonal PDE ∀p ∈ Ω, F * p (du(p)) = 1, ∀p ∈ ∂Ω, u(p) = σ(p).

(3)

We refer to [START_REF] Bardi | Optimal Control and Viscosity Solutions of Hamilton-Jacobi-Bellman Equations[END_REF] for this theory, including the concept of discontinuous solutions to eikonal equations, which is required for some of our problem instances [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] but is out of the scope of this paper. The eikonal PDE (3) involves the dual metric F * : Ω × E * → [0, ∞[, defined by

F * p (p) := sup{ p, ṗ ; ṗ ∈ E, F p ( ṗ) ≤ 1}. ( 4
)
Here and below, points p ∈ Ω are denoted with bold letters, tangent vectors ṗ ∈ E are distinguished with dots, and co-vectors p ∈ E * are decorated with hats. Once u is known, the minimal path joining a given q ∈ Ω from the corresponding optimal p ∈ ∂Ω is extracted by solving (backwards in time) the following Ordinary Differential Equation (ODE):

∀t ∈ [T -, T + ], γ(t) = dF * γ(t) (du(γ(t))), γ(T + ) = q, (5)
where T -:= σ(p) and T + := u(q). We denoted by dF * p (p) the differential of p ∈ E * → F * p (p), the dual metric at the point p ∈ Ω. See [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF] Appendix C.

We discuss in §1.2 possible choices of metric F that can be addressed with the HFM library, and in particular the singular metrics used to encode curvature penalization. The design choices underlying our numerical method for solving the PDE (3) are then presented in §1.3. Applications of minimal path methods to image processing and motion planning are briefly surveyed in §1.4, where we also present the outline of the paper.

Isotropic, Anisotropic and Non-Holonomic metrics

In the applications for which our software is intended, such as image processing and motion planning, the metric F encoding the geometry of the domain is usually data driven and spatially inhomogeneous. The specific structure of the expression of the metric is in principle dictated by the application, but it must also be compatible with the numerical method used for minimal path computation. We present in the following the metric models that can be addressed using our software, distinguishing three general classes: Isotropic, Anisotropic, and Non-Holonomic.

In each case, we provide the generic expression of the primal metric, as well as the dual metric which is defined by (4) and appears in the eikonal equation (3) and the geodesic backtracking ODE (5).

Isotropic metrics are locally proportional to the Euclidean norm, which is defined as usual on the embedding space E := R d ⊇ Ω. The dual metric turns out to be proportional to the Euclidean norm as well, but with an inverse ratio. For any point p ∈ Ω, any vector ṗ ∈ E and any co-vector p ∈ E *

The cost function c : Ω →]0, ∞[ appearing in these expressions, is given by the user and usually data driven. It is homogeneous to the inverse of a speed, and assumed to be continuous. The HFM software uses the classical discretization scheme [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF] to numerically solve isotropic eikonal equations on cartesian grids of arbitrary dimension. Anisotropic metrics define, at any given point, different norms for unit euclidean vectors in different directions. Riemannian metrics are the simplest and most common instance of Figure 2: A metric associates to each point a gauge, i.e. a generalized norm which may be nonsymmetrical and may take infinite values. The three first figures display unit balls (also referred to as Tissot's indicatrix) of the gauges of successively: an Isotropic metric, an Anisotropic Riemannian metric, and the non-holonomic sub-Riemannian metric associated with the Reeds-Shepp model on R 2 × S 1 . In the later case, the admissible directions of motion depend on the present angular coordinate θ, associated with the vertical axis, reflecting the fact that a car cannot move sideways. anisotropic metrics. They are determined by a field M : Ω → S ++ (E) of positive definite tensors, and take the form

F p ( ṗ) := ṗ M(p) , F * p (p) = p D(p) , (7)
where ṗ M := M ṗ, ṗ , and where D(p) := M(p) -1 denotes the dual tensor field. The HFM software uses the recently developed FM-VR1 numerical scheme [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], to solve Riemannian eikonal equations, in dimension d ∈ {2, 3}. Finsler metrics [START_REF] Zach | Globally Optimal Finsler Active Contours[END_REF][START_REF] Chen | Finsler Geodesics Evolution Model for Region based Active Contours[END_REF], i.e. non-Riemannian anisotropic metrics, are presently not supported by the HFM software. For these we refer to the earlier open source library [START_REF] Mirebeau | Anisotropic Fast Marching in ITK[END_REF], due to the first author, which implements the (two-dimensional) FM-ASR scheme [START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF].

Non-Holonomic control models forbid, at some points in space p ∈ Ω, some directions of motion ṗ ∈ E. The corresponding metrics are singular, in that they associate infinite costs F p ( ṗ) = +∞ to these forbidden directions. Sub-Riemannian metrics are the simplest instances of non-holonomic metrics, and can be regarded as generalized Riemannian metrics having some infinite eigenvalues. The HFM software is able to handle strongly anisotropic Riemannian metrics, with condition numbers 10, which can approximate sub-Riemannian metrics convincingly enough for many applications. We implement in particular the Reeds-Shepp car model, defined on the configuration space R 2 × S 1 of positions and orientations, which is adequate for modeling wheelchair-like vehicles [START_REF] Reeds | Optimal paths for a car that goes both forwards and backwards[END_REF] but is also fundamental in image vision [START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF]. Several generalisations, defined on the 5-dimensional domain R 3 × S 2 , see [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF][START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], are also present.

We address the computation of planar paths globally minimizing several curvature dependent energies, following [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF]. For that purpose, consider again the three-dimensional state space M := R 2 × S 1 of positions and orientations, which points are denoted p = (x, θ). A singular quasi-metric F is introduced, involving a data-driven cost function c : M →]0, ∞[ and a curvature penalty C : R →]0, ∞]. For any tangent vector ṗ = ( ẋ, θ), with unit physical velocity ẋ = 1 and arbitrary angular velocity θ, and any co-tangent vector (x, θ), one has

F (x,θ) ( ẋ, θ) := c(x, θ) • C( θ) if ẋ = n(θ), +∞ otherwise. F * (x,θ) (x, θ) = c(x, θ) -1 sup θ∈R x, n(θ) + θ θ C( θ) . (8)
We denoted n(θ) := (cos θ, sin θ), and adopted the convention that S 1 := R/(2πZ) has tangent space R. The quasi-metric F is extended by positive 1-homogeneity to the case ẋ = 1, and the dual metric expression follows from (4), see also Appendix A. Consider a path η : [0, T ] → R 2 , parametrized at unit euclidean speed, and abusively identify the unit vector η(s) with the corresponding direction in S 1 . Then the length of the orientation-lifted path (η, η) measured w.r.t. the above metric F equals the following curvature penalized energy

L F ((η, η)) = L 0 c(η(s), η(s)) C(η(s)) ds,
In addition, by construction of the metric F, any orientation lifting (η, η) : [0, T ] → R 2 × S 1 of the physical path η, distinct from the canonical one (η, η), satisfies L F ((η, η)) = +∞, hence cannot be optimal.

The state dependent cost c(x, θ) is arbitrary and user defined, but the curvature cost C( θ) must presently be chosen among the following instances, related to the Reeds-Shepp car without reverse gear [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF], the Euler-Mumford elastica curves [START_REF] Mumford | Elastica and computer vision[END_REF], and Dubins's car [START_REF] Dubins | On Curves of Minimal Length with a Constraint on Average Curvature, and with Prescribed Initial and Terminal Positions and Tangents[END_REF]. For any κ ∈ R, standing for the path curvature θ

C RS (κ) := 1 + (ξκ) 2 , C EM (κ) := 1 + (ξκ) 2 , C D (κ) := 1 if |ξκ| ≤ 1, +∞ else, ( 9
)
where the parameter ξ > 0 modulates the amount of curvature penalization and has the dimension of a turning radius. Note that the metric (8) associated with (9, left) defines a model, referred to as the Reeds-Shepp forward model, which distinguishes itself from the original sub-Riemannian Reeds-Shepp model [START_REF] Reeds | Optimal paths for a car that goes both forwards and backwards[END_REF][START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF] discussed in the previous paragraph by the absence of reverse gear, see the discussion in [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF]. Finally, we also provide the more parametrable models defined by the metrics of the form

F (x,θ) ( ẋ, θ) = c(x, θ) C ξ(x, θ)( θ -κ(x, θ)) if ẋ = n(θ). ( 10
)
As before, the metric is extended 1-homogeneously to the case where x is positively collinear with n(θ), and set to +∞ otherwise. The fields ξ : M →]0, ∞[ and κ : M → R, provided by the user, locally modulate the turning radius and imbalance of the car. See Figure 20.

An Eulerian and causal discretization of the eikonal equation

The HFM software relies on an Eulerian and causal discretization of the eikonal equation (3), on a cartesian grid, based on the papers [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF][START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF]. We briefly describe here these two design choices, and contrast them with the possible alternatives. Numerical solvers of the time optimal control problem (2) can be of either semi-Lagrangian or Eulerian nature. Semi-Lagrangian schemes are a discretization of Bellman's optimality principle, which expresses that the characterization (2, left) of the problem solution u still holds if the global domain boundary ∂Ω is replaced with the boundary ∂V (p) of an arbitrary (in practice chosen small) neighborhood V (p) ⊂ Ω of the current point p. A continued line of research has been devoted to this approach since the seminal work [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], see the paragraph on causality below. The HFM software does not follow this popular route, but relies in contrast on an Eulerian scheme, directly expressing the eikonal PDE (3) at the discrete level using finite differences, similarly to [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF][START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF]. This approach requires an exact or approximate representation of the dual metric in the following form: for any point p ∈ Ω, and any co-vector p ∈ E *

F * p (p) 2 ≈ 1≤i≤I α i (p) max{0, p, ėi } 2 , (11)
1≤i≤I α i (p) max{0, U (p) -U (p -h ėi )} 2 = h 2 , ∀p ∈ ∂X, U (p) = σ(p). (12)
where h > 0 is the discretization scale, and X and ∂X are subsets of the cartesian grid hZ d devoted to discretizing Ω and ∂Ω respectively. Note, crucially, that the numerical scheme (12, left) is a non-decreasing function of the finite differences (U (p) -U (p -h ė)) ė∈Z d . This property, referred to as monotony or degenerate ellipticity, is essential in establishing the well posedness and stability of the discretized problem, see [START_REF] Oberman | Convergent Difference Schemes for Degenerate Elliptic and Parabolic Equations: Hamilton-Jacobi Equations and Free Boundary Problems[END_REF]. The implementation of the HFM library is slightly more flexible than (12), allowing for e.g. axis dependent grid scales, the introduction of obstacles in the domain, and second order accuracy.

Causality is a property of the discretization scheme which reflects the deterministic nature of the original control problem, and enables it to be solved in a single pass using the fast marching algorithm/dynamic programming principle, see §3.1. In the context of Eulerian schemes, causality means that the numerical scheme is a function of the positive parts of the finite differences (U (p) -U (p -h ė)) ė∈Z d , which is indeed the case in (12, left). See Definition 2.1 below for a formal definition, and [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations[END_REF] for a counterpart in the context of semi-Lagrangian schemes. At the continuous level, u(p) can be regarded as the arrival time at p ∈ Ω of a front originating from the domain boundary ∂Ω at the given time σ, and propagating with a speed locally dictated by the metric F. The continuous counterpart of the causality property is the fact that the front arrival time at a given point only depends on the earlier arrival times, and not on the future ones. In the context of shortest paths on graphs, causality amounts to the positivity of the edge weights, which similarly enables Dijkstra's single pass algorithm. Note that monotone but non-causal discretization schemes can also be numerically solved, using a variety of iterative methods such as Fast-Sweeping [START_REF] Zhao | A fast sweeping method for eikonal equations[END_REF] or (adaptive) Gauss-Seidel iterations [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF][START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF]. Their practical usage can however be tedious, in particular with the strongly anisotropic and inhomogeneous metrics encountered in image processing, which lead to long and unpredictable computation times [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF].

Designing causal discretization schemes for anisotropic metrics is a non-trivial task, which has been the subject of a continued line of research in the semi-Lagrangian setting [KS98, JBT + 08, AM11], following the discovery [START_REF] Sethian | Ordered upwind methods for static Hamilton-Jacobi equations[END_REF][START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF] that it is related to a geometrical property of the discretization stencils. The works [START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF][START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF] assume a cartesian grid discretization and use, similarly to here but in the semi-Lagrangian setting, tools from lattice geometry to design reasonably small causal stencils for (strongly) anisotropic metrics. In the Eulerian setting, the design of causal schemes for non-isotropic metrics has in contrast been overlooked until the recent works [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF][START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF], which are at the foundation of the proposed software.

Applications to image processing and motion planning

Minimal path methods have numerous applications [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF], of which we can only give a glimpse. We focus here on applications which can directly benefit from the numerical methods provided within the proposed software, and in particular from paths globally minimizing a curvature dependent energy.

Motion planning. When preparing the motion of a automatic vehicle in a known environnement, or the movement of an articulated machine such as a robotic arm, using the path globally minimizing a well chosen energy has obvious advantages [START_REF] Latombe | Robot motion planning[END_REF]. A variety of methods for path optimization exist, such as local optimization from an initial guess, geodesic shooting, stochastic exploration, or fine geometrical descriptions [START_REF] Boissonnat | Shortest paths of bounded curvature in the plane[END_REF]. Approaches based on the eikonal PDE have uncomparable flexibility, and a guarantee of global optimality. Their main drawbacks are their cost for the complex models, due to the curse of dimension, and until recently the inability to handle non-holonomic constraints. The HFM library solves the latter issue, and for instance allows to introduce curvature penalization in the energy functional, so as to minimize the effort on the machine structure. Hard constraints can also be accounted for, such as a vehicles with a bounded turning radius. In the proof of concept work [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF], the HFM library is used to optimize a surveillance system for detecting an an enemy vehicle, subject to non-holonomic constraints, in a worst case scenario.

Image Segmentation. Numerous imaging applications involve paths, such as two dimensional region segmentation, tubular structure segmentation, or white fiber tracking in dMRI scans of the brain. Various path constructions exist, such as fiber following methods based on ODEs in the last case, often enhanced with stochastic perturbations [SKF + 14]. Another convenient selection principle is to select a path minimizing an adequate energy, which in practice is obtained by local optimization [START_REF] Kass | Snakes: Active contour models[END_REF], or by solving an eikonal equation [START_REF] Laurent | Global minimum for active contour models: A minimal path approach[END_REF]. The latter approach, chosen in this paper, has the appealing guarantee that globally optimal path is found, but it was historically more limited in the type of energies that could be minimized. The HFM library intends to alleviate this constraint, by enabling the efficient computation of globally optimal Riemannian geodesics, as well as paths minimizing a family of curvature dependent energies. We refer to [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF], for an in-depth overview of the uses of minimal paths in image processing, and devote the following paragraphs to works which specifically investigate the use of anisotropic (Riemannian or Finslerian) or non-holonomic (curvature penalizing) metrics in image processing.

Riemannian metrics were first applied in [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] to the segmentation of image regions and of tubular structures. The metric tensors are built in a preliminary step, based on a local filtering of the image, and are intended to guide the minimal paths along the image structures of interest; their design remains a challenging task, which is application dependent, requires expert knowledge, and remains the object of active research. Activity in the field [START_REF] Chen | Vessel extraction using anisotropic minimal paths and path score[END_REF][START_REF] Chen | Vessel tree extraction using radiuslifted keypoints searching scheme and anisotropic fast marching method[END_REF] was renewed after the development of fast Riemannian eikonal equation solvers [START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF]. Finslerian metrics allow path length to be measured differently depending on the direction, by e.g. augmenting a Riemannian metric with an asymmetric linear term [START_REF] Randers | On an Asymmetrical Metric in the Four-Space of General Relativity[END_REF]. They have important applications in image segmentation [START_REF] Melonakos | Finsler Active Contours[END_REF][START_REF] Chen | Finsler Geodesics Evolution Model for Region based Active Contours[END_REF]. An efficient discretization scheme exists for two dimensional finslerian eikonal equations [START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF][START_REF] Mirebeau | Anisotropic Fast Marching in ITK[END_REF], but for technical reasons2 it is not implemented in the HFM library.

Path energies featuring second order terms, such as curvature, are completely natural in image segmentation [START_REF] Kass | Snakes: Active contour models[END_REF]. Curvature penalization is for instance a good prior for region segmentation in noisy images, in particular if the objects of interest are convex [START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF]. Another use case is tubular structure segmentation in images of the retina, which feature complex overlays of vessels: curvature penalization helps eliminate "shortcuts" where the extracted path runs along the concatenation of several distinct intersecting vessels [START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF]. Finally, white fiber tractography in dMRI scans of the brain is a promising application field: the curvature penalty adds inertia to the paths, preventing them from getting lost at fiber crossings [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF]. Global optimization path energies featuring curvature and higher order tems, using dynamic programming, was first investigated in [START_REF] Strandmark | Shortest Paths with Curvature and Torsion[END_REF][START_REF] Liao | Globally Optimal Curvature-Regularized Fast Marching For Vessel Segmentation[END_REF]. We do believe that the approach used in the HFM library lies on firmer grounds, which are presented mathematically in §2.4 and which history is briefly described below. The first step was to introduce the configuration space of all positions (in the image domain) and orientations [START_REF] Pechaud | Extraction of tubular structures over an orientation domain[END_REF]. Among other things, this allows to design the cost function(s) involved in the path energy using an orientation sensitive filtering of the image, based on e.g. wavelets [START_REF] Duits | Image Analysis and Reconstruction using a Wavelet Transform Constructed from a Reducible Representation of the Euclidean Motion Group[END_REF] or Gabor filters. The second step is to relate the second order energy model of interest with a non-holonomic metric on the configuration space, see [START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF][START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF] for the Reeds-Shepp and Euler-Mumford models respectively. The final step is to design a causal discretization scheme so as to solve the resulting eikonal PDE (possibly slightly relaxed) more efficiently, which was addressed in [SBD + 15, DMMP16, Mir17a, Mir17b].

Outline. We discuss in §2 the models that can be addressed with our numerical method, and describe in detail our discretization of the related eikonal PDEs. Section §3 describes the fast marching algorithm, and the related methods of geodesic backtracking and sensitivity analysis. Numerical experiments are presented in §4.

Expression of the Hamiltonian

In this section, we describe the discretization of several generalized eikonal equations, using an original and specific representation of the corresponding Hamiltonian, as announced in the introduction (11). These implementation principles are at the foundation of our numerical method, but are invisible from the top-level interface of the HFM library. They are of interest to curious users and to those who would like to extend the core C++ library with additional models and features.

In the rest of this paper we denote by H p := (F * p ) 2 , p ∈ Ω, the squared dual metric, which is abusively referred to as the Hamiltonian. Note that, in classical mechanics, the Hamiltonian is defined as 1 2 H p , but the multiplicative factor 1/2 tends to introduce unnecessary clutter in both the description and the implementation of our numerical method. Our approach requires the Hamiltonian to be representable, exactly or approximately (in which case a consistency error is introduced), in the following general form. For any point p ∈ Ω and any co-vector p

∈ E * H p (p) ≈ max 1≤k≤K 1≤i≤I α ik p, ėik 2 + + 1≤j≤J β jk p, ḟjk 2 , ( 13
)
where a + := max{0, a} for any a ∈ R. The integers I, J, K are meta-parameters which are fixed for each class of metric. In contrast, the weights α ik , β jk ≥ 0 and the offsets ėik , ḟjk ∈ Z d implicitly depend on the current point p ∈ X of the domain, where 1 ≤ i ≤ I, 1 ≤ j ≤ J, and 1 ≤ k ≤ K. Let us mention that several strategies, discussed in §3.1.3, are used to limit the memory footprint of storing the local parameters α ik , β jk , ėik , ḟjk of the Hamiltonian. The finite-difference approximation HU (p) of a continuous function's Hamiltonian H p (du(p)), at a point p ∈ Ω, is obtained by inserting the following first order upwind approximations in the r.h.s. of ( 13)

dU (p), ė 2 + ≈ h -2 max{0, U (p) -U (p -h ė)} 2 , ( 14
) dU (p), ė 2 ≈ h -2 max{0, U (p) -U (p -h ė), U (p) -U (p + h ė)} 2 . ( 15
)
We denoted by U : X ∪ ∂X → R a discrete map, where X and ∂X are disjoint finite subsets of the cartesian grid hZ d of scale h > 0, devoted to approximating Ω and ∂Ω. By convention, if needed, U is extended by +∞ outside of its domain so as to implement outflow boundary conditions. The eikonal equation (3) is thus discretized in the form

∀p ∈ X, HU (p) = 1, ∀p ∈ ∂X, U (p) = σ(p). (16)
To be more specific, our implementation requires the discretization domain X ∪ ∂X to be box shaped. It also allows for various types of boundary conditions, devoted to e.g. the commonly used manifolds R 2 × S 1 and R 3 × S 2 , and for axis-dependent grid scales, such as a physical scale and an angular scale in the previous case. See §D for more detail on the discretization grid conventions.

Before turning to specific models, we formally define the properties of monotony and causality, already mentioned in the introduction §1.3. These properties enable solving the system (16) in a single pass using the fast marching algorithm, see §3.1. Definition 2.1. A numerical scheme on a finite set Z is a map

H : Z × R × R Z → R. It is said
• Monotone iff H is non-decreasing w.r.t. the second and (each of the) third variables.

• Causal iff H only depends on the positive part of the third variable(s).

Slightly abusively one denotes for any U : Z → R and any p ∈ Z HU (p) := H p, U (p), (U (p) -U (q)) q∈Z .

For notations to match between Definition 2.1 and (16), one must set Z := X ∪ ∂X, and define HU (p) := U (p) for all p ∈ ∂X. By design, these properties are obeyed by any numerical scheme obtained by inserting the upwind finite differences (14) in an expression of the form (13), namely

HU (p) := h -2 max 1≤k≤K 1≤i≤I α ik max 0, U (p) -U (p -h ėik ) 2 (17)
+ 1≤j≤J β jk max 0, U (p) -U (p -h ḟjk ), U (p) -U (p + h ḟjk ) 2 .

Dijkstra's algorithm

We show that, with a special choice of parameters, our discretization scheme (16) reduces to the classical problem of finding shortest paths on a graph with non-negative edge lengths. In addition, the fast marching algorithm reduces to Dijkstra's method in that case. Indeed, assume that I = 1 and J = 0, while K ≥ 1 remains arbitrary. Then HU (p) = 1 rewrites as

max 1≤k≤K U (p) -U (p -h ėk ) l k (p) 2 + = 1, equivalently U (p) = min 1≤k≤K l k (p) + U (p -h ėk ).
We denoted l k := 1/ √ α 1k and ėk := ė1k , which both (sometimes implicitly) depend on the current point p. Almost any optimal control problem can be addressed with Dijkstra's algorithm, by using sufficiently wide stencils, including those involving curvature or torsion dependent energies [START_REF] Strandmark | Shortest Paths with Curvature and Torsion[END_REF]. We advocate however for more clever designs of the weights and offsets, as discussed below for several metric structures of interest, so as to improve the accuracy of the numerical results.

Isotropic metrics

We recall in this subsection the Eulerian discretization of isotropic eikonal equations, which dates back to [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF]. As discussed in the general introduction, isotropic metrics are (by definition) locally proportional to the euclidean norm • on the ambient space E := R d : one has F p ( ṗ) = c(p) ṗ , for any point p ∈ Ω and any vector ṗ ∈ E, where c : Ω →]0, ∞[ is a cost function provided by the user. The dual metric reads F * p (p) = c(p) -1 p , for any co-vector p ∈ E * . The Hamiltonian is representable exactly in the form (13), using only d scalar products (I = 0, J = d, K = 1). Denoting by ėi = (0, • • • , 0, 1, 0, • • • , 0) the unit vector directed along the i-th coordinate axis, for any 1 ≤ i ≤ d, one has indeed

H p (p) = c(p) -2 p 2 = c(p) -2 1≤i≤d p, ėi 2 . ( 18
)
Note that, in this special case, the offsets ( ėi ) d i=1 do not depend on the current base point p ∈ Ω. For completeness, we write the full form of the discretization scheme, specializing (17) and recovering [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF]. For any U : X ∪ ∂X → R and any grid point p ∈ X

HU (p) = (c(p)h) -2 1≤i≤d max{0, U (p) -U (p -h ėi ), U (p) -U (p + h ėi )} 2 . ( 19
)
Several numerical schemes often exist for discretizing a given eikonal equation. Substituting

p, ėi 2 = p, ėi 2 -+ p, ėi 2
+ in (18) yields for instance the alternative discretization

HU (p) = (c(p)h) -2 1≤i≤d τ ∈{-1,1} max{0, U (p) -U (p -τ h ėi )} 2 . ( 20
)
This scheme is typically less accurate than (19) near the cut locus (the points of non-differentiability of the distance map), but has the advantage of being differentiable w.r.t. the entries of the discrete map U , which is welcome if the eikonal PDE solution is part of a larger optimization problem, see §4.2. Diagonal metrics are an elementary generalization of isotropic metrics, in which the propagation cost is possibly distinct along each coordinate axis. The metric, and dual metric, read:

F p ( ṗ) = 1≤i≤d c i (p) 2 ėi , ṗ 2 , F * p (p) = 1≤i≤d c i (p) -2 p, ėi 2 . ( 21
)

Riemannian metrics, and sub-Riemannian approximations

Riemannian metrics are the most common class of non-isotropic metrics. They are determined by a field M : Ω → S ++ (E) of positive definite tensors, and take the form F p ( ṗ) := ṗ M(p) , where we recall ṗ M := M ṗ, ṗ . In differential geometry, the natural distance on a submanifold of R n can be described by a Riemannian metric on a parametrization domain. In image segmentation, the Riemannian tensors may stem naturally from the data [FDHTV + 16], or be creatively designed based on some local image analysis [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF].

In order to discretize Riemannian eikonal equations, we introduce an adequate decomposition of positive definite tensors. The condition number of a positive definite tensor

D ∈ S ++ (E d ) is defined by Cond(D) := D D -1 . ( 22
) Proposition 2.2. Let D ∈ S ++ (E)
, where E := R d , and let d := d(d + 1)/2. Then there exists non-negative weights α i ≥ 0 and integer offsets ėi ∈ Z d , where

1 ≤ i ≤ d , such that D = 1≤i≤d α i ėi ⊗ ėi . (23)
Furthermore, this decomposition can be chosen such that ėi ≤ Cond(D) α , for all 1 ≤ i ≤ d , where α := d -1. In dimension d ≤ 3, one has the improved estimate α = 1. We assume in the following that the decomposition (23) is the one provided by Voronoi's first reduction of the quadratic form ṗ → ṗ, D ṗ , which enjoys these properties.

We refer to [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] for the proof of Proposition 2.2, and to [START_REF] Schürmann | Computational geometry of positive definite quadratic forms[END_REF] for more background on Voronoi's theory and the field of additive lattice geometry. Our numerical codes devoted to Riemannian metrics do compute the tensor decomposition (23) at each gridpoint, using a simple and efficient algorithm due to Selling [START_REF] Selling | Ueber die binären und ternären quadratischen Formen[END_REF], see Appendix B. This algorithm applies in dimension d ∈ {2, 3} only, but we hope to address the case d ∈ {4, 5} in future using other techniques. Our discretization scheme directly involves the offsets (e i ) d i=1 , hence their norm should be as small as possible. The chosen construction is proved in [START_REF] Mirebeau | Minimal stencils for discretizations of anisotropic PDEs preserving causality or the maximum principle[END_REF] to be optimal in this regard, in dimension d = 2, and sharp worst and average case estimates of the stencil radius are also established. The tensor decomposition (23) can be used to discretize PDEs other than eikonal equations, such as anisotropic diffusion [START_REF] Fehrenbach | Sparse non-negative stencils for anisotropic diffusion[END_REF]. Proposition 2.2 yields an exact representation of Riemannian Hamiltonians in the desired form (13), using d := d(d + 1)/2 scalar products (I = 0, J = d , K = 1). Indeed, for any point p ∈ Ω and any co-vector p ∈ E * , one has denoting by D(p) := M(p) -1 the inverse of the Riemannian metric tensor

H p (p) = p 2 D(p) = p, D(p)p = 1≤i≤d α i (p) p, ėi 2 . ( 24
)
The weights α i (p) ≥ 0 and the offsets ėi = ėi (p) ∈ Z d are those appearing in the decomposition (23) of D(p). Specializing (17) yields

HU (p) = 1≤i≤d α i (p) max{0, U (p) -U (p -h ėi ), U (p) -U (p + h ėi )} 2 . ( 25
)
In comparison with the classical discretization (19) of isotropic eikonal equations, the offsets are point dependent: e i = e i (p), and are slightly more numerous: d instead of d. Interestingly, if the Riemannian metric tensor is proportional to the identity matrix, M(p) = c(p) 2 Id at some point p ∈ Ω, then one can show that (24) reduces to the previous numerical scheme (19) (in particular d -d of the coefficients α i (p) vanish, while the others are equal to c(p) -2 ). Sub-Riemannian metrics can be regarded as degenerate Riemannian metrics, which metric tensors have some infinite eigenvalues. The inverse tensor D : Ω → S + (E * ) is well defined, positive semi-definite, but rank deficient. They define non-holonomic control models, in the sense that, at some points in space p ∈ Ω, some directions of motion ṗ ∈ E are forbidden, namely those outside of the span of D(p). Proper sub-Riemannian metrics should also obey a local controllability property, expressed in terms of commutators of vector fields, which is out of the scope of this paper, see [START_REF] Montgomery | A Tour of Subriemannian Geometries, Their Geodesics and Applications[END_REF]. The HFM library addresses sub-Riemannian metrics using a relaxation approach, involving a family F ε of Riemannian metrics which tensors M ε explode as some parameter ε → 0. This technique is viable numerically thanks to the good behavior of our numerical scheme with strongly anisotropic metrics. In practice, tensors of condition number Cond(M ε ) ≈ 10 offer a good compromise between (a) the sharpness of the sub-Riemannian relaxation, and (b) the size of the discretization stencils, see Proposition 2.2.

The Reeds-Shepp car model (with reverse gear) [START_REF] Reeds | Optimal paths for a car that goes both forwards and backwards[END_REF] is perhaps the most notorious example of a sub-Riemannian metric, also appearing in Petitot's description of the visual cortex V1 [START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF]. This model is posed on the manifold M := R 2 × S 1 of positions and orientations, and is described here for a positive relaxation parameter ε > 0. For any point (x, θ) ∈ M and any tangent vector ( ẋ, θ)

∈ T (x,θ) M = R 2 × R we define following [SBD + 15] F ε (x,θ) ( ẋ, θ) 2 = ( ẋ, θ) 2 Mε(x,θ) := c(x, θ) 2 n(θ), ẋ 2 + ε -2 n(θ) ⊥ , ẋ 2 + (ξ θ) 2 , ( 26
)
where c : Ω →]0, ∞[ is a given cost function, and ξ > 0 is a parameter homogeneous to a radius of curvature. We denoted by n(θ) := (cos θ, sin θ) the unit vector of orientation θ ∈ S 1 , and by n(θ) ⊥ the clockwise orthogonal vector. Lateral physical motions, in the direction of n(θ) ⊥ , see their cost strongly penalized by the ε -2 factor in (26), and are completely excluded in the limit sub-Riemannian model for which ε = 0. Our software also implements the higher dimensional generalization of this model on R 3 ×S 2 , first considered in [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF], as well as a "dual" variant introduced in [Mir17a] which penalizes non-planarity (akin to a torsion penalization) instead of curvature.

We can make explicit the tensor field

M ε : R 2 × S 1 → S ++ (R 3 ) such that F ε (x,θ) ( ẋ, θ) = ( ẋ, θ) Mε(x,θ) , for all (x, θ) ∈ R 2 × S 1 . The matrix M ε (x, θ) is block-diagonal with structure
(2 × 2, 1 × 1), and reads

M ε (x, θ) = c(x, θ) 2   n(θ) ⊗ n(θ) + ε -2 n(θ) ⊥ ⊗ n(θ) ⊥ ξ 2   .
The eigenvectors are (n(θ), 0), (n(θ) ⊥ , 0), and (0, 1), with corresponding eigenvalues c(x, θ) 2 (1, ε -2 , ξ 2 ). Computing the inverse tensors D ε (x, θ) := M ε (x, θ) -1 , we obtain the following expression of the Hamiltonian: for any co-vector (x, θ)

∈ T * (x,θ) M H ε (x,θ) (x, θ) = (x, θ) 2 Dε(x,θ) := c(x, θ) -2 x, n(θ) 2 + ε 2 x, n(θ) ⊥ 2 + ( θ/ξ) 2 . ( 27
)
From this point, the HFM library applies the previous Riemannian discretization strategy (24), with a positive relaxation parameter. In practice, choosing ε := 0.1 yields good results. Note that the matrix D ε (x,θ) of the sub-Riemannian Hamiltonian is rank deficient when ε = 0, hence Proposition 2.2 is not directly applicable to the limit model.

Curvature penalized models

This section is devoted to the computation of paths globally minimizing a curvature dependent energy, following [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF], which is one of the main novelties brought by the HFM library. For that purpose we reformulate these second order planar paths models as ordinary shortest path problems in the three dimensional domain M := R 2 × S 1 , and w.r.t. a singular metric F : T M → [0, ∞], see §1.2. Said otherwise, the introduction of a second order path derivative in the cost function comes at the price of (i) an extra dimension S 1 , and (ii) a strongly (formally infinitely) anisotropic metric F. This approach shares many similarities with the sub-Riemannian Reeds-Shepp model (26), the main difference being that the models considered in this section lack the ability to shift into reverse gear. We describe below the computation of the dual metric, and its approximation in the form (13) which is required for our discretization strategy.

The dual metric F * to the models of interest is described in (8) as the solution to a one dimensional optimization problem. More precisely, one has

F * (x,θ) (x, θ) = c(x, θ) -1 f (x, θ),
where x := x, n(θ) ∈ R, and where f : R 2 → [0, ∞[ is defined as the following maximum

f (x, θ) := sup θ∈R x + θ θ C( θ) , extremal when θC( θ) = (x + θ θ)C ( θ). ( 28
)
We specialize, in the next equation, the optimality condition (28, right) for the curvature costs C RS , C EM and C D associated to the Reeds-Shepp, Euler-Mumford, and (non-smooth) Dubins models, see (9). After simplification this condition respectively reads

θxξ 2 = θ, θ2 ξ 2 θ + 2 θξ 2 x -θ = 0, θξ ∈ {-1, 1, ∞}.
Recall that ξ > 0 is a parameter, homogenous to a radius of curvature, which modulates the intensity of curvature penalization. Solving for θ, and inserting the optimal value in (28, left), one obtains the following expression of f (x, θ) respectively We refer to Appendix A for the details of the computations, and to [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] for a different proof.

x2 + + ( θ/ξ) 2 , x + x2 + ( θ/ξ) 2 , max{0, x + θ/ξ, x -θ/ξ} (29)
Let us focus on the case (29, left) of the Reeds-Shepp forward model, which squared dual metric thus reads

H RS (x,θ) (x, θ) = c(x, θ) -2 x, n(θ) 2 + + ( θ/ξ) 2 . ( 30
)
The main difference with the classical Reeds-Shepp reversible model (with reverse gear), which Hamiltonian is obtained by setting ε = 0 in (27), is the positive part in the physical term x, n(θ) 2 + . This difference accounts for the lack of reverse gear in the present model. The discretization of this term is addressed in the next proposition, proved in [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF], which introduces (similarly to the reversible case) a relaxation parameter ε > 0. We denote by P n the orthogonal projection onto the hyperplane orthogonal to a unit vector n:

P n := Id -n ⊗ n.
Proposition 2.3. Let n ∈ E := R d and let ε > 0. Consider non-negative weights and offsets (α i , ėi ) I i=1 ∈ (R + × E) I , obtained e.g. by Proposition 2.2, such that for all x ∈ E * one has

x, n 2 + ε 2 P n x 2 = 1≤i≤I α i x, ėi 2 .
Assume that n, ėi ≥ 0 (otherwise replace ėi with its opposite), for all 1 ≤ i ≤ I. Then

x, n 2 + ≤ 1≤i≤I α i x, ėi 2 + ≤ x, n 2 + + ε 2 P n x 2 .
Denote by α ε i (θ) and e ε i (θ) the weights and offsets obtained by applying Proposition 2.3 to the vector n(θ) ∈ R 2 and a suitably small ε > 0, where 1 ≤ i ≤ I and I = d = d(d + 1)/2 = (2 × 3)/2 = 3, see Proposition 2.2. Then we approximate the Hamiltonian (30) of the Reeds-Shepp forward model in the form (13) with parameters I = 3, J = 1, K = 1, as follows

H RS (x,θ) (x, θ) ≈ c(x, θ) -2 1≤i≤I α ε i (θ) x, ėε i (θ) 2 + + ( θ/ξ) 2 .
For completeness, we write the full discretization scheme for this model, specializing (17). Let h > 0 be the gridscale and let M h := hZ 2 × (hZ/2πZ). We assume here that 2π/h is an integer3 , so that the discretization of the periodic circle S 1 makes sense. For any U :

M h → R and any (x, θ) ∈ M h H RS ε U (x, θ) :=(c(x, θ)h) -2 1≤i≤I α ε i (θ) max{0, U (x, θ) -U (x -h ėε i (θ), θ)} 2 + ξ -2 max{0, U (x, θ) -U (x, θ -h), U (x, θ) -U (x, θ -h)} 2 .
The Dubins Hamiltonian is discretized using a similar procedure. In view of (29, right) one has indeed

H D (x,θ) (x, θ) = c(x, θ) -2 max{( x, n(θ) + θ/ξ) 2 + , ( x, n(θ) -θ/ξ) 2 + }. (31)
Apply Proposition 2.3 to the three dimensional vectors v + (θ) := (n(θ), ξ -1 ) and v -(θ) := (n(θ), -ξ -1 ), and denote by α ε +i (θ), α ε -i (θ), ėε +i (θ) and ėε -i (θ) the corresponding weights and offsets, where 1 ≤ i ≤ I and I = (3 × 4)/2 = 6. Then we approximate (31) in the form (13) with parameters I = 6, J = 0, K = 2, as follows

H D (x,θ) (x, θ) ≈ c(x, θ) -2 max 1≤i≤I α ε +i (θ) p, ėε +i (θ) 2 + , 1≤i≤I α ε -i (θ) p, ėε -i (θ) 2 + .
The final PDE discretization scheme is as usual obtained by specializing (17). The Euler-Mumford elastica model is addressed using a similar procedure, for which we refer to [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] due to space constraints.

Implementation

In this section, we provide the implementation details for the Fast-Marching algorithm, used to numerically solve our Eulerian and causal discretization (16) of the eikonal PDE. We strongly rely on the specific form (13) of the Hamiltonian approximation. We also discuss geodesic backtracking, and sensitivity analysis.

In order to keep notations simple, we assume that the Hamiltonian H of the model of interest is representable in the form (13) with parameters J = 0 and K = 1, whereas I remains arbitrary. The following arguments can be easily adapted to the general case where J and K are arbitrary as well. For any discretization point p ∈ X, any co-vector p ∈ E * , and any discrete map U , the Hamiltonian and finite differences scheme are thus assumed to read

H p (p) = 1≤i≤I α i max{0, p, ėi } 2 , HU (p) := h -2 1≤i≤I α i (U (p) -U (p -h ėi )) 2 + . ( 32
)
The weights α i = α i (p) ≥ 0 and offsets ėi = ėi (p) ∈ Z d , where 1 ≤ i ≤ I, depend on the current point p ∈ X. Recall that the objective is to find U :

X ∪ ∂X →] -∞, ∞] obeying ∀p ∈ X, HU (p) = 1, ∀p ∈ ∂X, U (p) = σ(p), (33)
where X and ∂X are finite subsets of hZ d , and where the boundary data σ : ∂X →] -∞, ∞] is given. In the following, we denote by N := #(X ∂X) the cardinality of the discrete domain.

Fast marching

The fast marching algorithm is a generalization of Dijkstra's algorithm [START_REF] Wybe | A short introduction to the art of programming[END_REF] for the computation of shortest paths on graphs. Said otherwise, it is a specialization of the dynamic programming principle. The general structure is unchanged, see Algorithm 1, and in particular each point of the domain X ∪ ∂X is tagged as Accepted only once, after what the corresponding value of the computed solution U is frozen.

Algorithm 1 Fast marching Pre-compute the reversed stencils (V [p]) p∈X∪∂X , defined in (35). Tag points of ∂X as Trial, and other points as Far.

Initialize U : X ∪ ∂X →] -∞, ∞] to the value +∞ on X, and the boundary condition on ∂X.

While there remains Trial points. Find a Trial point q minimizing U . Tag q as Accepted, and call Post-process(q). For each neighbor p ∈ V [q], either Far or Trial If p is Far, then tag it as Trial, and call Pre-process(p). Update U (p), taking into account the value U (q) of the last Accepted point.

The collection of points tagged as Trial can be regarded as a propagation front. Two additional sub-routines, appearing as Pre-process(p) and Post-process(q) in Algorithm 1, and referred to as pre-processing and post-processing, are called when a point enters and leaves the front respectively. By default, some memory management is performed at these points, see §3.1.2 and §3.1.3, but additional tasks may optionally be plugged in as well. At the post-processing stage one may for instance check for stopping criteria, or slightly alter the value U (p) before it is frozen using the (formally) second order HAFMM numerical scheme [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF].

Stencils and reversed stencils

Our algorithm uses point dependent, adaptive stencils, hence we need to cautiously determine the dependency graph underlying the discretized problem (33). For that purpose, we introduce the forward stencil V(p) ⊂ X ∪ ∂X, of an arbitrary interior point p ∈ X, and its counterpart the reversed stencil V [q] ⊂ X, for all q ∈ X ∪ ∂X. The forward stencil V(p) at p ∈ X, collects all points q ∈ X ∪ ∂X whose value U (q) is involved in the definition HU (p), see (32, right).

V(p) := {q = p -h ėi (p); 1 ≤ i ≤ I, q is visible from p} ∩ (X ∪ ∂X), (34)
The intersection with X ∪ ∂X is meant to exclude all points q which fall outside of the (box shaped) discrete domain. Outside of this domain, the unknown function U is by convention extended by +∞, which amounts to implement outflow boundary conditions. The visibility constraint is a test ensuring that the segment [p, q] does not intersect any of the obstacles that may optionally be introduced within the domain. This condition is required because our adaptive stencils are often 5 to 10 pixels wide, hence they would otherwise jump over thin obstacles, such as walls in in building maps which are usually drawn 1 pixel wide. See Figure 6. The reversed stencil V [q] ⊆ X at q ∈ X ∪ ∂X, is defined by inverting the stencil connectivity graph:

V [q] := {p ∈ X; q ∈ V(p)}. (35)
Thus V [q] collects all points p ∈ X such that HU (p) depends on U (q). Forward and reverse stencils coincide in the classical discretization of isotropic fast marching (2.2), but differ in most other cases.

Valid obstacles

Incorrect obstacles

Figure 6: Obstacles are specified to the HFM software as boolean images. They must be "water tight": contiguous obstacle voxels must share a full (d -1)-dimensional face. Otherwise the numerical scheme will leak through, and the extracted minimal paths may jump across walls.

The numerical cost of constructing the stencils and the reversed stencils can be neglected in models for which these sets are independent of the current point (e.g. isotropic metrics §2.2), or depend only on a strict subset of its coordinates (e.g. the angular coordinate for curvature penalized models §2.4), by exploiting redundancy as discussed in §3.1.3 below. In other cases (e.g. Riemannian metrics §2.3), the construction of the direct stencils V(p) has complexity4 O(I) for each p ∈ X, thus O(N I) overall. Once this is done, the reversed stencils are computed by sorting the pairs {(p, q); p ∈ X, q ∈ V(p)} according to their second element, which costs O(IN ln N ).

Finally, let us mention that equations (34) and ( 35) do not reflect our implementation entirely faithfully, for two reasons. First, we actually store the offsets instead of the points, e.g. ėi (p) ∈ Z d instead of p -h ėi (p) in (34). Second, the visibility and out of domain tests are performed while dynamic programming is running, and not during the initial stencil construction as suggested in (34).

Elementary update

Consider a point p ∈ X tagged Trial, and which value U (p) must be updated w.r.t. the Accepted points, as specified in the last line of Algorithm 1. This means that the currently stored value U (p) must be replaced with the largest solution λ ∈ R to the following univariate quadratic equation:

i∈I α i (p) (λ -U (q i )) 2 = 1,
where

I := {1 ≤ i ≤ I; q i ∈ V(p), q i is Accepted}. ( 36
)
The neighbor points q i := p -h ėi (p), where 1 ≤ i ≤ I, are those appearing in the numerical scheme HU (p) at p, see (32, right). Note that the index set I excludes points that are out of the domain X ∪ ∂X, are not visible from p, or are tagged Trial or Far. One can show that, by construction of the fast marching algorithm 5 , equation (36) has two real roots, the smallest of which is a numerical artifact, while the largest denoted λ * satisfies by construction λ * ≥ U (q i ) for all i ∈ I, and defines the updated value of U (p). Some of our discretization schemes involve a rather large number of neighbors, for instance I = 12 for three dimensional Riemannian metrics, and I = 27 for the (default implementation of the) Euler-Mumford elastica model. The latter figure may increase if the user requests more accuracy, see the discretization details in [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF]. In order to maintain a O(1) complexity of the elementary update, we cache for each Trial point p ∈ X the coefficients (a, b, c) of the second degree polynomial aλ 2 -2bλ + c defined by (36). This cache is initialized as (0, 0, -1) when Pre-process(p) is called, is increased (component-wise) by α i (p) (1, U (q i ), U (q i ) 2 ) each time a new neighbor q i is Accepted, and is deleted when Post-process(p) is called.

The fast marching algorithm, similarly to Dijkstra's algorithm, must also maintain a queue of all Trial points, sorted by increasing values of U . For each point p ∈ X ∂X, queue maintenance operations include the insertion of the new key (when p is first tagged as Trial), a local re-ordering when the attached value U (p) is modified (at most I times), and eventually the key removal (when p is tagged as Accepted). The overall complexity, using a heap based sorted data structure, is thus O(IN ln N ).

Memory usage

The memory footprint of our algorithm is linear w.r.t. the number N of discretization points.

The proportionality constant does matter however, since the image resolution of e.g. dMRI medical data is often huge. In addition, several applications introduce extra domain dimensions, accounting for e.g. the orientation, grayscale, or radius of the extracted tubular structure [CMC16a, PKP09], which multiplies the number of discretization points. For this reason, the HFM software implements several memory optimizations, discussed in this paragraph. We use the O(•) notation6 to denote space complexity, counted in Bytes.

Our generalized fast marching solver uses wide stencils, which weights and offsets take significant memory space if stored indiscriminately. By default7 we use weights α i ∈ R of "double" type, hence occupying 8 Bytes8 , and offsets ėi ∈ Z d with "signed char" components, hence occupying d Bytes, where d is the domain dimension, and where 1 ≤ i ≤ I. The storage cost of the forward and the reversed stencils is thus respectively

O(IN 1 (8 + d)) O(IN 2 d)
where N 1 and N 2 are the number of actually stored stencils. An extensive storage (N 1 = N 2 = N ) of the numerical scheme stencils thus costs O(I(8 + 2d)N ), whereas the memory footprint of the input speed function and of the output value function is O(2 × 8N ). For the Euler-Elastica model, which uses quite large stencils (I = 27 and d = 3), the space complexity O(378N ) of our numerical method would thus vastly exceed the incompressible external cost O(16N ) of input and output. Fortunately, we have developed two data management strategies for stencil storage which make the former cost negligible for this specific model, and severely limit it for others.

• Data sharing. The stencil offsets ėi (p) of some models only depend on a strict subset p of the coordinates of the current point p ∈ X, e.g. none (p = ∅) for isotropic metrics §2.2, or the angular coordinate only (p = θ where p = (x, θ)) for curvature penalized models §2.4 such as the Euler-Mumford elasticae. In addition, the weights are proportional to (the inverse square of) some user provided cost functions c(p, l), p ∈ Ω, 1 ≤ l ≤ L, where L is one or is a small positive integer. In summary, the discretization stencil structure takes the form ėi = ėi (p),

α i = c(p, l i (p)) -2 α i (p). ( 37
)
We take advantage of this structure by storing only the offsets ėi (p) ∈ Z d , weights α i (p) ≥ 0 and cost index l i (p) ∈ 1, L (omitted if L = 1) associated with the points of the lower dimensional domain p ∈ X. We denoted 1,

L := [1, L] N. Thus N 1 = N 2 = #(X) N .
For instance, in the Euler-Mumford elastica model, N 1 = N 2 is the number of angular directions, which is typically around 60. Thus the huge footprint O(378N ) of extensive stencil storage for this model, mentioned above, is with our strategy reduced to a negligible constant cost O(378 × 60).

• Recomputation. For a number of models, either Riemannian or involving generalized curvature costs (10), the numerical scheme stencils lack redundancy, hence the previous data sharing strategy is not applicable. We thus implement an alternative strategy for memory usage reduction, which is to only store the direct stencils V(p) of the points currently tagged as Trial. Stencils must then be computed twice: (i) initially, so as to construct the reversed stencils V [p], and (ii) when a point label changes from Far to Trial. The storage cost of the forward stencils becomes negligible, since

N 1 ≈ N 1-1/d N .
The additional computational complexity introduced by these recomputations is moderate -actually, it is dominated by the incompressible cost of maintaining the priority queue of all Trial points. The reversed stencils V [p] still must be stored all, thus N 2 = N , but fortunately they use only a fraction of the memory space of the forward ones (namely d/(8 + d) where typically d ∈ {2, 3}).

For the Euler-Mumford elastica model (I = 27 and d = 3), generalized in the sense of (10), the cost of extensive stencil storage O(378N ) is reduced with this recomputation strategy to O(81N + 297N 2 3 ). In comparison, the cumulated cost of the input parameter fields ξ, κ, c : X → R, and of the output, is O(32N ).

Our algorithm also stores the coefficients of the second degree polynomials (36) associated with the Trial points, in addition to the stencils, which incurs the memory cost O(4N +24N 1-1/d ) with the chosen implementation9 . Finally, we chose to store the set I of active neighbors of each point, defined in (36), which costs O( ln 2 I N ), for use in geodesic backtracking and sensitivity analysis.

Geodesic extraction

Geodesic backtracking is done by solving an ODE (5) involving the distance map, which is numerically approximated in a preliminary step by solving the discretized eikonal PDE (33). Despite the apparent simplicity of numerical ODE integration, this task deserves caution. Naive implementations of ODE geodesic backtracking indeed suffer from artifacts such as paths (i) interrupted nearby obstacles before they reach the seed points, (ii) going past the seed points, or (iii) endlessly oscillating close to singularities of the distance map. These difficulties are particularly marked when using metrics which are (a) strongly inhomogeneous, (b) strongly anisotropic, leading to wide discretization stencils, and (c) non-holonomic and/or non-locally controllable, often leading to discontinuous (in addition to being non-differentiable) front arrival times.

The proposed software implements two robust numerical methods for geodesic backtracking. They produce similar results in practice, but are based on distinct principles that we discuss below. These methods presently lack a proper convergence analysis, in particular for nonholonomic models, hence they may be regarded as partly heuristic.

Modified Euler method using upwind gradients

Our first approach to geodesic backtracking is to solve the ODE (5) using an explicit second order integration scheme, and an upwind approximation of the geodesic flow direction. More precisely, consider the differential equation

γ(t) = V (γ(t)),
where

V (p) := dF * p (du(p)), (38)
for all p ∈ Ω. We use the integration scheme10 γ(t + δ) ≈ γ(t) + δV γ(t) + 1 2 δV (γ(t)) , often called the modified Euler method, the midpoint rule, or the Runge-Kutta-2 method. The timestep δ > 0 is locally adjusted so that |δV (γ(t))| equals a fraction of the gridscale, by default h/4. The stopping criterion involves, for robustness, a safety radius around the seed points.

The main difficulty for implementing this method lies in the (approximate) evaluation of the geodesic flow direction V (p). For instance, naively approximating the gradient du(p) using centered finite differences often yields unstable results. Instead, we use an upwind approximation of V (p), derived from the numerical scheme, namely

V (p) ≈ V (p) := h -1 1≤i≤I α i (U (p) -U (p -h ėi )) + ėi . (39)
Let us formally justify this expression. On the one hand one has

dH p (p) = 2 1≤i≤I α i p, ėi + ėi , d p H(du(p)) ≈ 2h -1 1≤i≤I α i (u(p) -u(p -h ėi )) + ėi .
Substituting the PDE solution u with its numerical approximation U , one obtains (39) as announced. It is well known that the geodesics on a manifold move at constant speed, a property which in our context reads F p (V (p)) = 1 (under some differentiability assumptions), and can be derived from (38, right) and Legendre-Fenchel duality. The following proposition, never published elsewhere, further motivates (39) by establishing a discrete counterpart of this property. It is stated in more formal and abstract terms, and proved, in Appendix §C.

Proposition 3.1. Under the assumptions of our discretization, one has F p ( V (p)) ≤ 1 for any p ∈ X. In addition, equality holds if U coincides with an affine function over the stencil {p} ∪ {p -h ėi } 1≤i≤I .

The HFM library sets V (p) = 0 on boundary points p ∈ ∂X, and extends V to the continuous domain Ω ⊃ X ∪ ∂X by bilinear interpolation. Outlier points, inside walls or attached to inconsistently large11 values of U , are for robustness excluded from the interpolation.

Diffuse geodesics and reverse algorithmic differentiation.

Our second backtracking method is based on the following principle. Consider a smooth field ζ : Ω → R, that should be regarded as a perturbation of the front speed. Define u ε : Ω →]-∞, ∞] for any sufficiently small ε ∈ R by

∀p ∈ Ω, F * p (du ε (p)) = 1 + εξ(p), ∀p ∈ ∂Ω, u ε (p) = σ(p).
Note that u 0 = u is the solution to the original unperturbed shorted path problem. Under suitable assumptions, the front arrival time u ε (p * ) at a point of interest p * ∈ Ω can be differentiated w.r.t. the parameter ε, and denoting by γ * the minimal geodesic associated with u 0 (p * ) one has at first order

u ε (p * ) = u 0 (p * ) + ε 1 0 F γ * (t) ( γ * (t)) ξ(γ * (t)) dt + o(ε).
This derivation follows from the envelope theorem, see [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF]. Our second geodesic backtracking method takes advantage of a discrete counterpart of this property. Consider the discrete solution U ε :

X ∂X →] -∞, ∞] to ∀p ∈ X, HU ε (p) = (1 + εξ(p)) 2 , ∀p ∈ ∂X, U ε (p) = σ(p).
Generically, U ε is differentiable w.r.t. the parameter ε, and one has at first order

U ε (p * ) = U (p * ) + ε q∈X ρ(q)ξ(q) + o(ε),
for some coefficients ρ : X → R. The assumption underlying our second backtracking method is that the bulk of the coefficients ρ is supported along the minimal geodesic γ * of the continuous problem. We use sensitivity analysis by reverse accumulation to compute ρ, see §3.3. This diffuse support is summarized into a unique path using an averaging procedure. However, path extraction is restarted if the support of ρ is spread over an excessively large domain. This happens if e.g. there is no uniqueness of the minimal geodesic, in which case the support of ρ is split into several components, which is detected using an ad-hoc criterion, see Figure 14 (right). The computation of the weights ρ is described in the following subsection.

Sensitivity analysis

Sensitivity analysis is a generic name for (semi-)automatic (meta-)programming techniques aimed at computing the first order differential of a function defined algorithmically as the composition of several elementary functions. This approach combines the formal differentiation of the elementary functions, with the numerical propagation of the so-called sensitivities, using the composition rules for derivatives. Denoting by F : R m → R n and p ∈ R m the given map and point of interest, the objective is to evaluate the Jacobian matrix dF p of size m × n. Interestingly, sensitivity analysis comes in two flavors, which time and space complexities differ.

They are referred to as forward and reverse, and are mostly adequate when m n and when m n respectively. The HFM library implements sensitivity analysis "by hand", since this is particularly simple and efficient for the addressed problem, instead of relying on meta-programming automatic differentiation techniques [START_REF] Griewank | Evaluating Derivatives. Principles and Techniques of Algorithmic Differentiation[END_REF]. Our approach allows enhancements in performance and convenience, but also has a few limitations: differentiation is only possible w.r.t. certain parameters, and some variants of the numerical scheme are unsupported (i.e. slightly inexact results will be obtained with the second order accurate method, or using time-dependent parameter fields).

To our knowledge, previous literature only considered forward differentiation of isotropic fast marching [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF]. Our contribution is thus two fold: we generalize these previous works to anisotropic fast marching, and we implement reverse sensitivity analysis for the first time in this context. The algorithms presented in this subsection are used in [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF] to solve two player zero-sum games, where the first player places a surveillance system and the second player wants to visit a target undetected.

Inputs and outputs of the differentiation methods

The Hamiltonian appearing in the eikonal equation is represented internally in the HFM software in a sum of squares form (32), involving weights and offsets, denoted by α i (p) ≥ 0 and ėi (p) ∈ Z d , where p ∈ X and 1 ≤ i ≤ I. Our differentiation techniques assume that the offsets remain constant, but that the weights are proportional to (the inverse square of) one or several12 user provided cost functions c(p, l), 1 ≤ l ≤ L, which themselves are subject to a linear perturbation εξ(p, l). See the discussion (37) on the shared stencil data structure. We also assume a linear perturbation εζ(q) of the boundary conditions σ(q), q ∈ ∂X. Summarizing, the perturbed cost functions and boundary conditions read c ε (p, l) := c(p, l) + εξ(p, l), σ ε (q) := σ(q) + εζ(q),

for all p ∈ X, all q ∈ ∂X, and all 1 ≤ l ≤ L. Following (37), the weights of the perturbed Hamiltonian read α ε i (p) = α i (p)c ε (p, l i (p)) -2 . We assume that the discretized problem solution U ε is indeed differentiable at the parameter ε = 0. Forward mode. The procedure input consists of two fields ξ : X× 1, L → R and ζ : ∂X → R. The output is the first order derivative µ = U 0 , which obeys for all p ∈ X ∪ ∂X

U ε (p) = U 0 (p) + εµ(p) + o(ε).
Reverse mode. The procedure input is a point p * ∈ X ∪ ∂X, and the output consists of two families of weights ρ : X × 1, L → R and π : ∂X → R such that one has at first order, for any perturbation fields ξ and ζ

U ε (p * ) = U 0 (p * ) + ε p∈X, l∈ 1,L ρ(p, l)ξ(p, l) + q∈∂X π(q)ζ(q) + o(ε).
One may also request the Taylor expansion of a weighted sum of front arrival times, p∈X∪∂X φ(p)U ε (p) for some given φ : X ∪ ∂X → R, instead of a single value U ε (p * ).

Algorithmic strategy

In order to describe the application of sensitivity analysis to discretized eikonal equations, we introduce slightly modified notations. The perturbed Hamiltonian is rewritten as

H ε p ( ṗ) ≈ h -2 q∈V(p) exp (-2α ε (p, q)) p -q, ṗ 2 + , (41)
where ε ∈ R is a small parameter. This description is made equivalent to (32, left) by introducing the discretization stencils V(p) := {p -h ėi (p); 1 ≤ i ≤ I} (see (34) for boundary conditions and obstacles), and defining α ε (p, q) := -1 2 ln α ε i (p) when q = p -h ėi (p). Denote by U ε : X ∪ ∂X → R the unique solution to the perturbed and discretized eikonal PDE

∀p ∈ X, h -2 q∈V(p) exp(-2α ε (p, q))(U ε (p) -U ε (q)) 2 + = 1, ∀p ∈ ∂X, U ε (p) = σ ε (p). ( 42
)
Our objective is to relate the first order Taylor expansions of the boundary conditions σ ε , weights α ε , and solution U ε . Before we begin, let us mention that the simplified numerical scheme (32) considered in this section makes U ε differentiable, but that it is only almost everywhere differentiable in the general case (17). That is because a ∈ R → max{0, a} 2 is everywhere differentiable, but not (a, b) ∈ R 2 → max{a, b} 2 . For all p ∈ ∂X, the boundary condition (42, right) yields the simple relation

U ε (p) = σ ε (p)
where the prime denotes differentiation w.r.t. the parameter ε. For interior points p ∈ X, we obtain differentiating (42, left) that

q∈V(p) ω ε (p, q) U ε (p) -U ε (q) -(U ε (p) -U ε (q)) α ε (p, q) = 0,
where ω ε (p, q) := exp(-2α ε (p, q))(U ε (p) -U ε (q)) + . Note, crucially, that ω ε (p, q) > 0 requires U ε (p) > U ε (q). This property directly comes from the causality of the PDE discretization (42, left), i.e. the fact that it only involves positive parts of finite differences. Thanks to this, U ε (p) can be expressed in terms of those U ε (q) associated to points q ∈ V(p) reached by the front propagation strictly earlier than p. More precisely

U ε (p) = q∈V(p) ω ε (p, q) U ε (q) + (U ε (p) -U ε (q))α ε (p, q) , (43)
where ω ε (p, q) := ω ε (p, q)/ q∈V(p) ω ε (p, q), for any p ∈ X and any q ∈ V(p).

Let us reformulate this procedure in linear algebra terms. Let µ := (U 0 (p n )) N n=1 collect the derivatives of the solution, where the points of X ∪ ∂X are sorted by increasing values U 0 (p 1 ) ≤ • • • ≤ U 0 (p N ). We regard the perturbations ξ and ζ, of the cost functions and the boundary conditions, as column vectors. Thus ζ collects the values (σ 0 (p)) p∈∂X , and ξ is linearly related to (α 0 (p, q)) p∈X q∈V(p) . Therefore (43) at the parameter ε = 0 can be rewritten in the form

µ = Lµ + Aξ + Bζ,
where L is a strictly upper triangular matrix (in view of (43) and recalling that ω ε (p, q) > 0 requires U ε (p) > U ε (q)), and A and B are known matrices. Forward automatic differentiation computes µ := (Id -L) -1 (Aξ + Bζ), given ξ and ζ. Reverse automatic differentiation is the adjoint procedure, computing ρ := φ(Id -L) -1 A and π := φ(Id -L) -1 B, given φ : X ∪ ∂X → R. The triangular structure of Id -L is of course leveraged for fast inversion, by recursive substitution, and actually the matrices A, B and L are never explicitly assembled by the HFM software.

Numerical experiments

This section is devoted to numerical experiments, which are reproducible using the series of notebooks that come alongside this publication. Part of these are Jupyter notebooks written in the Python R language13 , while others are designed for the Mathematica R software14 . The original source code15 also features additional examples. We illustrate the different metric models introduced in §2, such as isotropic, Riemannian, or non-holonomic metrics, as well as the algorithmic techniques presented in §3, such as forward and reverse automatic differentiation.

The numerical examples discussed in this section were designed with the intent to illustrate the functionalities of the HFM library. In contrast, section §5 shifts the focus from algorithmics to modeling, and involves examples and test cases coming from applications.

Base functionalities, with an isotropic metric

This paragraph serves as an introduction to using the HFM library, on a problem involving a basic Isotropic metric, and that could thus be addressed using a variety of other software packages. Features more specific to our software are considered in the next subsections. We address here a two-dimensional shortest path problem involving an isotropic metric on a domain Ω ⊂ R 2 , with the following value function

u(p) := min γ(1)=p γ(0)∈∂Ω σ(γ(0)) + 1 0 c(γ(t)) γ (t) dt (44)
Such optimization problems are specified to the HFM library using a dictionary, which entries are in this paper denoted by key:value. The exact syntax for constructing such an object depends in practice on the language used for interfacing, e.g. Python, Matlab R or Mathematica R . This dictionnary is then fed to the adequate executable, for instance MatlabHFM Isotropic2, where the prefix denotes the external interface ("Matlab" could be replaced with "Python", "Mathematica", or "File" for the the stand-alone file based executable), and the suffix denotes the metric model, here an isotropic metric on a two-dimensional domain, as decided for (44). The next step is to define a cartesian grid, in our example of size 2n × n where n = 100, and a box domain Ω 0 = [-1, 1] × [0, 1], enclosing the PDE domain Ω, and containing the seed points where finite boundary conditions are imposed, and possibly some obstacles. For that purpose we provide the keys dims:(2n, n), origin:(-1, 0) (the bottom left corner), and gridscale:h, where h := 1/n. A few adjustments to this construction might unfortunately be necessary, depending on the visualization software, see the technical discussion in §D.

We next introduce some starting points for the front propagation, by e.g. the key pair seeds:[(-0.5, 0.3), (0.5, 0.8)], and some boundary conditions at these points, e.g. seedValues:[0,0.5]. We also need to define the cost function, e.g. the constant cost:1, and opt in or out of the second order enhancement to the fast marching method, e.g. sndOrder:1. At this point, we have gathered enough data to run the fast marching algorithm, but the HFM library must also be instructed what to export in return. We here choose to request minimal geodesics towards a family of points, e.g. tips:[(0,0.6),(-0.9,0.5),(0.8,0.8)], as well as the numerical solution U to the eikonal equation and the upwind geodesic flow vector field (39), as specified by the key pairs exportValues:1 and exportGeodesicFlow:1. The outputs of this program are illustrated on Figure 10.

In order to enrich the problem, we introduce some obstacles in the domain and a position dependent cost function. They are specified using the key pairs walls:arrBool and cost:arrFloat (replacing the former cost:1), where arrBool and arrFloat are arrays of boolean and of positive floating point values, sized according to discretization grid. See Figure 11.

The front propagation can be stopped early in many applications, for instance as soon as the geodesic(s) of interest can be backtracked, saving substantial computation time: often up to 80% in applications to tubular structure segmentation, see §4.5 . Fortunately, the single pass nature of the fast marching algorithm makes it particularly simple to define and implement termination criteria. Recall that the points of the discretization grid are successively "accepted", and the corresponding values frozen, in the order

U (p σ(1) ) ≤ • • • ≤ U (p σ(n)
) of increasing values of the numerical solution U . We may instruct the HFM library to stop the front propagation when all or any of a set of points are accepted, which is enough for geodesic backtracking, using the keys stopWhenAllAccepted or stopWhenAnyAccepted respectively. See Figure 12 (left). One may also put a threshold on the number of accepted points, or on the values of the solution U .

Several side products may be computed simultaneously while the fast marching algorithm is running, for use within additional stopping criteria and/or for being returned to the user. For instance, the computation of the Voronoi diagram associated with the different seeds is triggered by assigning labels of them, with e.g. the pair seedFlags:[0,1]. See Figure 12. Identical labels may be assigned to several seeds, in which case the regions are merged. Optionally, front propagation can be stopped when the Voronoi regions meet, using VoronoiStoppingCriterion:'RegionsMeeting'. In that case the minimal geodesic from the seeds of the meeting regions is returned, in the form of two halves joining at the regions meeting point. Another side product of interest is the Euclidean length of the minimal geodesics, see Figure 12 (right).

The use of a time dependent cost function c = c(p, t), giving rise to the eikonal equation c(p, u(p)) ∇u(p) = 1, is illustrated on Figure 13. We use an explicit scheme, which is simpler than [START_REF] Vladimirsky | A fast implicit method for timedependent Hamilton-Jacobi PDEs[END_REF] but introduces an O(h 2 ) overall error in the solution (expected to be dominated by the discretization error, since the numerical scheme used is at best second order). Other classical models implemented in the HFM library are also illustrated on Figure 13, namely two and three dimensional, isotropic and diagonal metrics, corresponding to the executables ending with 'Isotropic3', 'Diagonal2', or 'Diagonal3'. Finally, we recall that features presented in this subsection are transversal and applicable to any minimal path model implemented in the HFM library.

Automatic differentiation

Automatic, or algorithmic, differentiation, is a meta-programming technique for evaluating the Jacobian of a numerical function. We refer to §3.3 for details on these methods, and only recall here their basic intent. Denote by u = F (c, σ) the mapping which associates the numerical PDE solution u to the numerical cost function c and boundary conditions σ. Denote by J be the Jacobian matrix of F at some fixed input (c, σ), assuming it is indeed differentiable there. The matrix J is usually non-sparse and excessively large. Fortunately, the technique of forward (resp. backward) differentiation lets us evaluate at a very reasonable cost any matrix-vector product µ = J • (ξ, ζ) with the Jacobian (resp. (ρ, π) = J T • φ with its adjoint).

In the context of the fast marching algorithm, the first reported uses of forward and backward differentiation respectively appear in [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF] and [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF]. For simplicity, these functionalities are illustrated here in the context of isotropic metrics, but they are applicable to all the minimal path models implemented in the HFM library, including those featuring dimension lifting and curvature penalization techniques [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF].

Forward differentiation is used in Figure 14 to compute the first order perturbation in the distance map u, constructed in §4.1, subject to a given perturbation of the input cost function (left), or of the boundary conditions (center). Conversely, backward differentiation yields the sensitivity of a specific pointwise value u(p * ) to arbitrary modifications in the cost and boundary conditions, which reveals the support of the minimal geodesic(s), as discussed §3.2.2 and illustrated in Figure 14 (right).

We next present an optimization problem in which one aims to maximize the distance between two points p 0 and p 1 , by adjusting the cost function c : Ω →]0, ∞[ of an isotropic metric. See [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF][START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF] for problems of similar nature, some of them more complex. The cost function is constrained by upper and lower bounds, and subject to an integral penalty term, all determined by positive constants α, β, γ. The problem, which turns our to be convex, reads as follows:

max α≤c≤β min γ(0)=p 0 γ(1)=p 1 1 0 c(γ(t)) γ (t) dt -γ Ω c(x)dx. (45)
We use the same domain and obstacles as in the previous section, on a 100 × 200 cartesian grid, but switch to the differentiable discretization scheme (20) for isotropic eikonal PDEs.

The optimization is solved in approximately 20 seconds using the LBFGS algorithm16 . For comparison, solving similar instances reportedly took hours in [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF]. The lower complexity of our implementation mainly comes from the use of backward differentiation, instead of forward differentiation in [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF]. The optimal cost function c for (45) has several interesting features.

For instance it creates "barriers" in the space between close obstacles. Another peculiarity is the fact that there are uncountably many distinct paths of minimal cost, w.r.t. c, joining the endpoints p 0 and p 1 . See Figure 15.

Riemannian metrics

Riemannian metrics are the most widely occurring class of non-isotropic metrics in applications.

Extending the fast marching algorithm to this context is a non-trivial task, which has been the subject of a continued line of research since the 2000's. Many of the proposed approaches have significant drawbacks, such as the use of excessively large stencils [KS98, SV01, AM11], the pre-condition that the condition number of the metric tensors be mild [JBT + 08], or the fact that they give up on the causality property [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF].

The author has introduced efficient and specific techniques for discretizing Riemannian eikonal PDEs on cartesian grids. Two flavors were developed, in the semi-Lagrangian and Eulerian framework respectively, referred to as17 the FM-LBR and FM-VR1. They are introduced in [START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF] and [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF] respectively, with open source code distribution and reproducible experiments in [START_REF] Mirebeau | Anisotropic Fast Marching in ITK[END_REF] and the present manuscript respectively. For those readers who may hesitate between the two approaches, let us say that the FM-VR1 is faster, especially in dimension 3, and more accurate if the second order enhancement of the numerical scheme is activated, but slightly less accurate otherwise, see [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF] for a more in-depth comparison. From a practical point of view, the FM-VR1 code is more feature packed, and has interfaces to more scripting languages, while the FM-LBR is tightly integrated within the Insight ToolKit (ITK R ) library.

We repeat in this subsection the numerical experiments presented in [START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF][START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], which can now be easily reproduced using the provided Python and Mathematica R notebooks. The test cases are all synthetic, but are inspired from problems arising in geometry processing, seismic imaging and tubular segmentation. In practice, instructing the HFM software to numerically solve a Riemmanian eikonal equation is done by selecting the appropriate executable, ending with 'Riemann2' or 'Riemann3', depending on the dimension18 , and providing the key-value pair metric:arrTensors, where "arrTensors" is an array containing the Riemannian metric tensors. (Alternatively, the inverse tensors may be supplied as dualMetric:arrTensors.)

A series of two dimensional test cases, considered in [SV01, Mir14, Mir17a], are described below. Level sets of each example are displayed in Figure 17. We recall that the condition number of a symmetric tensor is defined by Cond(M ) := M M -1 . The numerical techniques implemented in the HFM library are primarily intended Riemannian metrics which tensors have condition number 10, although more degenerate instances are often handled successfully.

1. (Geometry processing inspired) We compute the distance from the origin on a two-dimensional manifold, embedded in R 3 and defined by the following height map, where x, y ∈ [-1, 1] 2 z(x, y) := (3/4) sin(3πx) sin(3πy).

The Riemannian metric is given by M(x, y) = I + ∇z(x, y)∇z(x, y) T , and has condition number 6. See Figures 16 and17. A particularity of this test case is that the eigenvectors of the metric are quite often almost aligned with the coordinate axes. We thus rotate them by π 6 , in a second test case, to show that our numerical method is not sensitive to this bias, contrary to several others, see the discussion in [START_REF] Mirebeau | Anisotropic Fast-Marching on Cartesian Grids Using Lattice Basis Reduction[END_REF].

(Seismic imaging inspired)

We consider the distance-to-origin problem for the metric tensor M(x, y) with eigenvalue 0.8 -2 associated with the eigenvector (1, π 2 cos(4πx)), and eigenvalue 0.2 -2 associated with a perpendicular eigenvector, on the domain [-1, 1] 2 . The condition number of the Riemannian metric is thus bounded by 4.

(Tubular segmentation inspired)

We illustrate the robustness of our numerical scheme by choosing a Riemannian metric which is both discontinuous and extremely anisotropic. More precisely, the metric is isotropic euclidean sufficiently far from a predefined curve, in this case a spiral, and close to the curve it is highly anistropic with eigenvalues (1, 100 -2 ). The tangent to the curve is the eigenvector corresponding to the small eigenvalue.

Three dimensional counterparts of the "seismic imaging" and "tubular segmentation" inspired test cases are also presented in Figure 16 (right), and in the Python notebooks.

Planar curvature penalization

An important innovation of the HFM library is the ability to find planar curves globally minimizing an energy involving their curvature. For that purpose, the problem is lifted in the configuration space R 2 × S 1 of positions and orientations, following an idea introduced in [START_REF] Pechaud | Extraction of tubular structures over an orientation domain[END_REF]. A suitable non-holonomic metric is then introduced on this set, as proposed in [START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF][START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF] for the Reeds-Shepp and Euler-Mumford models. Lastly, a specific generalization of the fast marching method is required for efficiently solving these problems numerically, see [SBD + 15, DMMP16], and [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF][START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] which are at the foundation of our software.

The HFM library implements four curvature penalized path models, whose distinctive features are listed below. See §1.2 and §2.4 for a formal definition.

• The classical Reeds-Shepp sub-Riemannian model, encodes the motion constraints of a wheelchair, which is able to move both forward and backward. Distinctive feature: the trajectories feature cusps, at places where the vehicle changes from forward to reverse gear.

• The Reeds-Shepp forward variant is mathematically similar, except that the vehicle cannot go backwards. Distinctive feature: the trajectories feature in-place rotations, in particular at their endpoints and around obstacles.

• The Euler-Mumford elastica model describes the rest position of an elastic bar19 . Distinctive feature: the trajectories are smooth, and "feel natural".

• The Dubins car model has a bounded turning radius. Distinctive feature: the optimal trajectories are concatenations of straight lines and of arcs of circle.

These models are illustrated on Figure 18 in the special case of a constant cost function c ≡ 1.

On Figure 19 we consider both a position dependent cost function c = c(x), which is a common case in applications, and an orientation dependent cost function c = c(θ), which in particular is suitable for modeling a sailboat. The exectuable ending with 'ReedsShepp2', 'ReedsShepp-Forward2', 'Elastica2' or 'Dubins2' must be selected to solve these instances. The parameter ξ, homogeneous to a radius and modulating the intensity of curvature penalization see (9), must be assigned a value, such as xi:1. The cost function is set with the pair cost:arr where 'arr' is either a constant, a one-dimensional array if c = c(θ), a two-dimensional array if c = c(x), or a three dimensional array in the general case where c = c(x, θ), of suitable dimensions.

The HFM software also implements generalizations of the above four models, in which the curvature penalization parameter may vary over the domain ξ = ξ(x, θ), and an additional term κ = κ(x, θ), homogeneous to a curvature, introduces some asymmetry in the treatment of right and left turns. See (10) and Figure 20.

Computational efficiency is an important aspect of numerical techniques for minimal path extraction, which severely constrains their range of applications (response times of different orders of magnitude are required for e.g. real time control of a dynamical system, convenient interaction with a human within an image processing software, or offline data processing). On Figure 7, we compare the accuracy and execution time of four numerical methods for the classical Reeds-Shepp model, extending the work [SBD + 15]. Are displayed: the HFM library with (red) or without (blue) second order enhancement, the semi-Lagrangian FM-LBR implementation (yellow) [SBD + 15], and a non-causal scheme solved using an iterative numerical method (green) [START_REF] Bekkers | A PDE Approach to Data-Driven Sub-Riemannian Geodesics in SE(2)[END_REF]. A cartesian grid of size (4n, 4n, 2n) is used to discretize the domain [-2π, 2π] 2 × S 1 , and the cost function is set constant c ≡ 1. The exact solution to the eikonal equation, computed using a geodesic shooting method, is used for reference.

The first three compared numerical methods involve a relaxation parameter, here set to the default value ε = 0.1. Therefore one could expect the numerical error to plateau at sufficiently high grid scales n. However this phenomenon is not observed in our experiment, which only shows straight lines in log-log scale. This suggests that the discretization error introduced by the relaxation is not a limiting factor for our numerical scheme, at reasonable grid sizes.

At a given image resolution n, the most accurate numerical method is the FM-VR1 scheme with second order enhancement, which is implemented in the HFM library. The (slow) iterative method comes next, and is best among the first order schemes, followed by the FM-LBR scheme, and finally (a bit disappointingly) the first order instantiation of the FM-VR1. Given a Log-log scale. In both figures we compare an iterative approach (green), the FM-LBR (yellow), the FM-VR1 implemented in the HFM-library (blue), and the same with second order enhancement (red). The latter turns out to be the most efficient, see §4.4 for a discussion. computation time budget, the FM-VR1 with second order enhancement offers the best accuracy. The first order FM-VR1 comes second, compensating its lower accuracy resolution-wise w.r.t. the alternative methods by its much better complexity. The FM-LBR comes third, and the non-causal iterative method comes last due to its much longer computation times, which are superlinear w.r.t. the number of grid points.

Reeds-Shepp models in

R 3 × S 2
The HFM library can extract three dimensional curves globally minimizing a curvature dependent energy. For that purpose, the curves are lifted in the five dimensional configuration space R 3 × S 2 , and the second order path energy functional is rewritten in terms of a suitable non-holonomic metric, similarly to the planar case albeit in higher dimension. Numerically, an adequate eikonal PDE is solved, after what the optimal curves are backtracked. Incidentally, the workflow that we just described involves numerically solving a non-linear PDE with degenerate anisotropy on a five dimensional manifold with non-trivial topology, which may look like a fairly challenging problem. Convincing results are nevertheless obtained in our examples, which run within 5s to 40s on a standard laptop.

Three minimal path models are considered in this section, which are higher dimensional generalizations and variants of the Reeds-Shepp model. The original sub-Riemannian model distinguishes itself by the presence of cusps, where the path orientation is reversed. The forwardonly variant lacks cusps, but often features in place rotations, in particular at the minimal path endpoints, see [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF] for a discussion. Finally, a dual variant favors paths which may be non-smooth, but are embedded in smooth manifolds. It can be related with torsion penalization, rather than curvature penalization, see [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF]. See also [START_REF] Strandmark | Shortest Paths with Curvature and Torsion[END_REF] for a distinct attempt to implement genuine torsion penalization. Note that the three dimensional analogues of the Euler-Mumford and Dubins models, considered in the planar case §4.4, are not implemented in the HFM library at the time of writing.

The non-holonomic metrics considered involve a data-driven cost function c : R 3 ×S 2 →]0, ∞[, a parameter ξ > 0 homogeneous to a radius of curvature, and in our case a relaxation parameter20 ε > 0. For any point (x, n) ∈ R 3 × S 2 , and any tangent vector ( ẋ, ṅ), the metric associated with the Reeds-Shepp model, the forward variant, and the dual variant, respectively read

F ε (x,n) ( ẋ, ṅ) 2 = c(x, n) 2 n, ẋ 2 + ε -2 n × ẋ 2 + ξ 2 ṅ 2 , F ε,+ (x,n) ( ẋ, ṅ) 2 = c(x, n) 2 n, ẋ 2 + + ε -2 n, ẋ 2 -+ ε -2 n × ẋ 2 + ξ 2 ṅ 2 , F ε,× (x,n) ( ẋ, ṅ) 2 = c(x, n) 2 ε -2 n, ẋ 2 + n × ẋ 2 + ξ 2 ṅ 2 .
We denoted by n × ẋ the cross product of the three dimensional vectors n ∈ S 2 and ẋ ∈ R 3 . In the limit as ε → 0, the angular orientation n and the physical velocity ẋ are respectively constrained to be collinear, positively collinear, and orthogonal if the path cost remains bounded, due to ε -2 penalization terms.

Our implementation relies on a parametrization of the sphere S 2 using two Euler angles θ ∈ [0, π] and φ ∈ [0, 2π], with appropriate boundary conditions: n(θ, φ) := (cos θ, sin θ cos φ, sin θ sin φ).

The Reeds-Shepp model and the dual variant, but not the forward variant, also make sense if n belongs to the projective orientation space P 2 = S 2 /{-1, 1}, which is parametrized as above except that θ ∈ [0, π/2] only. The same angular spacing h φ = 2π/n φ is used for discretizing the angles θ and φ, hence the cartesian grid devoted to the angular S 2 (resp. P 2 ) space has dimension (n φ /2, n φ ) (resp. (n φ /4, n φ )). See Figure 8.

The HFM library executables ending with 'ReedsShepp3' or 'ReedsSheppForward3' must be selected to solve these problems. In the case of the Reeds-Shepp model, the projective angular space is imposed by the key-value pair projective:1, and the dual model is activated by dual:1. The curvature penalization and the relaxation parameter need to be set, e.g. xi:1 and eps:0.2. In the examples presented, the discretization grid size is defined as dims:[60, 40, 25, n θ , 24], where n θ = 6 for the projective angular space P 2 , and n θ = 12 in the case of S 2 , see Figure 8. The cost function is specified by cost:arr, where arr is either a constant, a two-dimensional array (c = c(n)), a three dimensional array (c = c(x)), or a five-dimensional array (c = c(x, n)), of suitable dimensions.

A test case inspired by tubular structure segmentation is presented in Figure 21. The cost function is small in the neighborhood of two curves, which respectively have low curvature and low torsion, and is large elsewhere. The common endpoints of the curves are respectively used as the seed of the front propagation, and the tip from which to backtrack a geodesic. The minimal paths selected by the Reeds-Shepp and the Reeds-Shepp dual model go along the low curvature and low torsion curve respectively, as was expected and desired. A second test case, inspired by motion planning, is displayed on Figure 22. The backtracked curves are in that case smooth and with well distributed curvature, up to the occasional and expected singularities: cusps for the Reeds-Shepp model, and in-place rotations for the forward variant. Computation times21 are respectively around 5s and 40s in the tubular segmentation and motion planning related experiments, despite the similar grid sizes, because the front propagation can be aborted earlier in the first case.

Selected applications

In this section, we present selected applications of the HFM library. The difference with the numerical experiments discussed the previous section lies in the intent and focus: modeling aspects are more emphasized here, and the datasets come from other application fields (although some are synthetical), instead of being specifically designed so as to illustrate the functionalities of the HFM library.

An interpretation of Poggendorff 's visual illusions

We present two visual illusions due to Poggendorff, and their explanation according to the work [START_REF] Franceschiello | Modelling of the Poggendorff Illusion via Sub-Riemannian Geodesics in the Roto-Translation Group[END_REF]. This research builds on the works [START_REF] Petitot | The neurogeometry of pinwheels as a sub-Riemannian contact structure[END_REF][START_REF] William H Bosking | Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex[END_REF], which have shown that the first layer V1 of the visual cortex could be regarded as a biological implementation of the manifold R 2 × P 1 , equipped with the sub-Riemannian structure defining the Reeds-Shepp model.

Consider an image displaying a curve ending at a point x * ∈ R 2 with tangent orientation θ * , and a curve starting at another point x * with tangent orientation θ * . Under some conditions, a human presented with this image will associate the two curves and infer a connection between them [START_REF] Duits | Association Fields via Cuspless Sub-Riemannian Geodesics in SE(2)[END_REF]. In everyday experience, a similar phenomenon helps us guess the shape of objects which boundary is partially occluded. The curve reconstructed by the visual system is, according to the works cited above, the (planar projection of the) minimal path joining the oriented endpoints p * := (x * , θ * ) and p * := (x * , θ * ) w.r.t. the Reeds-Shepp model. The parameter ξ, which determines the amount of curvature penalization in this model, depends on the scale at which the image is displayed, and is adjusted by hand in our experiments. The cost function is chosen constant c ≡ 1 in our experiments, although a more complex data-driven construction is considered in the original work [START_REF] Franceschiello | Modelling of the Poggendorff Illusion via Sub-Riemannian Geodesics in the Roto-Translation Group[END_REF].

Assume now that several endpoints (x * i , θ * i ) i∈I are present in the image, instead of a single one. Then the visual system infers a curve joining the source point (x * , θ * ) to the closest endpoint, again w.r.t. the distance defined by the Reeds-Shepp model, a process called perceptual grouping [START_REF] Bekkers | Nilpotent Approximations of Sub-Riemannian Distances for Fast Perceptual Grouping of Blood Vessels in 2D and 3D[END_REF]. The images created by Poggendorff put this physiological process in evidence by displaying a partially occluded straight line or circle. To the eye, the true endpoints of the geometrical figure seem bizarrely misaligned, and we would like to move one slightly along the boundary of the occluded region, see Figure 23. This is a visual illusion, highlighting the difference between (a) the Reeds-Shepp geodesic toward the closest point across the occluded region (with the correct orientation), and (b) the continuation by a straight segment or an arc of circle that would be "logical" in the presented image.

Motion planning

We illustrate on Figure 24 a two dimensional motion planning problem, namely finding the optimal path towards the exit of a building22 . The Isotropic, Reeds-Shepp, Reeds-Shepp forward, Elastica and Dubins models are compared, and their distinctive features are clearly visible, see §4.4.

The authors readily acknowledge that shortest paths (w.r.t. any reasonable metric) are a quite questionable approximation of natural paths, in view of the potential tradeoffs with safety, comfort or computational effort. For instance, the optimal computed trajectories are tangent to the wall corners, which would be hazardous in practice in case someone comes in the opposite direction, and anticipate obstacles before they come in the line of sight. Paths closer to e.g. natural pedestrian motion are obtained by penalizing motion close to the walls, see again Figure 24. See [START_REF] Cartee | Anisotropic Challenges in Pedestrian Flow Modeling[END_REF] and references therein on the topic of pedestrian crowd motion models based on eikonal equations. Optimal paths also do make sense in other motion planning contexts: if one is playing a race, or is trying to avoid a surveillance system [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF].

Tubular structure segmentation

Tubular structure segmentation in medical images is one of the main applications intended for the HFM library. Early versions, and related software of the first author, are used for that purpose in [CCM14, SBD + 15, CMC16a, DMMP16, CMC16b, MDS + 17, BCP17]. Vessel extraction in images of the retina is illustrated in Figure 25, using two popular models involving respectively radius and orientation lifting. Let us briefly recall the principles underlying this approach.

Consider an image, e.g. in grayscale and represented by a function I : U → [0, 1], displaying a family of possibly overlaid tubular structures. Segmentation methods based on minimal paths attempt to extract the centerlines of these tubular structures as geodesics w.r.t. a suitable metric, joining endpoints which are either prescribed by the user or automatically detected by another method. They are frequently combined with dimension lifting techniques: a cost function c : U × A →]0, ∞[ is defined over the product of the image domain U with an abstract parameter space A. In practice, the latter often represents a range of radii A = [r, R], of grayscale intensities A = [0, 1], or of angular orientations A = S 1 , see e.g. [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] Pechaud | Extraction of tubular structures over an orientation domain[END_REF]. The local cost c(x, a) is data-driven and built during a preliminary filtering of the input image, see e.g. [START_REF] Duits | Image Analysis and Reconstruction using a Wavelet Transform Constructed from a Reducible Representation of the Euclidean Motion Group[END_REF]. The quantity c(x, a) should be small iff a tubular structure having feature a is likely present at the position x in the image.

The purpose of dimension lifting is two-fold. First: robustly extract the additional feature of interest, in addition to the vessel centerline, for applications in medical diagnostic. Vessel radius and curvature, inferred from orientation in the latter case, are particularly important in this respect. Second: disentangle the overlaid tubular structures, by cleanly separating them in three dimensional space. This helps avoid or reduce the "shortcut" and "leak" type artifacts, which often plague minimal path methods for tubular structure extraction. Figure 25 illustrates radius lifting, combined with an isotropic metric model, and orientation lifting, combined with the Reeds-Shepp model.

Segmentation in 3D DMRI data

We illustrate three dimensional tubular structure segmentation in (simulated) dMRI data. This medical imaging technique measures diffusivity at each point x ∈ R 3 , and in each direction n ∈ S 2 . The input data is therefore intrinsically defined over the configuration space R 3 × S 2 , in contrast with §5.3 where dimension lifting was an artificially introduced image processing technique.

We use a digital phantom dataset that was constructed, using the dMRI simulation method [CDD + 14], for the ISBI 2013 reconstruction challenge, where it was used as a benchmark for tracking methods [DCRD + 14]. The dataset itself consists of a fairly challenging configuration of 27 simulated white matter bundles, inside a spherical volume. The entire volume has a size of 50 × 50 × 50 voxels. In the same volume, there are three spherical regions with isotropic diffusion.

We use a constrained spherical deconvolution approach [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution[END_REF] to estimate from the raw simulated dMRI data a fiber orientation distribution (FOD), that in each voxel indicates the probability of having a fiber in that direction. The density of this FOD is denoted by f : R 3 × S 2 → R + . Then following the approach in [START_REF] Duits | Optimal Paths for Variants of the 2d and 3d Reeds-Shepp Car with Applications in Image Analysis[END_REF], we use the FOD to construct the cost function as follows: for any p

= (x, n) ∈ R 3 × S 2 c(p) = 1 + σ f (p) ||f || ∞ p -1 , σ, p ≥ 1. (46)
The parameter σ influences how strongly the data should affect the optimal path, and the exponent p influences how 'sharp' the FOD profiles are. As a proof-of-concept experiment on this dataset, we consider the case where the seed point and end point of a bundle are given and we compute the shortest path connecting the two. The path should then at least stay inside the bundle. We use the Reeds-Shepp forward model, with parameters parameters σ = 3, p = 1, ξ = 0.5, ε = 0.1, n θ = 10 and n φ = 20, resulting in 200 points in the discretization of the sphere. The geodesics (in black) stay in almost all cases entirely inside the volume of the corresponding bundle, as is displayed in Fig. 9.

Conclusion

In this paper, we introduced the Hamiltonian Fast Marching (HFM) library, a state-of-the-art software for global path optimization with respect to non-isotropic metrics. Our numerical methods apply in particular to Riemannian metrics in dimensions two and three, possibly strongly anisotropic, using an adaptive discretization scheme based on lattice geometry techniques. We also globally optimize second order energies depending on the path curvature, by introducing a dimension lifted space equipped with a non-holonomic metric, suitably relaxed in the implementation. The classical models due to Reeds-Shepp, Euler-Mumford, and Dubins, are implemented, as well as numerous generalizations and variants which all have specific qualitative properties and use cases. A number of transversal methods are provided, such as backward algorithmic differentiation which is novel in this context. This software is provided in the form of a unified, optimized, and open source c++ code, with convenient interfaces for several major scripting languages.

The HFM library has already found several applications in the study of visual illusions, motion planning, and image segmentation. We hope that this publication, and the series of introductory notebooks that are provided alongside, will help popularize these techniques and their present and future applications. The HFM library is part of an ongoing research effort and will continue to evolve. Considered additions include higher dimensional Riemannian metrics and a better enforcement of non-holonomy constraints.

A Dual metric for curvature penalization

This appendix is devoted to formally deriving the expression of the dual metric associated with the Reeds-Shepp forward, Euler-Mumford, and Dubins model. It is announced in §1.2 and §2.4. We use a unified and systematic presentation that easily generalizes to e.g. the asymmetric variants (10), or the higher dimensional instantiations. See [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] for another approach.

For notational simplicity, let us fix a point (x, θ) ∈ R 2 × S 1 of the configuration space of these models, and denote F := F (x,θ) , c := c(x, θ) and n := (cos θ, sin θ). Our first step is to justify the semi-explicit dual metric expression (8, right). Indeed, one has

F * (x, θ) = sup ( ẋ, θ) =0 x, ẋ + θ θ F ( ẋ, θ) = sup ( ẋ, θ) =0 ẋ≥0
x, ẋn + θ θ

F ( ẋn, θ) = sup θ∈R x, n + θ θ c C( θ) = f ( x, n , θ) c .
We used successively (i) the definition (4) of the dual metric, (ii) the fact that F ( ẋ, θ) = ∞ unless the physical velocity is a non-negative multiple of the current orientation n, (iii) the 1-homogeneity of F , allowing to assume23 that ẋ = 1, and its explicit expression (8, left), and (iv) the definition (28) of the function f . Our next objective is to compute the supremum, denoted f (x, θ), of

θ ∈ R → (x + θ θ)/C( θ). (47)
For that purpose we recall the specific instantiations C RS , C EM and C D of the cost function C penalizing curvature, see (9). For notational simplicity, and up to rescaling θ, we assume ξ = 1.

C RS ( θ) = 1 + θ2 , C EM ( θ) = 1 + θ2 , C D ( θ) = 1 if θ ≤ 1, ∞ else.
In the Dubins case, the map (47) is linear on the interval [-1, 1], and vanishes elsewhere, hence the maximum is f (x, θ) := max{x -θ, x + θ, 0}. In the Reeds-Shepp and Euler-Mumford cases, the supremum of (47) is attained either as θ → ±∞, or for a finite value θ ∈ R obeying the optimality condition θC( θ) = (x + θ θ)C ( θ).

For the Reeds-Shepp curvature cost, observing that C RS ( θ) = θ/ 1 + θ2 we obtain

θ 1 + θ2 = (x + θ θ) θ 1 + θ2 ⇔ θ(1 + θ2 ) = (x + θ θ) θ ⇔ θ = x θ.
Thus ẋ = θ/x, provided x = 0, and therefore f (x, θ) = (x + θ θ x )/ 1 + ( θ x ) 2 = x2 + θ2 as announced (since θ → ±∞ is sub-optimal). In addition f (0, θ) = sup θ∈R θ θ/ 1 + θ2 = |x|, in the limit θ → sign( θ)∞, hence the previous expression still holds in the case x = 0. For the Euler-Mumford curvature cost, observing that C EM ( θ) = 2 θ we rewrite (48) as θ(1 + θ2 ) = (x + θ θ)2 θ ⇔ θ2 + 2 θ x θ -1 = 0 ⇔ θ = θ± = -

x θ ± 1 + ( x θ ) 2 , (49)
assuming θ = 0. Inserting these roots in (47) we obtain

x + θ θ±

C EM ( θ± ) = θ C EM ( θ± ) = θ 2 θ± = -θ θ∓ 2 = x ∓ x2 + θ2 2 , ( 50
)
where we used successively (i) the identity (48) defining the roots, (ii) the expression C EM ( θ) = 2 θ, (iii) the fact that θ+ θ-= -1, which follows from (49, center) and (iv) the explicit expression of θ± . Choosing "+" sign in (50, right) yields the largest value, hence f (x, θ) = (x+ x2 + θ2 )/2 (since θ → ±∞ is sub-optimal). In the limit case θ = 0, this expression becomes f (x, 0) = (x + |x|)/2, which is still correct, and is attained for θ = 0 if x > 0, and as θ → ∞ otherwise.

B Selling's algorithm

Our discretization of Riemannian eikonal equations strongly relies on a particular decomposition of positive definite tensors, see Proposition 2.2. In this Appendix, we describe its algorithmic computation in dimension d ∈ {2, 3}, as it is implemented in the HFM library. For that purpose, we need to introduce a few concepts from lattice geometry. Any superbase can be used to decompose any 2-tensor D, as follows

D = - 0≤i<j≤d e i , De j v ij ⊗ v ij , (51)
where v ij := e ⊥ k in dimension d = 2 and with {i, j, k} = {0, 1, 2}, and v ij := e k × e l in dimension d = 3 and with {i, j, k, l} = {0, 1, 2, 3}. See Lemma 4.4 in [START_REF] Mirebeau | Fast Marching methods for Curvature Penalized Shortest Paths[END_REF] for a proof. We are interested in non-negative tensor decompositions (23), which in the case of ( 51) is equivalent to a geometrical property of the superbase.

C Norm of the geodesic flow

This section is devoted to the proof of Proposition 3.1. Our first step is to rewrite the announced result in a more abstract framework, keeping only the minimal assumptions. Pixel area. The HFM library assumes, as shown above and in dimension two for simplicity, that the "pixel" of coordinates (i, j) occupies the square with corners (i, j)h and (i + 1, j + 1)h. However, some visualization software use define the "pixel" of index (i, j) as the square of side h centered at (i, j)h or at (i+1, j +1)h. This raises (minor) issues when geodesic paths are overlaid on images, and can be adjusted for by shifting the origin by (-0.5, -0.5)h or (0.5, 0.5)h.

Array ordering. When displaying the data contained in a two or three-dimensional array, some software choose to assign the array indices (i, j, k) to the coordinate axes (X, Y, Z) in a transposed manner: i → Y, j → X, k → Z. In addition, the array may be internally stored in row-major or column-major format. The arrayOrdering field accounts for these specificities, and must be set to 'RowMajor' (default value) for Python Mayavi R , 'YXZ RowMajor' for Python PyPlot, or' YXZ ColumnMajor' for Matlab R .

E Figures

Figure 1 :

 1 Figure 1: The four steps of minimal path computation. (I) Define a metric, here of Riemannian type and visualized as a family of ellipsoids. (II) Compute of a distance map, here from the center point, by solving the eikonal PDE (3). (III) Compute the geodesic flow direction p ∈ Ω → dF * p (du(p)). (IV) Backtrack the minimal geodesics, by solving the ODE (5).

1 0F

 1 γ(t) ( γ(t))dt, d F (p, q) := inf{L F (γ); γ ∈ Γ, γ(0) = p, γ(1) = q}. (1)

F

  p ( ṗ) := c(p) ṗ , F * p (p) = p /c(p).

Figure 3 :

 3 Figure 3: The stencil of a PDE numerical scheme is the collection of offsets, towards close neighbors usually, involved in its definition. We display the stencils used for the discretization of eikonal equations, associated with (left, center left) an isotropic metrics in dimension 2 and 3, (center right) an anisotropic Riemannian metric, (right) the Euler-Mumford elastica model.

Figure 4 :

 4 Figure 4: Dual unit ball { ṗ; ṗ D -1 ≤ 1}, and offsets associated by Voronoi's first reduction to a tensor D, involved in the decomposition (23). Left: several examples in dimension two. Right: a three dimensional case.

Figure 5 :

 5 Figure 5: Left: Applying Proposition 2.3 to a vector (black) with ε = 0.1, in dimension 3, yields the red offsets. Center: Discretization stencil for the Reeds-Shepp model, red offsets obtained by applying Proposition 2.3 in dimension d = 2, blue offsets related to the angular dimension. Right: Stencil for Dubins model, red and blue offsets obtained by applying Proposition 2.3 in dimension d = 3 to two distinct vectors.

  dH p (du(p)) = 2F * p (du(p)) dF * p (du(p)) = 2V (p), where we used first the identity H p (p) = F * p (p) 2 , and second the eikonal equation F * p (du(p)) = 1. Recalling that H p (p) = I i=1 α i p, ėi 2 + (left), and using a first order Taylor expansion (right), one obtains on the other hand

Figure 7 :

 7 Figure 7: Comparison of the resolution/accuracy (left) and accuracy/computation time (right) compromise of several numerical methods for the eikonal equation of the Reeds-Shepp model.Log-log scale. In both figures we compare an iterative approach (green), the FM-LBR (yellow), the FM-VR1 implemented in the HFM-library (blue), and the same with second order enhancement (red). The latter turns out to be the most efficient, see §4.4 for a discussion.

Figure 8 :

 8 Figure 8: Left: Discretization of the unit sphere S 2 , using 24 × 12 points, used for the presented numerical experiments in R 3 × S 2 . The next three pictures show the sets { ẋ ∈ R 3 ; F ε (0,n) ( ẋ, 0) ≤ 1}, with the relaxation parameter ε = 0.2 and for several orientations n ∈ S 2 , where F is successively the metric of the Reeds-Shepp model (needles), of the forward variant (half needles), and of the dual variant (plates).

Figure 9 :

 9 Figure 9: Visualization of the ISBI dataset from different viewpoints. The two left-most images show the full dataset, the other four images show a selection of the bundles. The black lines inside these bundles indicate the geodesics obtained with fast marching.

Definition B. 1 .

 1 A superbase of Z d is a (d+1)-tuple (e 0 , • • • , e d ) ∈ (Z d ) d+1 such that | det(e 1 , • • • , e d )| = 1 and e 0 + • • • + e d = 0.

i=0 e i 2 Di=0 e i 2 D + 1 2 i<j e i + e j 2 D

 2222 Definition B.2. A superbase (e 0 , • • • , e d ) of Z d is said D-obtuse, where D ∈ S ++ (R d ), iff M e i , e j ≤ 0 for all 0 ≤ i < j ≤ d. In dimension d ∈ {2, 3}, Algorithm 2 due to Selling [Sel74] takes a positive definite tensor D as input, and outputs a D-obtuse superbase. Algorithm 2 Selling's algorithm [Sel74] Input: A positive definite tensor D ∈ S ++ (R d ), an arbitrary superbase b = (e 0 , • • • , e d ). Output: A D-obtuse superbase. While there exists 0 ≤ i < j ≤ d such that e i , De j > 0 do If d = 2 then b ← (-e i , e j , e i -e j ), If d = 3 then b ← (-e i , e j , e k + e i , e l + e i ).Termination guarantee. Consider the energy E(e 0 , • • • , e d ), defined by 2 in dimension d = 2 and by 3 in dimension d = 3. This quantity decreases by 4 e i * , De j * > 0 at each iteration, where i * and j * are the indices appearing in the While condition, as can easily be checked by expressing the present and next superbase energy in terms of the scalar products ( e i , e j ) d i,j=0 . Since there exists only finitely many superbases which energy E is below any given bound, the algorithm must terminate.

Proposition C. 1 .

 1 Figure 10: Left: distance map u(p) := min{|p -p 1 |, σ 2 + |p -p 2 |}, where p 1 = (-0.5, 0.3), p 2 = (0.5, 0.8) and σ 2 = 0.5, numerically computed using isotropic fast marching with second order enhancement. Center: geodesic flow V (x) := ∇u(x) obtained numerically using (39). Right: verification that |∇V (x)| = 1, at all points of differentiability.

Figure 11 :

 11 Figure 11: Distance map obtained using isotropic fast marching (left), in the presence of obstacles (center), and of a non-constant cost function (right, level sets shown).

Figure 12 :

 12 Figure 12: Left: Front propagation stopped when all geodesics tips are accepted, which is enough for backtracking. Center: Voronoi regions associated with the two seeds. Right: The Euclidean length l(p) of the minimal path from each point p ∈ Ω, shown here, can be computed simultaneously with the distance map u(p). It is used as a stopping criterion in this experiment.

Figure 13 :

 13 Figure 13: Distance map associated with, respectively: a time dependent cost function (left), a three dimensional isotropic metric (center), or a diagonal metric (right).

Figure 14 :

 14 Figure 14: Left: First order term µ in the Taylor expansion u ε = u 0 + εµ + o(ε) of the distance map associated with the perturbed cost c ε = c 0 + εχ x<0 . Center: Likewise, for a perturbation of the boundary condition σ 2,ε = σ 2 + ε at the seed p 2 = (0.5, 0.8). Right: Sensitivity φ, with respect to the cost function, of the distance map u(p * ) at a point p * . The support of φ approximates the minimal geodesic(s) from p * to the seed points. Two geodesics are visible here since p * was chosen in the cut locus.

Figure 15 :

 15 Figure 15: Left: Cost function maximizing the distance between the seed (-0.7, 0.7) and tip p * = (0.8, 0.4), subject to bound constraints and a mean penalty, see (45). Center: Sensitivity of u(p * ) w.r.t. the cost function, at the optimum. Right: Distance map u from the seed, backtracked geodesics from the tip p * and nearby pixels, at the optimum.

Figure 16 :

 16 Figure 16: Left: Minimal geodesics on a parametrized surface in R 3 . Center: Top view of the same geodesics. Right: Level sets of the front, for a test case inspired by three dimensional tubular structure segmentation.

Figure 18 :

 18 Figure 18: From top to bottom: Reeds-Shepp, Reeds-Shepp forward, Euler-Mumford elastica, and Dubins car models. Left: Control sets, i.e. unit balls of the local metric. Center: Minimal paths in R 2 × S 1 . Right: Planar projection of the minimal paths.

Figure 19 :

 19 Figure 19: Left: Euler-Mumford elastica paths with a cost function c = c(x) depending only on the physical position x ∈ R 2 (smaller cost in the upper part). Center: Polar plot of a cost function c = c(θ), related with sailboat navigation, depending only on the orientation θ ∈ S 1 . Center: Minimal paths for the Reeds-Shepp forward model with the later cost function. The paths tack, as an actual sailboat would do, instead of going straight against the wind.

Figure 20 :

 20 Figure 20: Minimal paths for variants of the Dubins car. Left: The car has a position dependent turning radius ξ = ξ(x) = 1 + χ y<1/2 , smaller in the upper part of the domain. Right: The car has a distinct radius of curvature on the right and left.

Figure 21 :

 21 Figure 21: Left: Level set of the cost function, depending on the physical position only c = c(x), and which is small in the neighborhood of two curves, with respectively low torsion and low curvature. Center: The minimal path backtracked with the Reeds-Shepp model on R 3 × S 2 follows the low curvature curve. Right: The minimal path backtracked with the Reeds-Shepp dual model follows the low torsion curve. A point γ(t) = (x(t), n(t)) of the backtracked path is represented by an arrow with origin x(t) and direction n(t).

Figure 22 :

 22 Figure 22: Left: Motion planning with the three-dimensional Reeds-Shepp model, with configuration space R 3 ×S 2 . Note the cusps, where the vehicle orientation is reversed. Center: Likewise, but with the Reeds-Shepp forward variant. Right: orientations for the forward model, note the in place rotations at some of the path endpoints.

Figure 23 :

 23 Figure 23: Left: Consider a seed point p in the configuration space M := R 2 × S 1 , and a family of tips q 1 , • • • , q n ∈ M. According to the Petitot Citti-Sarti model, the visual cortex V 1 will infer a connection between the seed and the closest tip w.r.t. the sub-Riemannian Reeds-Shepp model (with parameter ξ to be determined). Center and right: This principle explains two visual illusions due to Poggendorff, where the two parts of an occluded straight line or circle appear to mismatch.

It uses a semi-Lagrangian discretization, in contrast with the Eulerian discretization chosen for the HFM library, see §1.3.

In the HFM software, two distinct scales hx and h θ are used, and the latter is specified via the integer n θ := 2π/h θ .

We rely on basis reduction techniques to compute the tensor decomposition of Proposition 2.2, see §B. Strictly speaking, their cost depends on the tensor condition number, but it grows so slowly that it can be regarded as constant for the applications of interest, see[START_REF] Phong | Low-dimensional lattice basis reduction revisited[END_REF] for a complexity analysis of similar methods.

Using the fact that points are tagged as Accepted sequentially, in the order of increasing values of U .

This notation is slightly abusive here, since there is no hidden proportionality constant.

The C++ code is templated over these types.

A Byte, the unit of computer memory capacity, consists of 8 bits.

The three floating point coefficients of this quadratic polynomial are stored in an array sized according to the number ≈ N 1-1 d of trial points. A table of unsigned integers stores, for each point p ∈ X ∪ ∂X, the corresponding array index if p is tagged Trial, or a dummy value otherwise.

The ODE is in fact solved backwards in time, but we omit this detail here so as to alleviate notations.

Such inconsistent values typically arise when the front arrival times u are discontinuous in the domain interior Ω, with e.g. the Dubins model.

L = 1 for isotropic or curvature penalized models, but L = d for diagonal models, see §2.2.

github.com/Mirebeau/HFM_Python_Notebooks

github.com/Mirebeau/HFM_Mathematica_Notebooks

github.com/Mirebeau/HamiltonFastMarching

It turns out empirically that the LBFGS algorithm still works with the classical non-differentiable scheme (19), although it requires twice as many iterations to reach the same accuracy.

Fast-Marching using Lattice Basis Reduction, and Fast-Marching using Voronoi's First Reduction

The current numerical implementations is indeed limited to dimensions 2 and 3.

Assuming that the cost function c is constant.

The formula numerically implemented in the HFM library is intermediate between the original singular model (ε = 0) and the relaxed model (with the given ε > 0), in the spirit of Proposition 2.3.

On a 2.7GHz processor using a single core.

The experiment uses a partial map of Museum Georges Pompidou in Paris

We need not consider the case ẋ = 0, since F was extended to the line {(0, θ); θ ∈ R} by its lower semi continuous envelope.

This work was partly supported by ANR research grant MAGA, ANR-16-CE40-0014

Consider δ = (δ 1 , • • • , δ n ) ∈ R n , and assume that H(δ) = 1 and that H is differentiable at δ. Then denoting

one has F( v) ≤ 1. In addition, equality holds if there exists v ∈ E * such that δ i = v, ėi for all 1 ≤ i ≤ I.

Before turning to the proof, we make the connection with Proposition C.1, which assumptions are scattered over the paper. The assumptions of positivity, convexity, and two homogeneity of the function H, and the equality (52) (left) follow from the general expression (13) (understood, obviously, with an = sign, instead of the informal ≈). The element δ ∈ R n gathers the finite differences δ i = U (p) -U (p -h ėi ) of the map U in Proposition C.1, and the identity H(δ) = 1 expresses that U is a solution to the discretized PDE (33) (left). Finally, the vector v defined by (53) generalizes the expression (39) which is limited to Hamiltonians of the specific form (32).

Proof of Proposition C.1. The result follows from a rather direct computation, presented below. Before that, we observe that (52, right, rhs) is always non-negative, by choosing p = 0 and observing that H(0) = 0 by homogeneity, hence the function F is well defined. Let p ∈ E * be arbitrary, and let α = ( p, ėi )

= H(δ) = 1.

We use in each line successively: (i) definitions (52, left) and (53), of v and H, (ii) the definition of α ∈ R n , (iii) Euler's identity dH(δ), δ = 2H(δ) which follows from the 2-homogeneity of H, (iv) the convexity inequality dH(δ), α -δ ≤ H(α) -H(δ), (v) the assumption H(δ) = 1. Taking the supremum of p ∈ E * we obtain F( v) 2 ≤ 1, thus F( v) ≤ 1 as announced. Finally, under the additional assumption, one can test against the co-vector p := v, so that α = δ. In that case (54) becomes an equality and therefore F * ( v) 2 ≥ 1, thus F * ( v) = 1 as announced.

D Axes ordering and pixel area

Software visualization packages have introduced an unfortunately large number of incompatible ways to display a two or three dimensional array of values. In this appendix, we discuss the conventions chosen in the HFM library, and how to interface with other common tools.

We assume that the HFM library is provided with the following parameters: