Jean-Marie Mirebeau

Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction *

We introduce a modification of the Fast Marching algorithm, which solves the anisotropic eikonal equation associated to an arbitrary continuous Riemannian metric M, on a two or three dimensional domain. The algorithm has a complexity O(N ln N +N ln κ(M)), where N is the discrete domain cardinality. The logarithmic dependency in the maximum anisotropy ratio κ(M) of the Riemannian metric allows to handle extreme anisotropies for a limited numerical cost. We prove the consistence of the algorithm, and illustrate its efficiency by numerical experiments. The algorithm relies on the computation at each grid point z of a special system of coordinates: a reduced basis of the lattice ZZ m , with respect to the symmetric positive definite matrix M(z) encoding the desired anisotropy at this point.

Introduction

The anisotropic Eikonal equation, or static Hamilton-Jacobi equation, is a Partial Differential Equation (PDE) which describes an elementary front propagation model: the speed of the front depends only on the front position and orientation. This PDE is encountered in numerous applications, such as motion planning control problems [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF], modeling of bio-medical phenomena [START_REF] Sermesant | An anisotropic multi-front fast marching method for real-time simulation of cardiac electrophysiology[END_REF], and image analysis [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF]. It was also recently used in the context of medical image analysis [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] for extracting vessels in two dimensional projections or three dimensional scans of the human body, and for performing virtual endoscopies. This application requires to solve a highly anisotropic generalized eikonal equation with a high resolution on a cartesian grid, at a computational cost compatible with user interaction. It is one of our key motivations.

This paper is devoted to the construction and the study of a new algorithm, Fast Marching using Lattice Basis Reduction (FM-LBR), designed to solve the anisotropic eikonal equation associated to a given Riemannian metric M, and able to handle large or even extreme anisotropies. The domain must be of dimension two or three, and discretized on a cartesian grid. The FM-LBR, as its name indicates, is a variant of the classical Fast Marching algorithm [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], an efficient method for solving the eikonal equation when the metric is isotropic (proportional at each point to the identity matrix). Lattice Basis Reduction [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] is a concept from discrete mathematics, used in the FM-LBR to produce sparse causal stencils for the discretization of the eikonal equation; it allows to benefit in an optimal way of the interplay between the Riemannian geometric structure of the PDE, and the arithmetic structure of the discretization grid. A similar technique is used in [START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF] to construct sparse non-negative stencils for anisotropic diffusion.

In order to illustrate the specificity of our approach, we need to introduce some notation. Denote by S + m the collection of m × m symmetric positive definite matrices, and associate to each M ∈ S + m the norm u M := u, M u on IR m . Consider a bounded open domain Ω ⊂ IR m , equipped with a Riemannian metric M ∈ C 0 (Ω, S + m ). We address the anisotropic eikonal equation: find the unique viscosity [START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF] solution D : Ω → IR of ∇ D(z) M(z) -1 = 1 for almost every z ∈ Ω, D = 0 on ∂Ω.

See the end of §1 for different boundary conditions, and algorithmic restrictions on the dimension m. Introduce the Riemannian length of a Lipschitz path γ : [0, 1] → Ω: length(γ) :=

1 0 γ ′ (t) M(γ(t)) dt, (2)
and denote by D(x, y) the length of the shorted path joining x, y ∈ Ω, also referred to as the Riemannian distance between these points. The PDE (1) admits an optimal control interpretation: D(x) = min{D(x, y); y ∈ ∂Ω} is the minimal distance from x ∈ Ω to the boundary. Consider a discrete set Z ⊂ IR m , which in our case will be a cartesian grid. Discretizations of (1) take the form of a fixed point problem: find d :

Z → IR such that d(z) = Λ(d, z) for all z ∈ Z ∩ Ω, d(z) = 0 for all z ∈ Z \ Ω. (3)
This formulation involves the Hopf-Lax update operator Λ(d, z) [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF][START_REF] Kushner | Numerical Methods for Stochastic Control Problems in Continuous Time[END_REF][START_REF] Gonzales | On Deterministic Control Problems: an Approximate Procedure for the Optimal Cost, I, the Stationary Problem[END_REF], which mimics at the discrete level Belmann's optimality principle associated to the optimal control interpretation of [START_REF] Alton | Fast Marching Methods for Stationary Hamilton-Jacobi Equations with Axis-Aligned Anisotropy[END_REF]. See Appendix B for details on this principle, the approximations underlying its discretization, and their accuracy. The definition of Λ(d, x) involves a mesh (or stencil ) V (x) of a small neighborhood of x ∈ Z ∩ Ω, with vertices on Z, and reads: Λ(d, x) := min y∈∂V (x)

xy M(x) + I V (x) d(y), [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] where I V denotes piecewise linear interpolation on a mesh V . (In this paper, a mesh in IR m is a finite collection T of m-dimensional non-flat simplices, which is conforming in the sense that the intersection T ∩ T ′ of any T, T ′ ∈ T is the convex hull of their common vertices.) Numerical solvers of the eikonal equation differ by (i) the construction of the stencils V (z), z ∈ Ω, and (ii) the approach used to solve the system (3), which is inspired by the algorithms of Bellmann-Ford or of Dijkstra used for computing distances in graphs, instead of continuous domains. The algorithm presented in this paper, FM-LBR, belongs to the category of Dijkstra inspired algorithms with static stencils, and among these is the first one to guarantee a uniform upper bound on the stencil cardinality, independently of the Riemannian metric M. The anisotropy ratio κ(M ) of a matrix M ∈ S + m , and the maximum anisotropy κ(M) of the Riemannian metric M, are defined as follows:

κ(M ) := max u = v =1 u M v M = M M -1 , κ(M) := max z∈Ω κ(M(z)). (5)
We denote by N := #(Z ∩ Ω) the cardinality of the discrete domain.

• Bellman-Ford inspired algorithms. The discrete fixed point problem (3) is solved via Gauss-Seidel iteration: the replacement rule d(z k ) ← Λ(d, z k ) is applied for k = 0, 1, 2, ... to a mutable map δ : Z ∩ Ω → IR, until a convergence criterion is met. In the fast sweeping methods, see [START_REF] Tsai | Fast sweeping algorithms for a class of Hamilton-Jacobi equations[END_REF] and references therein, the sequence of points (z k ) k≥0 enumerates repeatedly the lines and the columns of Z ∩ Ω. Alternatively this sequence is obtained via a priority queue in the Adaptive Gauss-Seidel Iteration (AGSI) of Bornemann and Rasch [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF]. The stencil V (z) of a point z ∈ Z ∩ Ω is usually the offset by z of a fixed stencil V given at the origin, such as those illustrated on Figure 1.

Fast sweeping methods have O(λ(M)N ) complexity when the metric M is isotropic (proportional to the identity at each point), but this result does not extend to anisotropic Riemannian metrics, see [START_REF] Zhao | A Fast Sweeping Method for Eikonal Equations[END_REF] for the proof and the expression of λ(M). The AGSI has complexity O(µ(M)N 1+ 1 m ), for arbitrary anisotropic Riemannian metrics, where µ(M) is a non explicit constant which depends on global geometrical features of the metric [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF]. The AGSI is a popular, simple, and quite efficient method, which is included for comparison in our numerical tests.

• Dijkstra inspired algorithms. The system (3) is solved in a single pass, non-iteratively, using an ordering of Z ∩ Ω determined at run-time. This is possible provided the Hopf-Lax update operator satisfies the so-called "causality property", see Proposition 1.1, which can be ensured if the stencil V (z) of each z ∈ Z ∩ Ω satisfies some geometrical properties depending on M(z), see Definition 1.2. The different Dijkstra inspired methods are characterized by the construction of the stencils V (z), in contrast with Bellman-Ford inspired methods which are characterized by the choice of the sequence (z k ) k≥0 . Solving the system (3) with a Dijkstra inspired algorithm has complexity O(µ(M)N ln N ), where µ(M) is an upper bound on the cardinality of the stencils (the number of simplices they are built of).

In the Ordered Upwind Method (OUM) of Sethian and Vladimirsky [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Vladimirsky | Static PDEs for Time-Dependent Control Problems[END_REF], the stencils are constructed at run-time; their cardinality is bounded by O(κ(M) m ) and drops to O(κ(M) m-1 ) as N → ∞. In contrast, the stencils are constructed during a preprocessing step and then static in the Monotone Acceptance Ordered Upwind Method (MAOUM) of Alton and Mitchell [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF]; their cardinality is bounded by O(κ(M) m ). The FM-LBR introduced in the present work uses a similar approach, except that the stencils cardinality is O(1), fully independent of the Riemannian metric M. The complexity estimates are thus O(κ(M) m N ln N ) for the OUM and the MAOUM (asymptotically O(κ(M) m-1 N ln N ) for the OUM), and O(N ln N + N ln κ(M)) for our approach the FM-LBR, where the second term in the complexity accounts for the stencil construction.

The above mentioned algorithms are consistent for the anisotropic eikonal equation associated to an arbitrary continuous Riemannian metric M ∈ C 0 (Ω, S + m ), in the sense that the discrete output d h of the algorithm executed on the grid Z h := hZZ 2 , of scale h > 0, converges to the viscosity solution D of the continuous problem (1) as h → 0. Some more specialized variants of the fast marching algorithm are only consistent for a restricted set of metrics, but can be executed nonetheless with an arbitrary anisotropic metric M; in that case the discrete system (3) may not be solved, and the numerical results are variable, see §3. For instance the original fast marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] is consistent if M(z) is proportional to the identity matrix for each z ∈ Ω, and more generally if M(z) is a diagonal matrix. In addition to these cases of isotropy and axis-aligned anisotropy, some variants are also consistent if the metric anisotropy κ(M) is smaller than a given bound κ 0 , see [START_REF] Jbabdi | Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography[END_REF] and Figure 1. Our numerical experiments §3 include for comparison one of these methods: Fast Marching using the 8 point stencil (FM-8, center left stencil on Figure 1), which is popular in applications [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] thanks to its short computation times and despite the lack of convergence guarantee for arbitrary metrics. Depending on the implementation [START_REF] Rasch | Remarks on the O(N ) Implementation of the Fast Marching Method[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], involving either a sorted list or a bucket sort, these methods have complexity O(N ln N ) or O(Υ(M)N ), where

Υ(M) := max z∈Ω M(z) max z ′ ∈Ω M(z ′ ) -1 .
In the applications for which our method is intended, one typically has ln(N ) < ∼ κ(M) ≤ Υ(M) ≪ N , in such way that the complexity O(N ln N + N ln κ(M)) of the proposed method is comparable to O(N ln N ) and smaller than O(Υ(M)N ). In summary, the FM-LBR combines the universal consistency (i.e. for any Riemannian metric) of the AGSI, OUM and MAOUM, with a quasi-linear complexity just as the original Fast Marching algorithm.

Remark. Each solver of the eikonal equation comes with a specific construction of the stencils V (z), z ∈ Ω. The highly efficient, but specialized, approach used in the FM-LBR limits its potential for generalization, see the end of §1. Static adaptive stencils also have a memory impact, see Remark 1.10. Since the Hopf-Lax update operator (4) depends on the stencils, the discrete solution d of (3) is scheme dependent, and so is its accuracy. See Appendix B for a heuristic accuracy analysis, and [START_REF] Mirebeau | On the Accuracy of Anisotropic Fast Marching[END_REF] for the case of constant metrics. Numerical experiments §3, on application inspired test cases, show that the FM-LBR accuracy is competitive.

1 Sparse causal stencils for the anisotropic eikonal equation

Our main contribution is the construction of discretization stencils for anisotropic eikonal equations, which have a uniformly bounded cardinality and preserve a structural property of the PDE: causality, inherited from its interpretation as a deterministic control problem. Discrete causality, allowing to solve the fixed point system (3) in a single pass, is the following property.

Proposition 1.1 (Causality property, J.A. Sethian, A. Vladimirsky, Appendix of [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF]). Let x ∈ IR m , and let V be a finite mesh of a neighborhood of x. Let M ∈ S + m and let us assume that yx, M (zx) ≥ 0 for any vertices y, z of a common face of ∂V (acuteness condition). Consider a discrete map d defined on the vertices of V , and the optimization problem

Λ := min y∈∂V x -y M + I V d(y),
where I denotes piecewise linear interpolation. Then the minimum defining Λ is attained on a

k-face [y 1 , • • • , y k ] of ∂V , with k ≤ m, such that Λ > d(y i ) for all 1 ≤ i ≤ k.
A mesh V satisfying the geometric acuteness condition of Proposition 1.1, is called a causal stencil, at the point x, and with respect to M .

Consider the directed graph G, associated to the fixed point problem [START_REF] Barth | Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains[END_REF], where for all x, y ∈ Z we place an arrow x → y iff the minimum defining d(x) = Λ(d, x) is attained on a face of the stencil ∂V (x) containing y. The causality property, states that the presence of an arrow x → y implies d(x) > d(y); in particular the graph G has no cycles, hence the system of equations (3) does not feature any dependency loop. The Fast-Marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] traces back theses dependencies, determining at run-time an ordering (x i ) N i=1 , of the discrete domain Z ∩ Ω, such that the distances (d(x i )) N i=1 are increasing. We reproduce this method for completeness, see Algorithm 1, but refer to its original introduction [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], or its use with alternative static adaptive stencils in the MAOUM [START_REF] Alton | Fast Marching Methods for Stationary Hamilton-Jacobi Equations with Axis-Aligned Anisotropy[END_REF], for the proof that it solves the discrete system (3). We denoted by Λ(d, x; b, y) the modification of the Hopf-Lax update operator (4) in which the minimum is only taken over faces (of any dimension) of ∂V (x) which vertices (i) contain y, and (ii) are all accepted. Regarding the FM-LBR complexity O(N ln N + N ln κ(M)), we refer for details to the classical analysis in [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Alton | Fast Marching Methods for Stationary Hamilton-Jacobi Equations with Axis-Aligned Anisotropy[END_REF] and simply point out that (i) each FM-LBR causal stencil costs O(ln κ(M)) to construct, (ii) maintaining a list of Ω∩Z, sorted by increasing values of the mutable map d, costs O(ln N ) for each modification of a single value of d, with a proper heap sort implementation, and (iii) the optimization problem defining the Hopf-Lax update (4), or its variant Λ(d, x; b, y), has an explicit solution: the minimum associated to each face of ∂V (x) is the root of a simple univariate quadratic polynomial, see Appendix of [START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF]. Memory usage is discussed in detail in Remark 1.10.

As announced, we limit our attention to PDE discretizations on cartesian grids, of the form

Z = hR(ξ + ZZ m ), (6)
where h > 0 is a scaling parameter, R a rotation, and ξ ∈ IR m an offset (in ( 6) and ( 7) we abusively apply geometric transformations not only to points, but also to sets of points and meshes). The use of an unbounded grid (6) is a mathematical artifact aimed to simplify the exposition; only points of Z close to Ω (precisely: such that V [y] = ∅) play an active role in Algorithm 1. We may construct a causal stencil at x ∈ Z with respect to M ∈ S + m , by suitably scaling, translating and rotating an R T M R-acute mesh T , defined below:

V = x + h R T . ( 7
)
Definition 1.

2. An M -acute mesh, where M ∈ S + m , is an m-dimensional mesh T covering a neighborhood of the origin, and such that the vertices v 0 ,

• • • , v m of any simplex T ∈ T satisfy (i) v 0 = 0, (ii) (v 1 , • • • , v m ) form a basis of ZZ m , and (iii) v i , M v j ≥ 0 for all 1 ≤ i < j ≤ m. The condition | det(v 1 , • • • , v m )| = 1 ensures that 0 (resp.
x) is the only grid point in the interior of the domain covered by T (resp. V ), hence that information does not "fly over" some grid points when solving [START_REF] Barth | Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains[END_REF]. Applying the next proposition to the classical stencils of Figure 1, we obtain that they are causal for Riemannian metrics of limited anisotropy. Proposition 1.3. Let T be an m-dimensional mesh satisfying the requirements of Definition 1.2, except (iii) (which does not make sense without a given matrix M ∈ S + m ). Let

κ(T ) := 1 + γ(T ) 1 -γ(T )
, where γ(T ) := min

T,(u,v) u, v u v , (8)
and where the minimum in γ(T ) is taken among all non-zero vertices u, v of a common simplex T ∈ T . The mesh T is M -acute for any M ∈ S + m such that κ(M ) ≤ κ(T ).

Proof. Let u, v be two non-zero vertices of a common simplex T ∈ T , and let M ∈ S + m . Let u ′ := u/ u and let v ′ := v/ v . By construction one has

u ′ + v ′ 2 = 2(1 + u ′ , v ′ ) ≥ 2(1 + γ(T )), u ′ -v ′ 2 = 2(1 -u ′ , v ′ ) ≤ 2(1 -γ(T )).
Let us assume for contradiction that u, M v < 0, which implies that

u ′ + v ′ M < u ′ -v ′ M . Observing that κ(M ) 2 = M M -1 ≥ u ′ -v ′ 2 M u ′ -v ′ 2 u ′ + v ′ 2 u ′ + v ′ 2 M > 1 + γ(T ) 1 -γ(T ) ,
we obtain that κ(M ) > κ(T ), which concludes the proof of this proposition.

Our construction of M -acute meshes relies on special coordinate systems in the grid ZZ m , adapted to the anisotropic geometry encoded by M ∈ S + m .

Definition 1.4 (Bases and superbases). A basis of ZZ

m is an m-plet (e 1 , • • • , e m ) ∈ (ZZ m ) m such that | det(e 1 , • • • , e m )| = 1. A superbase of ZZ m is an m + 1-plet (e 0 , • • • , e m ) such that e 0 + • • • + e m = 0, and (e 1 , • • • , e m ) is a basis of ZZ m . Definition 1.5. A superbase (e 0 , • • • , e m ) of ZZ m is said to be M -obtuse iff e i , M e j ≤ 0 for all 0 ≤ i < j ≤ m.
There exists for each M ∈ S + m , m ∈ {2, 3} at least one M -obtuse superbase [START_REF] Conway | Low-dimensional lattices. VI. Voronoi reduction of threedimensional lattices[END_REF]. The construction of 2D and 3D M -obtuse superbases relies on lattice basis reduction algorithms [START_REF] Lagrange | Recherches d'arithmétique[END_REF][START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] (hence the name of our numerical scheme) and has cost O(ln κ(M )), see §2. We arrive at the main contribution of this paper: the FM-LBR stencils, which are causal and of bounded cardinality.

Proposition 1.6 (The FM-LBR acute meshes). Let M ∈ S + m , and let (e 0 , • • • , e m ) be an Mobtuse superbase of ZZ m (if one exists). An M -acute mesh T (M ) is obtained by collecting the simplices of vertices

( k-1 i=0 b i ) m k=0 associated to all (m + 1)! permutations (b i ) m i=0 of (e i ) m i=0 . It has 2 m+1 -1 vertices.
Proof. Proof of the properties of the simplices. Let (b i ) m i=0 be a permutation of (e i ) m i=0 , and

let v k := k-1 i=0 b i , for all 0 ≤ k ≤ m. Clearly (i) v 0 = 0, and (ii) | det(v 1 , • • • , v m )| =
| det(b 1 , • • • , b m )| = | det(e 1 , • • • , e m )| = 1. Acuteness condition (iii): for any 1 ≤ k < l ≤ m one has v l = -m j=l b j by Definition 1.4, hence by Definition 1.5 v k , M v l = - 0≤i<k l≤j≤m e i , M e j ≥ 0.
Proof that T (M ) is a conforming mesh covering a neighborhood of the origin. Consider the m + 1-dimensional Kuhn simplices

T σ := {(λ 0 , • • • , λ m ); 0 ≤ λ σ(0) ≤ λ σ(1) ≤ • • • ≤ λ σ(m) ≤ 1}, associated to all permutations σ of {0, • • • , m}, which form a conforming mesh T 0 of [0, 1] m+1 . The linear map A : IR m+1 → IR m , defined by A(λ 0 , • • • , λ m ) := -m i=0 λ i e i , has a kernel generated by v 1 := (1, • • • , 1). It sends [0, 1] m+1 onto a neighborhood of 0 ∈ IR m ,
and transforms T 0 into an m-dimensional conforming mesh T M of this neighborhood, collapsing onto 0 the common edge [0, v 1 ] of the simplices T σ . Other vertices of [0, 1] m+1 have pairwise distinct images by A, since their difference is not proportional to v 1 ; hence T M has 2 m+1 -1 vertices (one less than T 0 ). Noticing the identity

-m i=0 λ i e i = m k=1 (λ σ(k) -λ σ(k-1) ) k-1 i=0 e σ(k)
, we find that the image of T σ by A is the simplex of T (M ) associated to the permuted superbase (e σ(i) ) m i=0 , hence T M = T (M ), which concludes the proof.

Obtuse superbases are similarly used in [START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF] to produce sparse non-negative stencils for Anisotropic Diffusion (AD-LBR scheme), with different results: 3D stencils have 12 non-zero vertices in [START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF], and 14 here. This illustrates the versatility of this concept, which can be used to design stencils satisfying various geometric properties: an acuteness condition for the FM-LBR (implying the scheme causality), and the non-negative decomposition of a tensor for the AD-LBR (guaranteeing the scheme monotonicity). Definition 1.7. A family of meshes (T (x)) x∈Ω , defined for all points of an open domain Ω ⊂ IR m is admissible iff there exists two constants 0 < r and R < ∞ such that: for each x ∈ Ω

• The mesh T (x) covers neighborhood of 0, and its vertices belong to ZZ m .

• (Boundedness)

The vertices e of T (x) satisfy e ≤ R.

• (Stability) There exists a basis B(x) of ZZ m which elements, and their opposites, are vertices of T (y) for all y ∈ Ω such that xy < r.

First Last v u v u v u v u v u v u v u v u v u v u v u v u v u First Last
Figure 3: The unit sphere {u ∈ IR 2 ; u M = 1}, an M -reduced basis (u, v), and the boundary of the FM-LBR mesh T (M ), for some M ∈ S + 2 of anisotropy ratio κ(M ) = 6, and eigenvector (cos θ, sin θ), θ ∈ [π/4, π/2], associated to the small eigenvalue.

A fortunate but non-trivial fact is that the FM-LBR family of meshes (T (M(x))) x∈Ω is admissible, see Proposition 1.8 and Figures 3 and4. This property implies the FM-LBR consistence, using Proposition 1.9 which is a minor extension of the convergence results in [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF].

∈ C 0 (Ω, S + m ). If m ≤ 3 then the FM-LBR family of meshes (T (M(x))) x∈Ω is admissible (with Boundedness constant R = C m κ(M), C 2 = 2, C 3 = 4). More generally, if m ≤ 4 then any family of meshes (T (x)) x∈Ω such that T (x) is M(x)-acute for all x ∈ Ω,
lim h→0 max z∈Z h ∩Ω | d h (z) -D(z)| = 0. ( 9
)
To sum up the FM-LBR strengths, this algorithm is universally consistent, has a competitive accuracy w.r.t. alternative methods, and its computational cost is almost unaffected by the Riemannian metric anisotropy. For fairness we discuss below the potential downsides of our original stencil construction, which limits the potential for generalization (efficiency is as often at the cost of specialization), and impacts memory usage.

• (Finsler metrics) The FM-LBR only applies to the anisotropic eikonal equation associated to a Riemannian metric, while other methods such as the AGSI, OUM, MAOUM [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF] can handle more general Finsler metrics. Indeed the structure • M(z) of the local norm at z associated to a Riemannian metric is required in the FM-LBR stencil construction Proposition 1.6, which involves an M(z)-obtuse superbase. Finsler metrics are in contrast defined by arbitrary asymmetric norms | • | z , depending continuously on z.

Constructing causal static stencils is an active subject of research in the case of Finsler metrics. A characterization obtained in [START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF], was used in [START_REF] Mirebeau | Efficient Fast Marching with Finsler Metrics[END_REF] to develop the FM-ASR (Fast Marching using Anisotropic Stencil Refinement), which is close in spirit and in efficiency to the FM-LBR, but has a different application range, since it handles Finsler metrics on two dimensional domains. It was also observed in [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF][START_REF] Alton | Fast Marching Methods for Stationary Hamilton-Jacobi Equations with Axis-Aligned Anisotropy[END_REF] that the canonical stencil (Figure 1, left) is causal for Finsler metrics featuring only axis-aligned anisotropy, in such way that the original fast marching algorithm can be applied.

• (Domain discretization) The FM-LBR requires a cartesian grid discretization. In contrast an important research effort [START_REF] Barth | Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains[END_REF][START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF][START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF][START_REF] Li | Numerical schemes for Hamilton-Jacobi equations on unstructured meshes[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF] has been devoted to the more difficult setting of meshed domains, with an unstructured set of vertices. These methods are natural candidates for the computation of distances on a manifold S ⊂ IR m , while the FM-LBR would require a (collection of) local chart(s) ϕ : Ω → S equipped with the metric M(z) := ∇ϕ(z) T ∇ϕ(z). The FM-LBR heavily relies on the cartesian grid arithmetic structure, through the concept of obtuse superbases.

• (Boundary conditions) The null boundary conditions chosen in the eikonal equation ( 1) can be replaced with Dirichlet data D 0 : ∂Ω → IR, of 1-Lipschitz regularity [START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF] with respect to the Riemannian distance D(•, •), see [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF]. In that event, Algorithm 1 requires to extend the boundary data D 0 to a ghost layer, containing all grid points y ∈ Z \ Ω such that the reverse stencil V [y] is non-empty. The FM-LBR uses large stencils, of euclidean radius O(hκ(M)) on a grid (6) of scale h, which complicates this extension in contrast with e.g. the AGSI [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF] of stencil radius O(h).

Outflow boundary conditions, on a portion Γ ∂Ω of the domain's boundary, are natural in applications, see §3. They are implemented by excluding in the definition (4) of Λ(d, x), faces of ∂V (x) containing an exterior vertex y ∈ Z \ Ω close to Γ. If x lies in a corner of Ω, and if the stencil V (x) is strongly anisotropic, then it may happen that all the vertices of ∂V (x) lie outside Ω, so that the solution of (3) satisfies d(x) = +∞. (In our experiments §3, this happened in the square domain's corners when test case 1 was rotated by an angle θ ∈ [0.56, 0.61] radians. These four infinite values were rejected when estimating numerical errors.) Despite these minor inconveniences, the FM-LBR behaves remarkably well in our numerical experiments, §3 and Appendix B, with these more general boundary conditions.

• (Dimension) The FM-LBR causal stencils construction, see Proposition 1.6, is limited to domains of dimension 2 and 3, because some matrices M ∈ S + 4 do not admit any M -obtuse superbase [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF]. An alternative construction of M -acute meshes of uniformly bounded cardinality is proposed in [START_REF] Mirebeau | On the Accuracy of Anisotropic Fast Marching[END_REF], for all M ∈ S + 4 ; this cardinality is not small unfortunately, with 768 simplices.

The FM-LBR however extends in a straightforward manner to Riemannian metrics M having a block diagonal structure, with blocks of size 1, 2 or 3, using the following construction. For i ∈ {1, 2} let m i be a positive integer, let M i ∈ S + m i , and let T i be an M i -acute mesh. Let m := m 1 + m 2 and let M ∈ S + m be the matrix of diagonal blocks M 1 , M 2 . An M -acute mesh T is obtained by collecting the m-dimensional simplices of vertices (0, 0), (u

1 , 0), • • • , (u m 1 , 0), (0, v 1 ), • • • , (0, v m 2 ), where the simplices of vertices 0, u 1 , • • • , u m 1 and 0, v 1 , • • • , v m 2 belong to T 1 and T 2 respectively.
Block diagonal metrics are not uncommon in the context of medical imaging [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], as they inherit the cartesian product structure of the fast marching domain: Ω = Ω 0 × Ω 1 , where Ω 0 is a physical domain of dimension ≤ 3, and Ω 1 is an abstract parameter domain of dimension ≤ 2.

Outline. Further insight on the FM-LBR is given in §2, including the proof of Propositions 1.8 and 1.9. Numerical experiments are presented in §3. In addition, Appendix A describes Point (II) can be avoided for the AGSI and the OUM when these methods are executed on a mesh with a trivial periodic structure, which is the case in our experiments. For the FM-LBR, Point (II) amounts to storing the non-empty reverse stencils V [y], at all points of Y := {y ∈ Z; V [y] = ∅}, since the direct ones can be recomputed individually on demand for a minor cost. The set Y the union of Ω ∩ Z and of a thin boundary layer (in our experiments, Y = Ω ∩ Z due to the use of outflow boundary conditions). The chosen data structure uses two tables: one of vectors (the differences xy, for x ∈ V [y], y ∈ Y , enumerated consecutively), and one of #(Y ) ≈ N integers (the start, for each y ∈ Y , of the description of V [y] in the previous table). We represent integers in 32bit format, and vector components in 8bit format, since these are small integers by construction.

Summing up, we find that the memory requirements of the FM-LBR are larger than those of the AGSI or the OUM (on a grid), by a factor ranging from 2 (metric and solution stored in double precision), to 9 (analytical metric, solution stored in single precision), through 3 (metric and solution stored in single precision), in two dimensions. Respectively, in three dimensions, from 2.6 to 24, through 4.3.

Analysis of the FM-LBR

We introduce in §2.1 the concepts of Lattice Basis Reduction. They are used in §2.2 to estimate the construction cost of the FM-LBR meshes T (M ), and to prove their admissibility in the sense of Definition 1.7, as announced in Proposition 1.8. We finally prove in §2.3 the announced convergence result Proposition 1.9.

Introduction to Lattice Basis Reduction

We briefly introduce the framework of (low dimensional) Lattice Basis Reduction, used in the next subsection to construct and demonstrate the properties of M -acute meshes. See [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] and references therein for more details on this rich theory, from which we use only one result: Theorem 2.2 stated below. We denote by b

1 ZZ + • • • + b k ZZ the sub-lattice of ZZ m generated by b 1 , • • • , b k ∈ ZZ m : b 1 ZZ + • • • + b k ZZ := {b 1 z 1 + • • • + b k z k ; z 1 , • • • , z k ∈ ZZ}.
If k = 0 then the above sum equals {0} by convention.

Definition 2.1. Let 1 ≤ m ≤ 4 and let M ∈ S + m . A basis (b 1 , • • • , b m ) of ZZ m is said M -reduced iff for all 1 ≤ k ≤ m: b k ∈ argmin{ z M ; z ∈ ZZ m \ (b 1 ZZ + • • • + b k-1 ZZ)}. ( 10
)
If M ∈ S + m is a diagonal matrix of coefficients (λ 1 , • • • , λ m ), and if 0 < λ σ(1) ≤ • • • ≤ λ σ(m)
for some permutation σ, then the permutation (e σ(i) ) m i=1 of the canonical basis of ZZ m is an M -reduced basis of ZZ m . See Figures 2,3,4, for some examples of M -reduced bases associated to non-diagonal matrices M . In dimension m ≥ 5 there exists matrices M ∈ S + m such that no basis of ZZ m satisfies the relations [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF], see [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] (these relations state that b i M equals the i-th Minkowski's minimum λ i (M )). Minkowski's reduction [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] is the natural generalization of Definition 2.1 in dimension m ≥ 5. Proof. The proof is contained in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF], and we only point out here the precise reference within the paper and the slight differences in notations. The algorithm described in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] takes as input a basis (b 1 , • • • , b m ) (here: the canonical basis of IR m ) of a lattice L (here: ZZ m ), and its Gram matrix with respect to some scalar product (here: the Gram matrix is M ). The algorithm outputs a greedy reduced basis of the lattice L, a notion which coincides with Minkowski's reduction if m ≤ 4 (Lemma 4.3.2 in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF]), which itself coincides with Definition 2.1 if m ≤ 4.

The main loop of the iterative algorithm is executed at most the following number of times (Theorem 6.0.5 in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF]):

O 1 + ln max 1≤i≤m b i M -ln min u∈L u M , hence O(1 + ln M 1 2 -ln M -1 -1 2 ) = O(1 + ln κ(M ))
times in our setting. The complexity of each of these iterations is dominated by a closest vector search, described in Theorem 5.0.4 in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF], which consists of the inversion of a k × k Gram matrix, where 1 ≤ k ≤ m -1, and an exhaustive search among O(1) candidate vectors. In terms of elementary operations (+, -, ×, /) among reals, each iteration of this algorithm thus has cost O(1), and the overall cost is the number of iterations O(1 + ln κ(M )).

Note that an important part of the discussion in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF] is devoted to the special case where the vectors (b 1 , • • • , b m ) have large integer coefficients, the Gram matrix is computed with respect to the standard euclidean scalar product, and the complexity of an elementary operation (+, -, ×, /) among integers is not O(1) but depends on the size of these integers. This more subtle notion of complexity, named bit complexity, is not relevant in our setting.

The two dimensional version of the algorithms mentioned in Theorem 2.2 dates back to Lagrange [START_REF] Lagrange | Recherches d'arithmétique[END_REF], and mimicks the search for the greatest common divisor of two integers. This algorithm uses only a pair (u, v) of (mutable) variables in ZZ 2 , initialized as the canonical basis of IR 2 . The pair (u, v) becomes an M -reduced basis at the end of the following loop, which takes at most O(ln κ(M )) iterations. Round denotes rounding to a closest integer.

Do (u, v) ← (v, u -Round( u, M v / v 2 M ) v), While u M > v M .
We end this subsection with a basic estimate of the norms and scalar products of the elements of M -reduced bases.

Proposition 2.3. Let 1 ≤ m ≤ 4, let M ∈ S + m and let (b 1 , • • • , b m ) be an M -reduced basis of ZZ m . Then for any 1 ≤ i ≤ m b i ≤ κ(M ),
and

b i M ≤ κ(M ) b 1 M . ( 11
)
For any integer combination

z ∈ b 1 ZZ + • • • + b i-1 ZZ + b i+1 ZZ + • • • + b m ZZ, of the basis elements distinct from b i , one has 2| b i , M z | ≤ z 2 M . (12)
Proof. Proof of 11. Let (e j ) m j=1 denote the canonical basis of IR m . By a dimensionality argument there exists 1

≤ j ≤ m such that e j / ∈ b 1 ZZ + • • • + • • • b i-1 ZZ.
By construction [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF] we obtain b i M ≤ e j M ≤ M 1 2 since e j = 1, hence b i ≤ κ(M ) as announced in [START_REF] Kushner | Numerical Methods for Stochastic Control Problems in Continuous Time[END_REF]. Observing that b 1 M ≥ M -1 -

b i 2 M ≤ b i + z 2 M = b i 2 M + 2 b i , M z + z 2 M , hence -2 b i , M z ≤ z 2 M . Likewise 2 b i , M z ≤ z 2
M , which concludes the proof.

Properties of M -acute meshes

The FM-LBR stencil construction has numerical cost O(1 + ln κ(M )).

The construction of an FM-LBR mesh T (M ), M ∈ S + m , has unit cost (for fixed m) given an M -obtuse superbase, see Proposition 1.6. We give below a unit cost construction of an M -obtuse superbase given an M -reduced basis, in dimension m ∈ {2, 3}, which by Theorem 2.2 can itself be obtained at cost O(1 + ln κ(M )).

Dimension 2. Let (b 1 , b 2 ) be an M -reduced basis. Up to replacing b 2 with its opposite, we may assume that b 1 , M b 2 ≤ 0. Using [START_REF] Lagrange | Recherches d'arithmétique[END_REF]

we obtain 2| b 1 , M b 2 | ≤ b 1 2 M , hence M b 1 , -b 1 -b 2 ≤ 0 and likewise M b 2 , -b 1 -b 2 ≤ 0, so that (b 1 , b 2 , -b 1 -b 2 ) is an M -obtuse superbase.
Dimension 3: we reproduce without proof the construction of [START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF]. Let (b 1 , b 2 , b 3 ) be the elements of an M -reduced basis, permuted and signed so that

| b 1 , M b 2 | ≤ min{-b 1 , M b 3 , -b 2 , M b 3 }. An M -obtuse basis is given by (b 1 , b 2 , b 3 , -b 1 -b 2 -b 3 ) if b 1 , M b 2 ≤ 0, and (-b 1 , b 2 , b 1 + b 3 , -b 2 -b 3 ) otherwise.

Radius of the FM-LBR stencils.

Assume that an FM-LBR acute mesh T (M ) is built from an M -obtuse superbase obtained as in the previous paragraph1 from an M -reduced basis (b i ) m i=1 . Then one easily checks that its vertices have the form e = m i=1 ε i b i , where ε i ∈ {-1, 0, 1} and (b i ) m i=1 is the used M -reduced basis. Hence the vertices norm e ≤ m i=1 b i ≤ mκ(M ) obeys the announced bound. An M -reduced basis of ZZ m contains by definition small vectors with respect to the norm • M : the smallest linearly independent ones with integer coordinates. As a result, the FM-LBR meshes T (M ) have a small radius with respect to the norm • M . In constrast, these meshes can be large from the euclidean perspective. This has consequences on the FM-LBR accuracy, see Appendix B.

Stability of M -acute meshes.

Consider the following distance on the set S + m of symmetric positive definite matrices: for all

M, N ∈ S + m dist(M, N ) := sup u =0 |ln u M -ln u N | . ( 13
)
This distance allows to compare the norms of vectors multiplicatively, in contrast with the classical operator norm which is tailored for additive comparisons. Indeed denoting α := dist(M, N ) and

β := M 1 2 -N 1 2
, one has for all u ∈ IR 2 such that u = 1

e -α ≤ u M / u N ≤ e α , and -β ≤ u M -u N ≤ β.
The next lemma establishes a lower bound on the • M norm of points with integer coordinates outside of an N -reduced mesh, when the matrices M, N ∈ S + m are close enough. Lemma 2.4. Let M, N ∈ S + m , with m ≤ 4. Let u 1 , • • • , u m be an arbitrary M -reduced basis of ZZ m , and let T be an N -reduced mesh. Consider a point z ∈ ZZ m which is not a vertex of T . Then there exists 1 ≤ l ≤ m such that

z ∈ u 1 ZZ + • • • + u l ZZ and z 2 M e 4 dist(M,N ) ≥ u l 2 M + u 1 2 M .
Proof. Since T covers a neighborhood of the origin, there exists a simplex T ∈ T and a real λ > 0 such that λz ∈ T . Denoting by v 1 , • • • , v m the non-zero vertices of T , there exists therefore non-negative reals

α 1 , • • • , α m ∈ IR + such that z = α 1 v 1 + • • • + α m v m . Since (v 1 , • • • , v m )
∈ u 1 ZZ + • • • + u l ZZ. There exists 1 ≤ i ≤ k such that v i / ∈ u 1 ZZ + • • • + u l-1 ZZ, hence v i M ≥ u l M ; other vertices satisfy v j M ≥ u 1 M , 1 ≤ j ≤ k,
since they are non-zero and have integer coordinates. Therefore, denoting δ := dist(M, N )

e 4δ z 2 M ≥ e 2δ z 2 N = e 2δ   1≤i≤k α 2 i v i 2 N + 2 1≤i<j≤k α i α j v i , N v j   ≥ e 2δ 1≤i≤k α 2 i v i 2 N ≥ 1≤i≤k α 2 i v i 2 M ≥ u l 2 M +   1≤i≤k α 2 i -1   u 1 2 M .
Observing that

α 2 1 + • • • + α 2 k ≥ 2,
since z is not a vertex of T , we conclude the proof.

We prove in the next corollary that the vertices of an M -acute mesh contain the elements of an N -reduced basis, and their opposites, when the matrices M, N are sufficiently close. This result, and the compactness of Ω, immediately implies the second point of Proposition 1.8.

Convergence of the FM-LBR

We prove in this section the uniform convergence of the discrete system (3) solutions towards the anisotropic eikonal PDE (1) solution, under the assumptions of Proposition 1.9 and with its notations. Following the steps of [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], we begin with a discrete Lipschitz regularity estimate for the maps d h : Z h → IR. Lemma 2.7. There exists constants h 0 > 0 and C 0 < ∞ such that for all 0 < h ≤ h 0 and all x, y ∈ Z h one has

| d h (x) -d h (y)| ≤ C 0 x -y . ( 15
)
Proof. We prove below that | d h (x)-d h (y)| ≤ C 1 h when x-y = h, in other words when x and y are neighbors on the grid Z h := hZZ m . This immediately implies

| d h (x) -d h (y)| ≤ C 1 x -y 1 , where (λ 1 , • • • , λ m ) 1 := m i=1 |λ i |, hence also (15) with C 0 := C 1 √ m.

If x /

∈ Ω and y / ∈ Ω, then d h (x) = d h (y) = 0, and the result is proved. Up to exchanging x and y, we may therefore assume that x ∈ Ω. Let B := B(x) be the basis corresponding to the property (Stability) of the admissible family of meshes (T (x)) x∈Ω , see Definition 1.7. Each element e of B is a vertex of T (x), hence satisfies e ≤ R, see Definition 1.7 (Boundedness). We abusively regard B as an m × m matrix which columns are the basis elements, and observe that

B ≤ R √ m and | det(B)| = 1. Therefore B -1 ≤ B m-1 /| det B| ≤ C 2 := (R √ m) m-1 .
Since we assumed xy = h, we have y = x + εhe j , for some ε ∈ {-1, 1}, 1 ≤ j ≤ m, and where (e j ) m j=1 denotes the canonical basis of IR m . Hence denoting by (b i ) m i=1 the elements of the basis B, and by (α ij ) m i,j=1 the coefficients of B -1 , we obtain y = x + εh m i=1 α ij b i . Let s := m i=1 |α ij |, and let (x k ) s k=0 be a finite sequence of points of Z h such that x 0 := x, x s := y, and x k+1x k ∈ {±hb i ; 1 ≤ i ≤ m} for all 0 ≤ k < s.

By Cauchy-Schwartz's inequality one has

s/ √ m ≤ ( m i=1 α 2 ij ) 1 2 ≤ B -1 ≤ C 2 , hence x k - x ≤ hRC 2
√ m for any 0 ≤ k ≤ s. We choose the upper grid scale bound h 0 so that this constant smaller than the radius r involved in property (Stability) of Definition 1.7. If x k ∈ Ω, then by property (Stability) a grid point x l , with |k -l| ≤ 1 and 0 ≤ l ≤ s, is a vertex of the stencil V h (x). Hence inserting x l in (4) we obtain

d h (x k ) = Λ h (d h , x k ) ≤ x k -x l M(x) + d h (x l ) ≤ hRC 3 + d h (x l ), with C 3 := max{ M(z) 1 2 ; z ∈ Ω}. If x k / ∈ Ω then d h (x k ) = 0, hence obviously d h (x k ) ≤ d h (x l ).
Exchanging the roles of of k and l we obtain

| d h (x k ) -d h (x l )| ≤ hRC 3 , hence by the triangular inequality | d h (x) -d h (y)| ≤ hRC 3 s ≤ hRC 3 C 2 √ m, which concludes the proof.
The rest of the proof is only sketched, since it amounts to a minor adaptation of Theorem 11 in [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF]. The only difference lies in the presence, in [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], of a canonical interpolation operator on Ω (piecewise linear interpolation on a prescribed mesh). Denote by D h the bilinear interpolant2 of d h on the grid Z h , and observe that | D h (x) -D h (y)| ≤ KC 0 xy , for all x, y ∈ IR 2 and all 0 < h ≤ h 0 , where K is an absolute constant depending only on the interpolation scheme, and C 0 , h 0 , are the constants of Lemma 2.7. Note also that supp(D h ) ⊂ {z + e; z ∈ Ω, e ≤ h √ m}, hence by Lipschitz regularity D h is bounded uniformly independently of h, and therefore by Arzelà-Ascoli's theorem the family (D h ) 0<h≤h 0 is pre-compact. Considering an arbitrary converging sub-sequence D h(n) , with h(n) → 0 as n → ∞, one observes that the limit is supported on Ω, and applying the arguments of Theorem 11 in [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF] that it is a viscosity solution of the eikonal PDE [START_REF] Alton | Fast Marching Methods for Stationary Hamilton-Jacobi Equations with Axis-Aligned Anisotropy[END_REF]. Uniqueness of such a solution [START_REF] Lions | Generalized solutions of Hamilton-Jacobi equations[END_REF] implies the pointwise convergence D h (x) → D(x), as h → 0, for all x ∈ Ω. Finally the announced uniform convergence (9) follows from the uniform KC 0 -Lipschitz regularity of D h , for all 0 < h ≤ h 0 .

Numerical experiments

We compare numerically the FM-LBR with two popular solvers (AGSI, FM-8) of the eikonal equation, which enjoy a reputation of simplicity and efficiency in applications, and with the recent and closely related MAOUM. The Adaptive Gauss Seidel Iteration3 (AGSI) [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF] produces numerical approximations which are guaranteed to converge towards the solution of the continuous anisotropic eikonal equation as one refines the computation grid 4 , for an arbitrary continuous Riemannian metric M. Fast Marching using the 8 point stencil (FM-8, stencil illustrated on Figure 1, center left) does not offer this convergence guarantee, but has a quasi-linear complexity O(N ln N ), in contrast with the super-linear complexity O(µ(M)N 1+ 1 m ) of the AGSI. Fast Marching using Lattice Basis Reduction5 (FM-LBR) aims to offer the best of both worlds: a convergence guarantee, and fast computation times6 . We also implemented the Monotone Acceptance Ordered Upwind Method (MAOUM) [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF], a Dijkstra inspired method using static stencils, like the FM-LBR. The difference between these two methods is that the MAOUM stencil7 V (z) at a grid point z ∈ Ω ∩ Z, is isotropic, only depends on the anisotropy ratio κ(M(z)), and its boundary has cardinality O(κ(M(z)) m-1 ); in contrast the FM-LBR stencil is anisotropic, aligned with the ellipse defined by M(z), and of cardinality O(1). The MAOUM stencils were precomputed and stored in a look-up table, resulting in a complexity O(κ(M)N ln N ) for this algorithm in 2D.

We consider application inspired test cases, which violate some of the simplifying assumptions used in our convergence analysis Proposition 1.9: they involve outflow boundary conditions, non-trivial Dirichlet boundary conditions in Appendix B, and discontinuous Riemannian metrics in cases 3 and 4. Their exact technical description is given in Remark 3.1.

The first test is a distance computation on a parametrized surface, considered in [START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF]. As shown on Figure 6, the FM-LBR is the fastest in terms of CPU time8 , but is less accurate than the AGSI or the FM-8. Rotating this test case by the angle θ = π/6, and conducting the same experiment, shows a different story: the numerical errors of the AGSI and the FM-8 increase by a factor larger than 5, while the FM-LBR, unaffected, is now the most accurate method, see Figure 6. The FM-LBR cuts L ∞ and L 1 numerical errors by 40% in comparison with the AGSI and the MAOUM, and CPU time by 85%, while the FM-8 produces even larger errors.

Figure 6 shows that the FM-LBR offers the best accuracy for more grid orientations θ than its alternatives. The maximal error and averaged error with respect to θ are also in favor of the FM-LBR. The strong dependence of the AGSI and the FM-8 accuracy on the test case orientation is puzzling, and contrasts with the more consistent behavior of the MAOUM and the FM-LBR. The author's heuristic and personal interpretation of this phenomenon, which is open to debate and put into question by a reviewer, is that this test case is for θ = 0 dominated by (close to) axis-aligned anisotropy. The AGSI and the FM-8, which are based on small and fixed stencils, benefit from this special configuration; the FM-8 also works well for θ = π/4, because its stencil includes the four diagonals.

The second benchmark, discussed in [START_REF] Vladimirsky | Fast methods for static Hamilton-Jacobi Partial Differential Equations[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF], is inspired by seismic imaging. There is no bias here towards axis-aligned anisotropy. As shown on the table Figure 6 a smaller CPU time and offers a better accuracy than its alternatives. Note that one can also construct configurations in which anisotropy is not axis-aligned and the AGSI is more accurate than the FM-LBR. In some cases, the accuracy advantage of the AGSI even grows unboundedly as the anisotropy ratio of the metric tends to infinity. A heuristic analysis and prediction of this phenomenon is presented in Appendix B, where it is illustrated with a fifth test case.

The third [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF] and fourth test cases are relevant benchmarks if one's objective is to use fast marching methods for the segmentation of tubular structures, in medical images or volume data respectively. The FM-LBR reveals its full potential, and stands out as the only practical option, in this more difficult setting which involves a discontinuous and highly anisotropic metric. The Riemannian metric tensor M(z) is the identity matrix, except at points z on the neighborhood of a curve Γ, where M(z) has a small eigenvalue δ2 0 associated to an eigenvector tangent to Γ, and the other eigenvalue is 1 on the orthogonal space. The shortest path joining the point (1, -1) (resp. (0, 0, 3)) to the origin is extracted via "gradient descent on the Riemannian manifold (Ω, M)", see Appendix A for details:

γ ′ (t) = -M(γ(t)) -1 ∇ D(γ(t)). ( 16
)
By construction of the Riemannian metric M, traveling close and tangentially to the curve Γ is cheap. This is reflected by the level lines of D, and by the allure of the minimal path, see Figures 5, 7 and 10. Heuristically, this path joins the neighborhood of the curve Γ in straight line, almost orthogonally, and then follows it. The alignment of the minimal path with the direction of anisotropy, observed in this test case, is not an uncommon phenomenon. The FM-LBR presumably benefits a lot from this behavior in terms of accuracy, since its stencils typically provide a good "angular resolution" in the direction of anisotropy, see Figures 3,4,7. Since in addition the stencil radii remain rather small for most anisotropy orientations, see [START_REF] Mirebeau | On the Accuracy of Anisotropic Fast Marching[END_REF] for details, usually most updates for points in the "fast band" come from the fast band when the FM-LBR is run on these examples. When the fast band is missed however, accuracy degrades and zig-zag artifacts appear in the extracted path, see Figure 11. Third test case [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], in 2D. The different method's performance is illustrated on Figure 8, except for the MAOUM which showed a poor accuracy, presumably due to the huge stencils it generated. The CPU time/resolution curve of the AGSI shows a stronger slope than for the one pass solvers, presumably reflecting its super-linear complexity O(µ(M)N Qualitatively, the approximate geodesic has the right behavior for a resolution as low as 170 × 170 with the FM-LBR, and 1000 × 1000 with the AGSI. This is presumably never the case for the FM-8, which is not consistent here. minutes to 2.5 seconds (!). As illustrated on Figure 9, the better accuracy of the FM-LBR in this test case effectively translates into a better extraction of minimal paths. The reason for the, unrivaled, performance of the FM-LBR in this specific test case is partly elucidated in Appendix B.

Fourth test case, in 3D. CPU time was 105s for the FM-LBR, while the AGSI took 480s and failed to recover the minimal path presented on Figure 10 (center) (a straight line joining the two endpoints was obtained instead). The FM-LBR is capable of addressing a large scale (more than 10 millions grid points), strongly anisotropic (κ(M) = 50) three dimensional shortest path problem, with a good accuracy and within reasonable CPU time on standard laptop computer.

Remark 3.1. The four application inspired test cases. Outflow boundary conditions, except for D(0) = 0; hence the solution is the Riemannian distance to the origin: D(z) = D(z, 0). Numerical errors with respect to a 4000 × 4000 reference solution, bi-linearly interpolated, obtained in cases 1,2, with the AGSI, and in case 3 with the FM-ASR [START_REF] Mirebeau | Efficient Fast Marching with Finsler Metrics[END_REF].

1. (Geometry processing [START_REF] Vladimirsky | Fast methods for static Hamilton-Jacobi Partial Differential Equations[END_REF], κ(M) ≃ 5.1) Compute the Riemannian distance from the origin (0, 0, 0) on parametric surface of height map z(x, y) := (3/4) sin(3πx) sin(3πy). Riemannian metric: M(x, y) = Id +∇z(x, y)∇z(x, y) T . Coordinates (x, y) restricted to the square [-0.5, 0.5] 2 , discretized on a 292 × 292 grid, or in a second step to this square rotated by the indicated angle θ.

2. (Seismic imaging [START_REF] Vladimirsky | Fast methods for static Hamilton-Jacobi Partial Differential Equations[END_REF][START_REF] Sethian | Ordered Upwind Methods for Static Hamilton-Jacobi Equations: Theory and Algorithms[END_REF], κ(M) = 4) The metric M(x, y) has the two eigenvalues 0.8 -2 , 3. (Tubular segmentation [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], κ(M) = 100) Define the curve Γ(t) = t(cos ω 0 t, sin ω 0 t), t ∈ [0, 1]. Set M(z) = Id, except if there exists 0 ≤ t ≤ 1 and 0 ≤ r ≤ r 0 such that z = Γ(t) + r(cos ω 0 t, sin ω 0 t). In that case M(z) has the eigenvalues δ 2 0 and 1, the former with eigenvector Γ ′ (t). Parameters: ω 0 := 6π, r 0 := δ 0 := 0.01. Domain: [-1, 1] 2 , grid sizes n × n with 120 ≤ n ≤ 1200.

4. (Tubular segmentation, κ(M) = 50, 3D) Define the curve Γ(t) = (cos ω 0 t, sin ω 0 t, t), with ω 0 := (5/2)π. Set M(z) = Id, except if there exists t, λ, µ ∈ IR such that z = Γ(t) + (λ cos ω 0 t, λ sin ω 0 t, µ) and λ 2 + µ 2 ≤ (r 0 /2) 2 . In that case M(z) has the eigenvalues δ 2 0 and 1, the former with eigenvector Γ ′ (t) and the latter with multiplicity 2. Parameters:

ω 0 := (5/2)π, δ 0 = r 0 = 0.02. Domain: [-1.1, 1.1] 2 × [0, 3], grid size 200 × 200 × 272.

Conclusion

The FM-LBR, introduced in this paper, combines the Fast Marching algorithm with a concept from discrete geometry named Lattice Basis Reduction. It has the following strongpoints. (I, Convergence) The FM-LBR is consistent for the anisotropic eikonal equation associated to any continuous Riemannian metric, of arbitrary anisotropy. (II, Complexity) It has a numerical cost comparable to classical isotropic Fast Marching, independently of the problem anisotropy.

(III, Accuracy) The accuracy of the FM-LBR is competitive in general, and striking in test cases, related to tubular segmentation in medical images, where the Riemannian metric has a pronounced anisotropy close to and tangentially to a curve. These qualities come at the price of the specialization of the FM-LBR: (i) the Riemannian metric may not be replaced with a more general Finsler metric, see [START_REF] Mirebeau | Efficient Fast Marching with Finsler Metrics[END_REF] for an adaptation to this setting in 2D, (ii) the domain needs to be discretized on a cartesian grid, and (iii) of dimension 2 or 3. Hopefully these requirements are met in many applications, and future work will be devoted to the application of the proposed algorithm in the context of medical image processing. i=1 and piecewise linear path γ in the second test case at resolution 100 × 100, using the FM-LBR. Right: In hard test cases, combining strong anisotropy, metric discontinuity, and low grid resolution, the extracted path exhibits zig-zag artifacts (detail of the third test case at resolution n = 200). estimate on x i+1x i , and rescaling by a factor h, we obtain as announced [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF]. Case k = 1: choosing λ = 1, j = 1, and observing that v = v 1 , we obtain u ′ ≤ u + 1 × vv 1 = u . The second case begins with a lemma.

i x i 1 u i u i 1 Γ v x i v 1 n v 2 n x 0 x 1 x 2 x 3
Lemma A.2. Let v 1 , v 2 ∈ IR 2 , and let v := α 1 v 1 + α 2 v 2 , with α 1 , α 2 > 0. Let µ := √ 2, and let w 1 := v 1 µ -v 2 /µ, w 2 := v 2 µ -v 1 /µ.
Then one can choose λ ∈ IR + and 1 ≤ j ≤ 2, such that v jλv is positively proportional to any of the following vectors: w 1 , w 2 , -w 1 or -w 2 , with a proportionality constant 0 < ν ≤ µ.

Proof. In the case of w 1 , choose j := 1, λ := 1/(α 1 +µ 2 α 2 ), so that ν = 1/(µ+µ -1 α 1 /α 2 ) ≤ 1/µ. In the case of -w 1 , choose j := 2, λ := 1/(α 2 + µ -2 α 1 ), so that ν = 1/(µ -1 + µα 2 /α 1 ) ≤ µ. The cases of w 2 and -w 2 are similar.

Case k = 2. We use the notations of the lemma. Denoting by A the matrix of lines w 1 , w 2 , there exists 1

≤ k ≤ 2 and ε ∈ {-1, 1} such that √ 2 u, εw k ≥ u, w 1 2 + u, w 2 2 = Au ≥ A -1 -1 u . (22)
The norm A -1 -1 is estimated as follows.

| det A| = (µ 2 -µ -2 )| det(v 1 , v 2 )| ≥ µ 2 -µ -2 , since det(v 1 , v 2 ) is a non-zero integer. On the other hand w j ≤ (µ + µ -1 ) max{ v 1 , v 2 } ≤ C 0 R with C 0 = 2(µ + µ -1 ), hence A ≤ w 1 2 + w 2 2 ≤ √ 2C 0 R. Finally A -1 -1 = | det A|/ A ≥ C 1 /R where C 1 := (µ -µ -1 )/( √ 2C
0 ). We next choose λ and j, using the previous lemma, so that v jλv = νεw k , with ε, k as in [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF], and where 0 < ν ≤ µ := √ 2. Hence, using (21): 3 } as announced, which concludes the proof.

u ′ 2 ≤ u -(v j -λv) = u -νεw k 2 = u 2 -2ν u, εw k + ν 2 w k 2 ≤ u 2 -2ν A -1 -1 u + ν 2 w k 2 ≤ u 2 -2C 1 u /R + C 2 0 R 2 . If u ≥ C 2 R 3 , with C 2 := C 2 0 /(2C 1 ), then u ′ ≤ u . If u is below this bound, then choosing λ = 0 in (21) yields u ′ ≤ u + v 1 ≤ C 2 R 3 + 2R. Thus u ′ ≤ max{ u , CR

B Stencil size, discretization errors and metric regularity

Experience in the discretization of Partial Differential Equations (PDEs) tells that robust and accurate numerical schemes are usually based on small, localized and isotropic stencils. The FM-LBR, which achieves causality in the eikonal equation by the use of long range, sparse and highly anisotropic stencils, seems to violate this principle. We give in this section a heuristic analysis of its accuracy, which explains its excellent performance in the third and fourth tests (once a source of puzzlement for the author and the reviewers), but also exposes some weaknesses.

Let us emphasize that the computational domain Ω is equipped with two geometries.

• The extrinsic Euclidean geometry, inherited from the embedding Ω ⊂ IR m .

• The intrinsic Riemannian geometry, given by the Riemannian metric M on Ω, which is part of eikonal PDE structure.

The AGSI stencils are small, localized, with respect to the extrinsic Euclidean distance. The FM-LBR stencil at a point z ∈ Ω is built from an M(z)-reduced basis, see §2.2, which consists of the smallest linearly independent vectors of ZZ m in the local norm • M(z) . Hence this stencil is by construction small and localized in the twisted perspective of the intrinsic Riemannian

The interpolation error ( 25) is harder to estimate, yet in favor of the FM-LBR stencil let us mention the intrinsic 1-Lipschitz regularity of the solution: for all x, y ∈ Ω |D(x) -D(y)| ≤ D(x, y).

In the special case of a constant metric, where [START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF] is exact and all error comes from (25), the FM-LBR stencil offers the best accuracy, see [START_REF] Mirebeau | On the Accuracy of Anisotropic Fast Marching[END_REF].

we observe that this map is unimodular: ∇ d(z) M(z) -1 = 1. The value d(z) can also be regarded as the geodesic distance from z to a point at infinity in the direction of -Du. The characteristic curves of the solution are parallel straight lines, of direction Du. Although the present discussion is about accuracy rather than CPU time, let us mention that thanks to this special property, and as pointed out by A. Vladimirsky, the AGSI converges in a single pass in the numerical tests below (provided its priority queue is suitably initialized: the upwind boundary points z must be sorted by increasing values of their scalar product z, Du with the characteristics direction Du). As a result, and in contrast with §3, the AGSI CPU time is here only half of the FM-LBR CPU time. Our fifth and last numerical test, involves a metric depending on three parameters κ ≥ 1, θ, ϕ ∈ IR: denoting e θ := (cos θ, sin θ), and (x, y) ⊥ := (y, -x), M(z; κ, θ, ϕ) := M (κ, θ) exp( z, e ϕ ), with M (κ, θ) := κe θ e T θ + κ -1 e ⊥ θ (e ⊥ θ ) T .

For each matrix M = M (θ, ϕ), given by its condition number κ, and anisotropy direction θ, we consider different unit vectors u = e ϕ , given by their angle ϕ with respect to the x-axis. For ϕ = θ the FM-LBR is favored, since K 0 /K 1 takes its maximal value √ κ. For ϕ = -π/4, the AGSI is favored, since K 0 /K 1 is close to its minimal value 1/ √ κ, and the interpolation error [START_REF] Tsai | Fast sweeping algorithms for a class of Hamilton-Jacobi equations[END_REF] vanishes. The domain is the [-2, 2] 2 square discretized on a 100 × 100 grid. The FM-LBR stencil is included in a square of (2w + 1) × (2w + 1) pixels, where w respectively equals 2, 3, and 5 for the three different pairs (κ, θ) in Table 1. The boundary condition ( 29) is applied on a the upwind part of the square boundary, on a layer of width w. The numerical tests presented in Table 1 are typical of the authors experience, and tend to agree with the above heuristical error analysis. Note however that the FM-LBR performance is unexpectedly bad The accuracy advantage of the AGSI is larger than a factor 6 for the last set of parameters in Table 1, and it may of course grow unboundedly as the anisotropy ratio κ tends to infinity, for suitable angles θ, ϕ. Anisotropy does therefore, sometimes, play against the FM-LBR accuracy.

Figure 1 :

 1 Figure 1: Some classical stencils used in the discretization of two dimensional (left) or three dimensional (right) eikonal equations. The meshes are M -acute, a property implying discrete causality, see Definition 1.2, for matrices M which are diagonal or of anisotropy ratio κ(M ) bounded by respectively 1, 1 + √ 2, 1, ( √ 3 + 1)/2 (from left to right).

Figure 2 :

 2 Figure 2: Connectivity of the 2D (left) and 3D (center left, interior edges omitted, missing boundary edges are symmetric w.r.t. the origin) FM-LBR meshes T (M ), see Proposition 1.6 (note that -e i = j =i e j ). Mesh T (M ) (center right) associated to M ∈ S + 3 of eigenvalues 3 2 , 1, 1, and eigenvector (3, 1, 2) for the first eigenvalue. Right: unit ball for the norm • M .

Proposition 1 . 8 .

 18 Let Ω ⊂ IR m be open and bounded, and let M

  satisfies the (Stability) property. Proposition 1.9. Let Ω ⊂ IR m be an open bounded set, equipped with a Riemannian metric M ∈ C 0 (Ω, S + m ), and an admissible family (T (z)) z∈Ω of meshes. For all h > 0 let Z h := hZZ 2 , and for all z ∈ Z h ∩ Ω consider the stencil V h (z) := z + hT (z). Then the solutions d h : Z h → IR ∪ {+∞} of the discrete system (3) converge uniformly as h → 0 to the viscosity solution D of the eikonal equation (1):

Figure 4 :

 4 Figure4: The unit sphere {u ∈ IR 2 ; u M = 1}, an M -reduced basis (u, v), and the boundary of the FM-LBR mesh T (M ), for some M ∈ S + 2 of anisotropy ratio κ(M ) ranging from 1 to 15, and eigenvector (cos(3π/8), sin(3π/8)) associated to the small eigenvalue.

Theorem 2 . 2 (

 22 Nguyen, Stelhé, 2009). There exists an algorithm which, given a matrix M ∈ S + m as input, 1 ≤ m ≤ 4, produces an M -reduced basis of ZZ m and has the numerical cost O(1 + ln κ(M )).

1 2 since b 1 ≥ 1 , 1 2

 111 and recalling that b i M ≤ M , we obtain the second announced estimate.Proof of[START_REF] Lagrange | Recherches d'arithmétique[END_REF]. Remark that bi + z / ∈ b 1 ZZ + • • • + b i-1 ZZ,since otherwise the basis element b i would be a linear combination of b 1 , • • • , b i-1 , b i+1 , • • • , b m . Definition 2.1 thus implies:

  form a basis of ZZ m , these coefficients are integers.Up to reordering the vertices v 1 , • • • , v m , we may assume that α 1 , • • • , α k are positive, and that α k+1 , • • • , α m are zero, for some 1 ≤ k ≤ m. Let l be the smallest integer such that z

Corollary 2 . 5 . 3 implies b l 2 M 2 Mb l 2 M ≤ b l 2 M + b 1 2 M

 25322222 Let M, N ∈ S + m , with m ≤ 4, be such that dist(M, N ) < ln(1 + κ(M ) -2 )/4. (14) Let (b 1 , • • • , b m ) be an M -reduced basis of ZZ m , and let T be an N -reduced mesh. Then b 1 , • • • , b m and -b 1 , • • • , -b m are vertices of T .Proof.Let 1 ≤ l ≤ m. We have b l / ∈ b 1 ZZ + • • • + b l-1 ZZ, since (b i ) m i=1 is a basis, and Proposition 2.e 4 dist(M,N ) < b l + κ(M ) -2 .Hence b l is a vertex of T , and likewise -b l , by Lemma 2.4.We finally estimate the radius of the FM-LBR meshes T (M ), see Proposition 1.6, in terms of the condition number of the matrix M ∈ S + m . This concludes the proof of Proposition 1.8. Corollary 2.6. Let M ∈ S + m , with m ∈ {2, 3}. Then any vertex e of T (M ) satisfies e ≤ C m κ(M ), with C 2 := 2 and C 3 := 4. Proof. Given m linearly independent vertices (b i ) m i=1 of T (M ), one can express any other vertex under the form e = m i=1 α i b i . A simple check by exhaustive enumeration shows that m i=1 |α i | ≤ C m (these coefficients are independent of M ), hence e ≤ C m max{ b i ; 1 ≤ i ≤ m}. Applying Corollary 2.5 with M = N , we find that the vertices of T (M ) contain an M -reduced basis (b i ) m i=1 , which by Proposition 2.3 satisfies b i ≤ κ(M ) for all 1 ≤ i ≤ m. This concludes the proof.

Figure 5 :

 5 Figure 5: Level lines of the solutions of the two dimensional test cases.

Figure 7 :

 7 Figure 7: Reference solution for the third test case (left). The Riemannian metric M is anisotropic only on a thin band along a spiraling curve, wide of a few grid points. Detail at resolutions n × n, where n equals 200 (center left), 500 (center right) and 1000 (right).

8 Figure 8 :

 88 Figure 8: CPU Time (left, in seconds), L ∞ error (center), and averaged L 1 error (right) of the FM-LBR, FM-8 and AGSI, at several resolutions ranging from 120 to 1200 (log-log scale).

Figure 9 :

 9 Figure9: Visual comparison of the accuracy of three algorithms, at three resolutions, in the 2D test case. Qualitatively, the approximate geodesic has the right behavior for a resolution as low as 170 × 170 with the FM-LBR, and 1000 × 1000 with the AGSI. This is presumably never the case for the FM-8, which is not consistent here.

Figure 10 :

 10 Figure10: Results of the FM-LBR in the fourth, 3D, test case. Iso-surface {d(z) = 2} (left), and shortest path joining the points (0, 0, 0) and (3, 0, 0) (center). Detail of the discrete points (represented by small cubes), in the neighborhood of the curve Γ(t) = (cos ω 0 t, sin ω 0 t, t), for which the Riemannian metric is not euclidean (right).

  x

Figure 11 :

 11 Figure 11: Left: Notations for the minimal path computation; the contour of the stencil V (x i ) is shown dotted. Center left: Grid points x 0 , • • • , x r ∈ Ω∩Z, corrections u 0 , • • • , u r ∈ IR m shown as arrows, and piecewise linear path γ. Center right: Grid Points (x i ) ri=1 and piecewise linear path γ in the second test case at resolution 100 × 100, using the FM-LBR. Right: In hard test cases, combining strong anisotropy, metric discontinuity, and low grid resolution, the extracted path exhibits zig-zag artifacts (detail of the third test case at resolution n = 200).

  The following illustrative example was proposed by A. Vladimirsky: let M ∈ S + m , let u ∈ IR m , and let M : IR m → S + m be the Riemannian metric defined byM(z) := M exp( u, z ). (28)Assume for normalization that det(M ) = 1, so that M = M identically, and D(x, y) = x-y M for all x, y ∈ IR m . Then, withD := M -1 : dist(M(x), M(y)) = | u, xy | ≤ u xy , dist(M(x), M(y)) = | u, xy | ≤ u D xy M = u D D(x, y).The Lipschitz regularity constants are therefore K 0 = u , andK 1 = u D . The discretization(24) is hence likely more accuracte with the AGSI stencil if u ≤ u D , and with the FM-LBR stencil otherwise. Defining for all z ∈ IR m d(z) := exp( u, z )/ u D ,

  Algorithm 1 The Fast Marching algorithm, with static stencils adapted to a metric. Input: The values M(z) of a Riemannian metric, for all z ∈ Z ∩ Ω. Construct causal stencils V (z), with respect to M(z), at all points z ∈ Z ∩ Ω. Construct the reversed stencils, defined by V [y] := {z ∈ Z ∩ Ω; y is a vertex of V (z)}. Initialize a (mutable) table d : Z → IR, to +∞ on Z ∩ Ω, and to 0 elsewhere. Initialize a (mutable) boolean table b : Z → {trial, accepted} with b(y) ← trial iff V [y] = ∅. While there remains a trial point (i.e. y ∈ Z such that b(y) = trial) do Denote by y a trial point which minimizes d, and set b(y) ← accepted.

	For all x ∈ V [y], set d(x) ← min{d(x), Λ(d, x; b, y)}. Output: The map d : Z → IR.

  (this interval is enough, since the dependence in θ is π/2-periodic and even). In average over theta, CPU times are 0.21s, 0.20s, 1.37s, 1.31s, L ∞ errors 5.16, 7.64, 6.86, 8.57 and averaged L 1 errors 1.34, 2.58, 1.95, 2.40 for the FM-LBR, FM-8, AGSI and MAOUM respectively. All errors are multiplied by 100, for better readability.

					L error
	CPU time L ∞ error L 1 error	FM-LBR FM-8 AGSI MAOUM First test 0.19 0.19 1.01 1.28 3.99 1.47 1.62 8.80 1.13 0.53 0.51 2.33	14 12 2 4 6 8 10	FM 8 AGSI MAOUM FM LBR
	CPU time L ∞ error L 1 error	First test, rotated by π/6 0.20 0.21 1.44 1.31 5.52 12.5 9.45 8.56 1.46 3.42 2.51 2.52	0.2 FM 8 L 1 error 5 4	0.4	0.6	0.8	Θ
	CPU time L ∞ error L 1 error	0.076 2.90 1.03	Second test 0.079 0.77 3.03 3.67 1.30 1.40	0.36 7.66 2.3	3 2 1		AGSI FM LBR	Θ
						0.2	0.4	0.6	0.8
							, the FM-LBR takes

MAOUM

Figure 6: Tables of CPU time in seconds, L ∞ error an averaged L 1 error (left). Accuracy, in the first test rotated by an angle θ ∈ [0, π/4]

Table 1 :

 1 9 in case * , and unexpectdly good in cases †. Metric z → M(z; κ, θ, ϕ), see[START_REF] Zhao | A Fast Sweeping Method for Eikonal Equations[END_REF], on the [-2, 2] 2 square discretized on a 100 × 100 grid. The most accurate algorithm, among the AGSI and the FM-LBR, can in most cases be predicted by the Lipschitz constants ratio K 0 /K 1 , at least when it is far from 1. The star points out an exception. CPU time approximatively 10 ms for the AGSI, and 20 ms for the FM-LBR. Numerical errors multiplied by 100 for better readability.

	κ	θ	ϕ	K 0 /K 1 = u / u D FM-LBR AGSI FM-LBR AGSI L 1 error L ∞ error
			π/3	1.73	2.78	6.96	40.2	60.5
	3	π/3	π/6	1	2.80	3.07	17.1	18.5
			-π/4 3π/8	0.59 3.16	3.95 3.74	1.45 9.45	37.4 38.4	13.7 79.5
	10 3π/8	1.48	1	1.34	2.03	7.44	11.4
			-π/4 π/8	0.34 5.48	2.92 3.89	0.83 6.62	27.3 94.3 *	7.82 61.0
	30 π/8	0.57	1	0.91 †	2.55	4.62 †	12.9
			-π/4	0.20	3.24	0.49	29.4	4.51

This implementation is used in our numerical experiments. See Corollary

2.6 for a proof with no assumption.

Other interpolation schemes could be used, such as piecewise linear interpolation on a trivial periodic mesh, provided one can control the Lipschitz regularity constant and the support of the interpolated function.

As suggested in[START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], the stopping criterion tolerance for the AGSI iterations is set to 10 -8 .

The grid is triangulated with a trivial periodic mesh, for the AGSI and the MAOUM. As a result the AGSI uses a 6 point stencil.

The FM-LBR stencil V (z), at each grid point z ∈ Ω ∩ Z, is built[START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF] from the M(z)-reduced mesh T (M(z)) of Proposition 1.6, except if the matrix M(z) is detected to be exactly diagonal. In that case we use the standard 4 vertices neighborhood in 2D (resp.

vertices in 3D), which is an M(z)-reduced mesh, see Figure1(left and center right). This modification has little impact on accuracy or CPU time, but avoids to pointlessly break symmetry.[START_REF] Conway | Low-dimensional lattices. VI. Voronoi reduction of threedimensional lattices[END_REF] Note that the FM-LBR memory requirement is higher than that of the AGSI and FM-8, see Remark 1.10.

The stencils for the MAOUM are here built using the ComputeUpdateSet stencil construction routine described and used in all the numerical experiments of[START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF]. The paper[START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF] also outlines sufficient conditions (called δ-NGA or DRB) for anisotropic stencils to be causal, but no explicit anisotropic stencil construction.

[START_REF] Gonzales | On Deterministic Control Problems: an Approximate Procedure for the Optimal Cost, I, the Stationary Problem[END_REF] All timings obtained on a 2.4Ghz Core 2 Duo, using a single core. Timings of the FM-LBR include the stencil construction, which typically accounts for 25%.

). The L ∞ and L 1 error curves suggest that the FM-8 is not consistent in this test case, contrary to the FM-LBR and the AGSI. At the resolution 1000 × 1000, typical in image analysis, the FM-LBR cuts the L ∞ error by 80% and the L 1 error by 75% with respect to the AGSI, while reducing CPU time from 11

The FM-LBR in that case produces large errors close to the downwind boundary, presumably due to its large stencils. Removing a 5 pixel band on the boundary yields in this case the L ∞ errors: 32.3 (FM-LBR) and 45.8 (AGSI), in favor of the FM-LBR as predicted from the Lipschitz constants ratio.

Acknowledgement. The author thanks Pr A. Vladimirsky for constructive discussions on the choice of test cases and the FM-LBR accuracy, and Pr P. Q. Nguyen for pointing out the notion of obtuse superbase of a lattice.

We propose a method for the computation of minimal paths, which trades high order accuracy for robustness, and never gets stuck if the eikonal solver is Dijkstra inspired. It takes advantage of the specific form (3) of the discretization of the eikonal equation, and does not rely on black box routines. It is also parameter free: there is not interpolation order or gradient step to adjust.

This work was partly supported by ANR grant NS-LBR ANR-13-JS01-0003-01.

A Robust extraction of minimal paths

Obtaining the shortest path joining two given points is essential in motion planning control problems [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF], as well as in the envisioned application to tubular structure centerline extraction [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF]. This involves solving the Ordinary Differential Equation (ODE) [START_REF] Mirebeau | Efficient Fast Marching with Finsler Metrics[END_REF], a task less trivial than it seems. The author is conscious that a comparison of minimal paths, as on Figure 9, reflects the properties of the ODE solver (and the time spent adjusting its sometimes numerous parameters), as much as those of the eikonal solver. This is done nevertheless due to the importance of minimal paths in applications. Eikonal solvers based on the discrete fixed point problem [START_REF] Barth | Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains[END_REF], such as the FM-LBR, FM-8 and AGSI, provide at each grid point x ∈ Ω ∩ Z an estimate d(x) of the distance D(x), and in addition an estimate v(x) of the direction and orientation of the distorted negative gradient -M(x) -1 ∇ D(x), of the form:

where the integer 1 ≤ k ≤ m, the positive coefficients (α j ) k j=1 and the vertices (z j ) k j=1 of ∂V (z) are the barycentric coordinates of the point y ∈ ∂V (x) achieving the minimum in the Hopf-Lax update operator (4), in the face of ∂V (x) containing y and of minimal dimension.

From this point, a typical approach to solve ( 16) is to extend the values of d or v to the continuous domain Ω via an interpolation procedure, and then to use a black box ODE solver or a Runge Kutta method. Note that the accuracy usually expected from these high order methods is mostly doomed, since the discretization (3) of the eikonal equation is only first order, and since the vector field M -1 ∇ D is discontinuous both at "caustics" and discontinuities of M. A more significant issue is that computations frequently get stuck, despite the use of state of the art and/or commercial interpolation methods and ODE solvers, see e.g. [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF] Algorithm 2 Minimal path computation, starting from a given grid point x 0 ∈ Ω ∩ Z. Initialisation: x ← x 0 , u ← 0 (mutable variables). While x is not an initial source point, do

Denote by z 1 , • • • , z k the grid points appearing in the expression (17

The successive iterations of Algorithm 2 generate grid points x 0 , • • • , x r ∈ Z, and small correcting offsets u 0 , • • • , u r ∈ IR m . If the causality property holds, see Proposition 1.1, then the values (d(x i )) r i=0 are strictly decreasing, so that the algorithm is guaranteed to terminate. The piecewise linear path γ : [0, r] → Ω, parametrized so that γ(i)

for non-integer t ∈ [0, r], where the constants (λ i ) r-1 i=0 are the minimizers in the second step of the while loop, and the vector v(z), z ∈ Ω ∩ Z, is defined in [START_REF] Mirebeau | Adaptive and anisotropic finite element approximation : Theory and algorithms[END_REF]. The particularity of our path extraction method is that the direction field v is not evaluated on the curve γ, but at the points (x i ) r i=0 which remain close as shown by the next proposition (the involved exponent 3 seems to be either over-estimated or a rare worst case scenario, in view of the experiments Figure 11). 2 ), and an admissible family (T (z)) z∈Ω of meshes with "Boundedness" constant R. (See Definition 1.7, R κ(M) for the FM-LBR). Consider a cartesian grid Z of scale h, equipped with the corresponding stencils, see (6) and [START_REF] Fehrenbach | Sparse Non-Negative Stencils for Anisotropic Diffusion[END_REF], and solve the discrete system (3). Let γ ∈ C 0 ([0, r], Ω) be a path extracted with Algorithm 2 and parametrized as in [START_REF] Nguyen | Low-dimensional lattice basis reduction revisited[END_REF]. Then for all t ∈ [0, r] γ(t) -

where C is an absolute constant (i.e. independent of Ω, (T (z)) z∈Ω , M, or the path origin).

The rest of this appendix is devoted to the proof. Consider a fixed 0 ≤ i < r, and observe that for all t ∈ [i, i + 1[ one has, since γ is linear on this interval and since γ(i) = x i + u i :

In order to avoid notational clutter, we denote x := x i , x ′ := x i+1 , u := u i , u ′ := u i+1 and v := v(x i ). Let k ∈ {1, 2} and let z j , 1 ≤ j ≤ k, be the vertices of ∂V (x) appearing in the expression (17) of the discrete negative gradient v.

We assume without loss of generality that the grid is Z := ZZ 2 , so that v j := z j -x is a vertex of ∂T (x), for all 1 ≤ j ≤ k. Hence v j ≤ R by (Boundedness) and if k = 2 then det(v 1 , v 2 ) is a non-zero integer. In particular x ′x = v j ≤ R, for some 1 ≤ j ≤ k. By construction (second step of the while loop in the path computation), we have for any λ ∈ IR + and any

We prove in the following an upper bound of the form u ′ ≤ max{ u , CR 3 }, which by an immediate induction argument implies u i ≤ CR 3 for any 1 ≤ i ≤ r. Using [START_REF] Peyré | Geodesic Methods in Computer Vision and Graphics[END_REF], our previous distance. We refer to [START_REF] Mirebeau | On the Accuracy of Anisotropic Fast Marching[END_REF] for a quantitative estimate of the size of the FM-LBR stencil, from this perspective, in average over all orientations of the discretization grid. To better reflect the shapes of these stencils, which (as much as possible) adapt their orientation and aspect ratio to the metric M, but keep a constant volume (they are built of (m + 1)! simplices of volume h m /m!), we introduce a normalized metric M: for all z ∈ Ω

so that det( M(z)) = 1 identically. For all x, y ∈ Ω we denote by D(x, y) (resp. D(x, y)) the Riemannian distance [START_REF] Alton | An Ordered Upwind Method with Precomputed Stencil and Monotone Node Acceptance for Solving Static Hamilton-Jacobi Equations[END_REF] on Ω associated to M (resp. M).

Bellman's optimality principle applied to the solution D of the eikonal equation ( 1), reads:

The Hopf-Lax update operator (4) reflects this identity at discrete level, up to two approximations:

where V (x) denotes the stencil at x (which is given under the form of a triangulation of a neighborhood of x), I V the linear interpolation operator on V , and y ∈ ∂V . If one ignores the issue of the scheme causality, then the choice of stencil should be dictated by the local regularity properties of the quantities that are approximated on it. Discretization [START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF] amounts to approximating the metric M with the constant M(x) on the stencil V (x). Such piecewise constant approximation errors are controlled by Lipschitz regularity constants. The natural distance on S + m is defined in [START_REF] Li | Numerical schemes for Hamilton-Jacobi equations on unstructured meshes[END_REF], but the distance on Ω should be chosen appropriately so as to reflect the geometry of the computational stencils. Indeed, the AGSI (resp. FM-LBR) stencil at a point x ∈ Ω is heuristically not much different from ball centered at x and of radius the grid scale h, defined with respect to the euclidian distance (resp. Riemannian distance D). The AGSI stencil should therefore be preferred if the metric M has a small Lipschitz regularity constant K 0 with respect to the extrinsic Euclidean distance:

On the other hand, the FM-LBR stencil is more suitable for metrics which have a small Lipschitz regularity constant K 1 with respect to the intrinsic (up to the normalization ( 23)) distance D:

In the fifth numerical example below, the most accurate of these two methods can indeed be guessed from the ratio K 0 /K 1 . Regularity conditions of the form [START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF] arise naturally in the study of anisotropic mesh generation, see Part III of [START_REF] Mirebeau | Adaptive and anisotropic finite element approximation : Theory and algorithms[END_REF]. The Riemannian metric involved in the third and fourth numerical tests of this section, inspired by applications to tubular structure segmentation [START_REF] Benmansour | Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement[END_REF], varies slowly in the direction of the eigenvector associated to the small eigenvalue of M(z) (the direction of the tube), but quickly (in fact discontinuously) in the orthogonal direction. Thus, although discontinuous, it is heuristically not far from satisfying [START_REF] Vladimirsky | Label-setting methods for Multimode Stochastic Shortest Path problems on graphs[END_REF], which explains the exceptional performance of the FM-LBR on these specific examples.