Jean-Marie Mirebeau

Fast Marching methods for Curvature Penalized Shortest Paths

We introduce numerical schemes for computing distances and shortest paths with respect to several planar paths models, featuring curvature penalization and data-driven velocity: the Dubins car, the Euler/Mumford elastica, and a two variants of the Reeds-Shepp car. For that purpose, we design monotone and causal discretizations of the associated Hamilton-Jacobi-Bellman PDEs, posed on the three dimensional domain R 2 × S 1 . Our discretizations involve sparse, adaptive and anisotropic stencils on a cartesian grid, built using techniques from lattice geometry. A convergence proof is provided, in the setting of discontinuous viscosity solutions. The discretized problems are solvable in a single pass using a variant of the Fast-Marching algorithm. Numerical experiments illustrate the applications of our schemes in motion planning and image segmentation.

Introduction

In this paper, we develop numerical schemes for computing distance maps and globally minimal paths with respect to data driven costs depending on the local path position, orientation, and curvature. We address a variety of models including two variants of the Reeds-Shepp car [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], the Euler-Mumford elastica [START_REF] Euler | Methodus inveniendi Lineas curvas maximi minive proprietate gaudentes[END_REF][START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF], and the Dubins 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]. Their qualitative features differ widely: depending on the model, minimal paths may or may not be smooth, and the associated distance may or may not be continuous. For that purpose, we discretize generalized eikonal equations, also called first order static Hamilton-Jacobi-Bellman PDEs, with a unified approach relying on a tool from lattice geometry named Voronoi's first reduction of quadratic forms, also considered in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF]. Our discretizations are monotone and causal, hence can be solved using the fast marching algorithm/dynamic programming principle, with complexity O(N ln N ) where N is the number of points in the discrete domain. Numerical experiments, presented in §5 and in [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF], illustrate the potential applications of our methods in motion planning control problems and medical image segmentation tasks.

In the models of interest to us, the cost of a smooth path x : [0, T ] → R 2 , parametrized at unit speed, takes the form T 0 α(x(s), ẋ(s)) C(|ẍ(s)|) ds.

(1)

One of the considered models, the Reeds-Shepp model with reverse gear, also applies to piecewise smooth paths with cusps, whose total cost is the sum of the costs of the smooth sections.

Our objective is to numerically compute globally minimal paths for this energy, constrained to a subdomain, with prescribed endpoints and tangents at these endpoints, and regarding the terminal time T as a free parameter. We denoted by α : R 2 × S1 →]0, ∞[ a continuous cost function, depending on the path position and orientation, and usually data driven in applications. The path curvature, |ẍ(s)| in (1), denoted by κ in (2) below, is penalized using a second cost function C : R → R + , for which we consider the following three instantiations:

C RS (κ) := 1 + |ξκ| 2 , C EM (κ) := 1 + |ξκ| 2 , C D (κ) := 1 if |ξκ| ≤ 1, +∞ otherwise.
(2)

The costs C RS , C EM and C D correspond respectively to the Reeds-Shepp car, the Euler-Mumford elastica, and the Dubins car models. The parameter ξ > 0 has the dimension of a length and should be interpreted as a typical radius of curvature. Our approach easily extends to models where ξ : R 2 × S 1 →]0, ∞[ is a continuous and positive function of the path position and orientation, similarly to the cost α. The Reeds-Shepp car (with or without reverse gear), the Euler-Mumford elastica and the Dubins car are classical path models involving increasingly strong penalizations of curvature. Their qualitative properties are strikingly distinct, as illustrated on Figure 1 and discussed below.

• The Reeds-Shepp cost C RS (κ) = 1 + |ξκ| 2 is used to model slow vehicles, typically wheelchairs. Following [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF] we consider two models based on this cost, referred to as the Reeds-Shepp reversible (RS±, as originally considered in [START_REF] Reeds | Optimal paths for a car that goes both forwards and backwards[END_REF]) and Reeds-Shepp forward (RS+) models, and where the vehicle respectively may, or may not, shift into reverse gear. Minimal paths for the reversible and forward models distinguish themselves by the presence of cusps and of in place rotations of the path orientation, respectively, see Figure 1. The latter happen at the path endpoints and sometimes at the corners of obstacles, and are admissible since the curvature cost C RS (κ) only grows linearly asymptotically as κ → ∞. 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 and the description of a semi-lagrangian PDE discretization of the Reeds-Shepp models, which is different from the one considered in this paper 1 .

The special case α ≡ 1 of a constant cost has been extensively studied by analytic means. In particular, an optimal synthesis of minimal paths is presented in [START_REF] Sachkov | Cut locus and optimal synthesis in the sub-Riemannian problem on the group of motions of a plane[END_REF] for the Reeds-Shepp reversible model. In addition, [START_REF] Duits | Association Fields via Cuspless Sub-Riemannian Geodesics in SE(2)[END_REF][START_REF] Boscain | Curve cuspless reconstruction via sub-Riemannian geometry[END_REF] study in detail the endpoint configurations for which the minimal path for the Reeds-Shepp reversible model is cuspless (such path are smooth, and are also minimal for the Reeds-Shepp forward model).

• The Euler-Mumford cost C EM = 1 + |ξκ| 2 has the physical interpretation of the bending energy of an elastic bar [START_REF] Euler | Methodus inveniendi Lineas curvas maximi minive proprietate gaudentes[END_REF], when the data driven cost is identically constant α ≡ 1 and the terminal time T is fixed. The relevance of this model for image processing and segmentation was first outlined in [START_REF] Mumford | Elastica and computer vision[END_REF]. Contrary to earlier works of the author [START_REF] Chen | Global Minimum for Curvature Penalized Minimal Path Method[END_REF][START_REF] Chen | Global Minimum for a Finsler Elastica Minimal Path Approach[END_REF], the PDE discretization introduced in this paper for this model obeys a causality property which makes the Fast-Marching algorithm applicable.

• The Dubins cost C D penalizes euclidean path length only, unless curvature exceeds the threshold ξ -1 , in which case the path is rejected, see (2). Minimal (relaxed) paths for this cost are known, when α ≡ 1, to be concatenations of straight lines and of circular segments of radius ξ. This description is used in [START_REF] Boissonnat | Shortest paths of bounded curvature in the plane[END_REF] to design exact polynomial time solvers for the minimal Dubins path problem, in the presence of smooth obstacles. In contrast, our PDE approach is approximate by nature, but it can accommodate non-constant costs α, and can easily be extended to variants of the model involving e.g. position dependent bounds on the radius of curvature or additional state variables.

In the rest of this introduction, we present the classical mathematical theory underlying the addressed problems, fixing in the process important notation for the rest of the paper, and provide additional motivation for our work. More precisely, we show in §1.1 how the second order planar path models of interest (1) can be cast as a first order path model in a three dimensional domain equipped with a singular metric. We then discuss in §1.2 the exit time optimal control problem and the corresponding static Hamilton-Jacobi-Bellman (HJB) Partial Differential Equation (PDE). We describe in §1.3 some potential applications, in the field of medical image processing, of the numerical methods developed in this paper. We finally present in §1.4 the paper outline and list of contributions.

Dimension lifting in an orientation domain

Our implementation of curvature penalization requires to lift paths into the configuration space M := R 2 ×S 1 of positions and orientations. Other strategies have been proposed in the literature, see Remark 1.1. We use the identification S 1 ∼ = R/(2πZ) of the unit circle with a periodic interval, and denote points of the configuration space by p = (x, θ) ∈ M. The tangent space to M is independent of the base point, and denoted E := R 2 ×R. Vectors are denoted ṗ = ( ẋ, θ) ∈ E, and co-vectors p = (x, θ) ∈ E * . The unit vector of orientation θ ∈ S 1 is denoted n(θ) := (cos θ, sin θ).

A (local) metric on M is a function F : M × E → [0, ∞], which is convex and 1-homogeneous in its second argument. The cost functions C = C RS , C EM , C D define by homogenization three metrics F = F RS+ , F EM , F D as follows: for all p = (x, θ) ∈ M and all ṗ = ( ẋ, θ) ∈ E with the additional convention that F p (0) = 0, and that F p ((0, θ)) = lim λ↓0 λ C(| θ|/λ). A fourth metric F RS± , corresponding to the Reeds-Shepp model with reverse gear, is defined like F RS+ up to the constraint which is replaced with the unsigned collinearity requirement ẋ = ẋ, n(θ) n(θ).

We prove in Appendix B that (3) does indeed define a convex lower semi-continuous function w.r.t. ṗ, as this is not entirely obvious from the definition. The "unit balls" in each tangent space w.r.t. to the local metric, are referred to as the control sets. Their graphical representation, see Figure 2, provides important geometric intuition and is called Tissot's indicatrix. The control sets are denoted B(p) ⊆ E, and defined for all p ∈ Ω by

B(p) := { ṗ ∈ E; F p ( ṗ) ≤ 1}. (4)
The length of a Lispchitz path γ : [0, 1] → R 2 × S 1 , w.r.t. a metric F, is defined as length F (γ) := 1 0 α(γ(t)) F γ(t) ( γ(t)) dt, where γ(t) := d dt γ(t),

and where α : M →]0, ∞[ is a continuous cost function, previously mentioned and fixed throughout this paper. A bounded domain Ω ⊆ M is also fixed throughout the paper, and to each (local) metric F is associated a quasi-distance d F defined for all p, q ∈ Ω by

d F (p, q) := inf{length F (γ); γ ∈ Lip([0, 1], Ω), γ(0) = p, γ(1) = q}. (6)
Note that one may have d F (p, q) = d F (q, p) (lack of symmetry), or d F (p, q) = ∞ (lack or global controllability), or lim q→p d F (p, q) = 0 (lack of local controllability). For the metrics F RS+ , F RS± , F EM and F D considered in this paper, the infimum (6) is attained whenever

d F (p, q) < ∞, see Appendix A of [CMC15] or Appendix A of [DMMP16].
Remark 1.1 (Alternative approaches to curvature penalized minimal paths). To the knowledge of the author, two alternative methods have been proposed to compute curvature penalized minimal paths via dynamic programming. Paths are approximated in [START_REF] Strandmark | Shortest Paths with Curvature and Torsion[END_REF] with collections of non-superposable short splines, each determined by three or four control points with integer coordinates, and the cost assigned to a path is the sum of the costs of the spline approximants. No convergence analysis is presented, and the numerical results do not address the models considered in this paper. The author remains doubtful that this method is appropriate for models whose minimal paths feature singularities such as cusps (Reeds-Shepp reversible) or in place rotations (Reeds-Shepp forward), or are subject to a hard constraint on the radius of curvature (Dubins).

Another approach [START_REF] Liao | Globally Optimal Curvature-Regularized Fast Marching For Vessel Segmentation[END_REF] consists in using the original fast-marching scheme designed for euclidean distance computations [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF], but with the following addition: each time a point is added to the propagated front, a local backtracking is performed to estimate the curvature of the geodesic reaching this point, and the front propagation cost is locally adjusted as a result. The method uses a two dimensional value map u : R 2 → R, instead of the three dimensional one u : R 2 × S 1 → R used in this paper, hence the physical projections of the paths extracted with this method never cross each other. This contradicts the observed behavior, see Figure 1, hence this approach cannot compute all curvature penalized minimal paths.

The eikonal PDE formalism

The objective of this paper is to numerically solve the following optimal control problem: find the shortest path from the domain boundary ∂Ω to any point in Ω. The value function u : Ω → [0, ∞] for this problem reads for all q ∈ Ω u(q) := inf

p∈∂Ω d F (p, q). (7)
A useful and natural generalization of this problem is to introduce a penalty σ : ∂Ω →] -∞, ∞] depending on the path origin p, but for simplicity we stick to simplest instance (7) from a theoretical standpoint.

The function u associated to the Reeds-Shepp reversible metric is continuous; indeed, this model is sub-Riemannian, hence locally controllable by Chow's theorem [START_REF] Montgomery | A Tour of Subriemannian Geometries, Their Geodesics and Applications[END_REF]. In contrast, the function u can be discontinuous along ∂Ω for non locally controllable models, such as the Reeds-Shepp forward, Euler-Elastica and Dubins models. In fact, u may even be discontinuous in the interior of Ω, in the Dubins case, as well as in the Reeds-Shepp case for some domain shapes, see Proposition 3.3. The level sets of u are illustrated on Figure 14.

Despite its potential discontinuities, the function u : Ω → R is a (discontinuous) viscosity solution to a generalized eikonal equation:

∀p ∈ Ω, H p (du(p)) = 1 2 α(p)
( ṗ) := 1 2 F p ( ṗ) 2 . ( 9
)
Another equivalent expression involves the control sets (4): for any point p ∈ Ω, and any co-vector p

∈ E * H p (p) = 1 2 sup{ p, ṗ 2 + ; ṗ ∈ B(p)}. ( 10
)
The explicit expressions of the Hamiltonians H RS± , H RS+ , H EM and H D associated to the models of interest are provided in §2, where we also provide monotone and causal discretizations of the Hamilton-Jacobi-Bellman PDE (8). Once the solution u to (8) is computed, minimal paths from ∂Ω to any given q ∈ Ω, are extracted as solutions to the following ODE (solved backwards in time). For all t ∈ [0, T ]

γ(t) = dH γ(t) (du(γ(t))) , (11)
with the terminal condition γ(T ) = q, where we denoted T := u(q). The Hamiltonian H is differentiated w.r.t. the second variable p in the ODE (11), see e.g. Appendix B of [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF].

Applications to image processing

Minimal path methods are a major tool in image processing [START_REF] Peyré | Geodesic methods in computer vision and graphics[END_REF]. One of our general objectives is to provide fast, accurate and proven numerical schemes, able to extract minimal paths w.r.t. new models, that can be used to enhance these methods. However, the present paper is primarily focused on the algorithmic aspects of minimal path computation, similarly to [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][START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], and our experiments in §5 thus only involve synthetic data. Applications to real data will be published elsewhere in collaboration with experts in the field, as they previously were [SBD + 15, CMC17, CMC16b, DMMP16]. Nevertheless, and at the request of a reviewer, we briefly discuss three medical image processing tasks for which the algorithms introduced in this paper are relevant.

The segmentation of tubular structures, in two and three dimensional (medical) image data, is a common and challenging task. Kass et al [START_REF] Kass | Snakes: Active contour models[END_REF] suggested to extract these features as paths minimizing an energy, referred to as the snake model, featuring in particular a second order term, akin to a curvature penalization. Unfortunately, numerical implementations of this method often gets stuck in local minima of the non-convex snake energy, which makes it sensitive to the image noise and the path initialization. In order to address these robustness issues, Cohen et al [START_REF] Laurent | Global minimum for active contour models: A minimal path approach[END_REF] considered a simpler energy model, locally proportional to Euclidean path length, and for which global minimization is achievable using the standard Fast Marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF]. This second method faces difficulties in two dimensional images displaying overlays of vessels, such as retinal background images, or three dimensional images featuring almost crossing tubular structures. For this reason, dimension lifting techniques have been introduced, based on computations within an extended domain featuring an abstract extra coordinate, accounting for e.g. the extracted vessel local radius [START_REF] Li | Vessels as 4-D curves: Global minimal 4-D paths to extract 3-D tubular surfaces and centrelines[END_REF] or orientation [START_REF] Pechaud | Extraction of tubular structures over an orientation domain[END_REF]. The objective of the present paper, and of the works [START_REF] Bekkers | A PDE approach to data-driven sub-Riemannian geodesics in SE (2)[END_REF] and [SBD + 15, CMC17, DMMP16] (involving the author), is to combine the strengths of these different approaches: we extract tubular structures using curvature penalized paths, globally minimizing an energy, and obtained via computation in the orientation lifted domain R d × S d-1 . See the numerical experiments in §5.3.

Tractography and white matter fiber segmentation in diffusion-weighted MRI. The extraction of brain connectivity maps also has motivated the development of a variety of minimal path models and algorithms. The individual white matter fibers are not visible in MRI scans, in contrast with the tubular structure segmentation tasks considered above, but a macroscopic fiber Orientation Distribution Function (ODF) is computed from the signal. The ODF can be summarized in a diffusion tensor, which interestingly is a meaningful quantity from the geometrical viewpoint, giving rise to a natural adjunct Riemannian metric [FDHTV + 16]. A variety of numerical approaches have been developed to extract global brain connectivity measures from the local fiber ODF. Ordinary differential equations, sometimes enhanced with probabilisitic models [SKF + 14], are commonly used for fiber tracking. Minimal path methods have also been used, sometimes with ad-hoc modifications. For instance, the classical Fast Marching algorithm [START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF][START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF] is often augmented with a local backtracking used to dynamically adjust the front propagation speed depending on local direction [START_REF] Parker | Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging[END_REF] or curvature [CSR + 05, LRr13] of the shortest paths. A drawback of these heuristical enhancements is that the obtained "distance map" is not the viscosity solution to a proper eikonal equation, which limits its reliability, see also Remark 1.1 (second point). This motivates the introduction in [JBT + 08] of a semi-Lagrangian fast marching method based on a 26-points stencil, and able to compute genuine shortest paths w.r.t. (mildly) anisotropic Riemannian metrics2 . Curvature penalized minimal paths, such as those considered in this paper, are also applicable to fiber tractography, see for instance the (synthetic) numerical examples in [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF].

Object segmentation in two dimensional images, for instance organ delineation in medical data, can be addressed by extracting the object boundary as a path minimizing an energy, similarly to tubular structure segmentation [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Laurent | Global minimum for active contour models: A minimal path approach[END_REF]. A specificity of this application is that objects contours are naturally oriented, e.g. counter-clockwise; hence they may be extracted as minimal paths w.r.t. a non-symmetric Finsler function designed to favor paths having the object at their right and the background at their left, see [START_REF] Zach | Globally Optimal Finsler Active Contours[END_REF]. (See also [START_REF] Melonakos | Finsler Active Contours[END_REF] for more discussion on Finsler active contours.) In a similar spirit, [START_REF] Chen | Finsler Geodesics Evolution Model for Region based Active Contours[END_REF] relies on the divergence theorem and an asymmetric Finsler metric, to globally minimize the Chan-Vese region segmentation energy [START_REF] Chan | Active contours without edges[END_REF]. Finally, curvature penalized shortest paths, such as those considered in this paper, are also applicable to object segmentation, see [START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF].

Contributions and outline

In section §2 we introduce new discretizations of the generalized eikonal PDEs (8) associated to the four models of interest. For that purpose, we use an original reformulation of the Euler-Mumford Hamiltonian, and a construction based on lattice geometry. Our numerical schemes are monotone and causal, enabling the use of the single pass Fast-Marching algorithm. Their stencils, carefully designed, are sparse, anisotropic, and have a reasonably small radius.

Section §3 is devoted to the convergence analysis. We prove in §3.1 that the discretized PDEs admit uniformly bounded solutions. We establish in §3.2 the convergence of the discrete solutions upper and lower continuous envelopes towards sub-and super-solutions of (8). We discuss in §3.3 the continuity properties of the value function u (7), which are related to the pointwise convergence of the discrete solutions towards u, depending on the model and on the domain geometry.

Section §4 is devoted to the proof of a key ingredient of our PDE discretization schemes, namely Proposition 2.2 below. It is based on tools from lattice geometry known as Voronoi's first reduction of quadratic forms, obtuse superbases, and Selling's algorithm.

Numerical experiments presented in §5 illustrate the potential of our approach in motion planning and image segmentation tasks. We also compare, for validation, our minimal geodesics with those obtained using a shooting method based on Hamilton's ODE of geodesics.

Discretization

In this section, we construct finite differences discretizations of the HJB PDEs (8) associated to the different models of interest in this paper. For that purpose we derive the explicit expression of the relevant Hamiltonian, and construct an approximation of a specific form. The Reeds-Shepp, Euler-Mumford, and Dubins models are respectively addressed in §2.1, §2.2 and §2.3.

Our numerical schemes obey two fundamental properties, referred to as monotony and causality, which are recalled in the next definition, see [START_REF] Oberman | Convergent Difference Schemes for Degenerate Elliptic and Parabolic Equations: Hamilton-Jacobi Equations and Free Boundary Problems[END_REF] and [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF] respectively. Monotony implies comparison principles, which are used in the convergence analysis §3. Causality allows to solve the discretized PDE in a single pass using dynamic programming, a generalization of Dijkstra's algorithm and the Fast-Marching method, which guarantees short computation times. Definition 2.1. A (finite differences) scheme on a finite set X is a continuous map F : X × R × R X → R. The scheme is said:

• Monotone, iff F is non-decreasing w.r.t. the second and (each of the) third variables.

• Causal, iff F only depends on the positive part of the third variable.

To the scheme is associated a function R X → R X still (abusively) denoted by F, and defined by

(FU )(x) := F(x, U (x), (U (x) -U (y)) y∈X ), (12)
for all x ∈ X and all U ∈ R X .

Let us fix a grid scale h > 0 of the form 2π/k for some positive integer k, and introduce the grid M h ⊆ M := R 2 × S 1 , the finite discrete domain Ω h , and the (formal) discrete boundary ∂Ω h defined by

M h := hZ 2 × (hZ/2πZ), Ω h := Ω ∩ M h , ∂Ω h := M h \ Ω h . ( 13
)
Our discretizations of the HJB PDE (8) are presented as follows

∀p ∈ Ω h , HU (p) = 1 2 α(p) 2 , ∀p ∈ ∂Ω h , U (p) = 0, (14)
For their design, we first construct an approximation of the local hamiltonian of the model of interest, H p : E * → [0, ∞[ where p ∈ Ω, under the following form

H p (p) ≈ max i∈I j∈J α ij p, ėij 2 + + k∈K β ik p, ḟik 2 , (15)
where I, J and K are arbitrary finite sets. Here and in the rest of this paper, we denote by a + := max{a, 0} and a -:= max{-a, 0} the positive an negative parts of a scalar a ∈ R. The chosen weights are always non-negative: α ij , β ik ≥ 0, and the offsets are integral: e ij , f ik ∈ Z 3 . They depend on the base point p, in our case on the angular coordinate θ ∈ S 1 of p = (x, θ) ∈ M only. The PDE (8) is discretized using upwind finite differences, forward or symmetric

HU (p) =h -2 max i∈I j∈J α ij max{0, U (p) -U (p -h ėij )} 2 + k∈K β ik max{0, U (p) -U (p -h ḟik ), U (p) -U (p + h ḟik )} 2 .
The pattern of offsets { ėij , ḟik , -ḟik ; i ∈ I, j ∈ J, k ∈ K} is referred to as the stencil of the numerical scheme. Our PDE discretizations rely on point dependent, sparse and strongly anisotropic stencils, with reasonable radius, see Figure 3. They are designed using the following result, whose proof presented in §4 relies on Voronoi's first reduction, a tool from discrete geometry characterizing the interaction of a positive quadratic form with an additive lattice [START_REF] Schürmann | Computational geometry of positive definite quadratic forms[END_REF]. Similar techniques are used for anisotropic diffusion PDEs in [START_REF] Fehrenbach | Sparse non-negative stencils for anisotropic diffusion[END_REF], for Monge-Ampere equations in [START_REF] Mirebeau | Minimal stencils for discretizations of anisotropic PDEs preserving causality or the maximum principle[END_REF], and for eikonal PDEs associated to Riemannian, sub-Riemannian and Rander metrics in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF].

More precisely, the next proposition shows how the (squared) positive part of a linear form p → p, ṅ + can be approximated using positive parts of linear forms p → p, ė + associated to few, small, integral vectors ė ∈ Z d . This result allows in our numerical schemes to approximate the directional derivative du(p), ṅ + at a point p ∈ M using finite differences of the form

ε ė( ṅ), ė ∈ Z d , such that for all p ∈ R d p, ṅ 2 + ≤ ė∈Z d ρ ε ė( ṅ) p, ė 2 + ≤ p, ṅ 2 + + ε 2 ṅ 2 p 2 .
Furthermore the support { ė ∈ Z d ; ρ ε ė( ṅ) > 0} has at most 3 elements in dimension d = 2 (resp. 6 elements in dimension d = 3), and is contained in a ball of radius C WS /ε, where C WS is an absolute constant. In addition

ė∈Z d ρ ε ė( ṅ) ė 2 = ṅ 2 (1 + (d -1)ε 2 ).
In practice, we usually choose the relaxation parameter ε = 1/10 and obtain a support which is 5 or 6 pixels wide, see Figure 3. Our numerical method thus belongs to the category of Wide-Stencil schemes (hence the subscript to the constant C WS ).

The Reeds-Shepp car models

This section is focused on the discretization of the Reeds-Shepp forward model, postponing the discussion of the reversible model to Remark 2.7. The metric F RS+ of the Reeds-Shepp forward model is obtained by inserting the curvature cost expression C RS (κ) := 1 + (ξκ) 2 in the generic expression (3). Thus for all p = (x, θ) ∈ M and all ṗ = ( ẋ, θ) ∈ E

F RS+ p ( ṗ) = ẋ 2 + |ξ θ| 2 if ẋ = ẋ n(θ), (16)
and F RS+ p ( ṗ) = +∞ otherwise. The control set B RS+ (p) is defined as the unit ball of the metric F RS+ p , for any p = (x, θ) ∈ M, see (4). Hence it is an half ellipse, or a half disk if ξ = 1, as illustrated on Figure 2.

B RS+ (p) = {( ẋ, θ) ∈ E; ẋ 2 + |ξ θ| 2 ≤ 1, ẋ = ẋ n(θ)}.
The Lagrangian of the Reeds-Shepp model is defined as the half square of the metric (16). Hence denoting by P n ( ẋ) := ẋn, ẋ n the component of a vector ẋ orthogonal to a direction n, one has

2L RS+ p ( ṗ) = ( n(θ), ẋ 2 + + ∞ n(θ), ẋ 2 -) + ∞ P n(θ) ( ẋ) 2 + |ξ θ| 2
, where, slightly abusively, we use infinite coefficients with the convention 0 × ∞ = 0. The Hamiltonian, also presented in [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF], reads for all p = (x, θ) ∈ Ω, and all p = (x, θ)

∈ E * 2H RS+ p (p) = x, n(θ) 2 + + | θ/ξ| 2 . ( 17
)
This expression follows from the piecewise quadratic and separable structure of the Lagrangian, and from two basic lemmas on the Legendre-Fenchel dual f * of a function f , which are recalled below without proof.

Lemma 2.3 (Legendre-Fenchel dual of a separable sum). Let (e 1 , • • • , e d ) be an orthogonal basis of R d , and let f 1 , • • • , f d : R →] -∞, ∞] be proper, convex and lower semi-continuous. Then

∀x ∈ R d , f (x) := 1≤i≤d f i ( e i , x ) ⇒ ∀x ∈ R d , f * (x) = 1≤i≤d f * i ( e i , x ). Lemma 2.4 (Legendre-Fenchel dual of a quadratic function). Let a, b ∈ [0, ∞]. Then ∀x ∈ R, f (x) = 1 2 (a 2 x 2 -+ b 2 x 2 + ) ⇒ ∀x ∈ R, f * (x) = 1 2 (a -2 x 2 -+ b -2 x 2 + ).
We propose the following discretization scheme H RS+ ε,h for the Hamitonian H RS+ , which depends on a relaxation parameter ε ∈]0, 1] and on the grid scale h. Proposition 2.2 is instantiated in dimension two to provide the weights ρ ε ė(n(θ)). The offsets appearing in this expression are illustrated on Figure 3. For any discrete map U :

M h → R and any p = (x, θ) ∈ Ω h 2H RS+ ε,h U (p) := h -2 ė∈Z 2 ρ ε ė(n(θ)) max{0, U (x, θ) -U (x -h ė, θ)} 2 (18) + (ξh) -2 max{0, U (x, θ) -U (x, θ -h), U (x, θ) -U (x, θ + h)} 2 .
Proposition 2.5. The discretization scheme H RS+ ε,h is monotone and causal, for any ε ∈]0, 1] and any h > 0. It is supported on 6 points at most, at distance at most C WS h/ε from p. Furthermore if U coincides with a linear function on these points, then

H RS+ p (dU (p)) ≤ H RS+ ε,h U (p) ≤ H RS+ p (dU (p)) + ε 2 2 dU (p) 2 .
Proof. The monotony and causality properties are clearly satisfied, see also Proposition 3.6 for a general class of schemes which obey these properties. The support cardinality and radius estimates follow from Proposition 2.2. If U is locally linear around p = (x, θ) ∈ M h , then denoting dU (p) = (x, θ) and observing that

U (x, θ) -U (x -h ė, θ) = x, ė and U (x, θ) - U (x, θ -h) = θh, one obtains H RS+ ε,h U (p) = ė∈Z 2 ρ ε ė(n(θ)) x, ė 2 + + ξ -2 max{0, θ, -θ} 2 ,
and the announced estimate follows by (17) and Proposition 2.2.

Proposition 2.5 easily extends to three dimensional Reeds-Shepp forward model, posed on R3 × S 2 . We refer to [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF] for more discussion on this extension and numerical experiments 3 , and we focus our attention on two dimensional path models in this paper. Nevertheless we describe our three dimensional discretization strategy in the following lines, since it is only briefly evoked in [START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF]. For that purpose, the two dimensional unit sphere is parametrized as n(θ, ϕ) := (cos θ, sin θ cos ϕ, sin θ sin ϕ), where (θ, ϕ) ∈ A 2 := [0, π] × [0, 2π] are similar to the Euler angles. The Reeds-Shepp forward metric reads in this context:

F RS+ p ( ṗ) := ẋ 2 + |ξ θ| 2 + | φξ sin θ| 2 if ẋ = ẋ n(θ, ϕ)
and F RS+ p ( ṗ) = +∞ otherwise, for all p = (x, θ, ϕ) ∈ R 3 × A 2 and all ṗ = ( ẋ, θ, φ) ∈ R 3 × R 2 . Using, as above, the separable and piecewise quadratic structure of the Lagrangian, we obtain

2 H RS+ p (p) = x, n(θ, ϕ) 2 + + | θ/ξ| 2 + | φ/(ξ sin θ)| 2 for all p = (x, θ, φ) ∈ (R 3 ) * × (R 2 ) * .
The rectangle A 2 is equipped with the adequate boundary conditions (θ, 0) ∼ (θ, 2π), (0, ϕ) ∼ (0, ψ), (π, ϕ) ∼ (π, ψ), for all θ ∈ [0, π], ϕ, ψ ∈ [0, 2π], and discretized using a cartesian grid of scale h = π/k for some positive integer k. Proposition 2.2 is instantiated in dimension three to provide the weights ρ ε ė(n(θ, ϕ)).

H RS+

ε,h U (p) := h -2 ė∈Z 3 ρ ε ė(n(θ, ϕ)) max{0, U (p) -U (x -h ė, θ, ϕ)} 2 + (ξh) -2 max{0, U (p) -U (x, θ -h, ϕ), U (p) -U (x, θ + h, ϕ)} 2 + (ξh sin θ) -2 max{0, U (p) -U (x, θ, ϕ -h), U (p) -U (x, θ, ϕ + h)} 2 .
Proposition 2.6. The discretization scheme H RS+ ε,h is monotone and causal, for any ε ∈]0, 1], h > 0. It is supported on 11 points at most, at distance at most C WS h/ε from p. Furthermore if U coincides with a linear function on these points, then

H RS+ p (dU (p)) ≤ H RS+ ε,h U (p) ≤ H RS+ p (dU (p)) + ε 2 2 dU (p) 2 .
We do not state the proofs Proposition of 2.6 and of the following remark, since they are entirely similar to that of Proposition 2.5.

Remark 2.7 (The reversible Reeds-Shepp model). The metric of the Reeds-Shepp reversible model has the same expression as (16), except for the modified constraint: ẋ = n(θ), ẋ n(θ). As a result, the Lagrangian and Hamiltonian read, with the same conventions as above

2L RS± p ( ṗ) = n(θ), ẋ 2 + ∞ P n(θ) ( ẋ) 2 + |ξ θ| 2 , 2H RS± p (p) = x, n(θ) 2 + | θ/ξ| 2 .
The discretization scheme (18) can be adapted by appropriately modifying its first line:

2H RS± ε,h U (p) := h -2 ė∈Z 2 ρ ε ė(n(θ)) max{0, U (x, θ) -U (x -h ė, θ), U (x, θ) -U (x + h ė, θ)} 2 + • • • .
This scheme supported on 9 points, and for a linear U on these points one has the identity

H RS± ε,h U (p) = H RS± p (dU (p)) + ε 2 2 dU (p) 2 .
A similar remark applies to the three dimensional Reeds-Shepp reversible model.

The Euler-Mumford elastica model

The metric F EM of the Euler-Mumford elastica model is obtained by inserting the curvature cost C EM (κ) := 1 + (ξκ) 2 in the generic expression (3). Thus for all p = (x, θ) ∈ M and all ṗ = ( ẋ, θ) ∈ E

F EM p ( ṗ) = ẋ + |ξ θ| 2 ẋ if ẋ = ẋ n(θ),
and

F EM p ( ṗ) = +∞ otherwise. Observing that ẋ + |ξ θ| 2 ẋ ≤ 1 ⇔ ẋ 2 + |ξ θ| 2 ≤ ẋ ⇔ ( ẋ -1/2) 2 + |ξ θ| 2 ≤ 1/4
we obtain that the control sets B EM are ellipses, or disks if ξ = 1. Note, however, that the origin 0 of the tangent space E is not in their center but on their boundary, see Figure 2.

B EM (p) = {( ẋ, θ) ∈ E; ( ẋ -1/2) 2 + |ξ θ| 2 ≤ 1/4, ẋ = ẋ n(θ)}, = 1 2 ((a + 1)n(θ), b/ξ) ; a, b ∈ R, a 2 + b 2 ≤ 1 .
Lemma 2.8. The Euler-Mumford elastica Hamiltonian reads, for all p = (x, θ) ∈ M and all

p = (x, θ) ∈ E * 2H EM p (p) := 1 4 x, n(θ) + x, n(θ) 2 + | θ/ξ| 2 2 . ( 19
)
Proof. The announced result follows from (10) and from the computation

2 sup ṗ∈B EM (p) p, ṗ = x, n(θ) + sup a 2 +b 2 ≤1 a x, n(θ) + b θ/ξ = x, n(θ) + x, n(θ) + | θ/ξ| 2 .
The expression (19) suggests to regard the Euler/Mumford model as a degenerate case of Rander geometry [START_REF] Randers | On an Asymmetrical Metric in the Four-Space of General Relativity[END_REF]. This approach is considered numerically in [START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF]. Unfortunately, existing numerical schemes for eikonal equations involving Rander metrics are either limited to two dimensions [START_REF] Mirebeau | Efficient fast marching with Finsler metrics[END_REF], or lack causality [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF] and thus cannot be solved using the Fast-Marching algorithm, which significantly impacts the numerical cost and the flexibility of their implementations. In this paper we advocate for a different approach, based on a second expression of the Hamiltonian H EM , in integral form, which to our knowledge is original.

Proposition 2.9. For any point p = (x, θ) ∈ M and any co-vector p ∈ E * , one has

2H EM p (p) = 3 4 π/2 -π/2 p, (n(θ) cos ϕ, ξ -1 sin ϕ) 2 + cos(ϕ) dϕ.
Proof. The first step of the proof, left as an exercise to the reader, is to show that

π/2 -π/2 (cos(ϕ -ψ)) 2 + cos ϕ dϕ = 1 3 (1 + cos ψ) 2 , ( 20
)
for any ψ ∈ R. The second step is the claim that for any ė0 , ė1 ∈ E one has

π/2 -π/2 p, cos(ϕ) ė0 + sin(ϕ) ė1 2 + cos ϕ dϕ = 1 3 p, ė0 2 + p, ė1 2 + p, ė0 2 . ( 21
)
Indeed, if p, ė0 = p, ė1 = 0 then there is nothing to prove. Otherwise, up to rescaling p, we may assume that p, ė0 2 + p, ė1 2 = 1, thus p, ė0 = cos ψ and p, ė1 = sin ψ for some ψ ∈ R. Therefore p, cos(ϕ) ė0 + sin(ϕ) ė1 = cos ϕ cos ψ + sin ϕ sin ψ = cos(ϕ -ψ) for any ϕ ∈ R, hence (21) follows from (20). Choosing e 0 := (n(θ), 0) and e 1 := (0 R 2 , ξ -1 ) we conclude the proof.

In order to discretize the Euler-Mumford Hamiltonian, we consider a second order consistent quadrature rule on the interval [-π/2, π/2] with cosine weight. Quadrature rules on the interval [-1, 1] for the uniform cost, such as the Clenshaw-Curtis or Fejer rules [START_REF] Féjer | On the infinite sequences arising in the theories of harmonic analysis, of interpolation, and of mechanical quadratures[END_REF], are for instance easily adapted to our needs thanks to the identity

1 -1 f (t) dt = π/2 -π/2 f (sin ϕ) cos ϕ dϕ.
More precisely, let K be a positive integer, and let (α

k , ϕ k ) ∈ (R + × [-π/2, π/2]) K be such that for any twice continuously differentiable f : [-π/2, π/2] → R 1≤k≤K α k f (ϕ k ) - π/2 -π/2 f (ϕ) cos(ϕ) dϕ ≤ C K 2 sup |f (t)|; t ∈ [-π/2, π/2] , (22)
where C is independent of f and K. Note that choosing f ≡ 1 on [-π/2, π/2], one obtains

1≤k≤K α k = π/2 -π/2 cos ϕ dϕ = 2. ( 23
)
We propose the following discretization of the Euler-Mumford Hamiltonian

2H EM ε,K,h U (p) := 3 4 h -2 0≤k≤K α k ė∈Z 3 ρ ε ė(n(θ) cos ϕ k , ξ -1 sin ϕ k ) max{0, U (p) -U (p -h ė)} 2 .
It is based on the three dimensional instantiation of Proposition 2.2, applied to the vectors (n(θ) cos ϕ k , ξ -1 sin ϕ k ), 0 ≤ k ≤ K with the relaxation parameter ε ∈]0, 1].

Proposition 2.10. The discretization scheme H EM ε,K,h is monotone and causal, for any ε ∈]0, 1], K ≥ 1, h > 0. It is supported on at most 6K + 1 points, at distance at most C WS h/ε from p. Furthermore if U coincides with a linear function on these points, then

H EM p (dU (p)) -CK -2 dU (p) 2 ≤ H EM ε,K,h U (p) ≤ H EM p (dU (p)) + C(ε 2 + K -2 ) dU (p) 2 , with C = C 0 max{1, ξ -2 } for some absolute constant C 0 .
Proof. Monotony, causality, and the stencil cardinality and radius estimates can be proved similarly to Proposition 2.5. Let p = (x, θ) ∈ Ω be fixed, and let p = dU (p). Let v(ϕ

) := (n(θ) cos ϕ, ξ -1 sin ϕ) for all ϕ ∈ [-π/2, π/2]. Using Proposition 2.2 and observing that v(ϕ) ≤ max{1, ξ -1 } one obtains 0 ≤ 8 3 H EM ε,h U (p) - 1≤k≤K α k p, v(ϕ k ) 2 + ≤ 1≤k≤K α k ε 2 max{1, ξ -2 } p 2 ,
where we used the non-negativity of the weights (α k ) K k=0 . Therefore using (23)

-C 1 p 2 K -2 ≤ 8 3 H EM ε,h U (p) - π/2 -π/2 p, v(ϕ) 2 + sin ϕ dϕ ≤ C 1 p 2 K -2 + C 2 ε 2 p 2 , ( 24
)
where C 2 = 2 max{1, ξ -2 }, and where we applied (22) to the function ϕ → p, v(ϕ) 2 + , whose second derivative makes sense as an L ∞ function and is bounded by C 1 p 2 . Here C 1 = C 1 (ξ) denotes the absolute constant of (22) times an upper bound for the partial second derivative ∂ 2 /∂ϕ 2 of the trigonometric polynomial q, v(ϕ) 2 = q, (n(θ) cos ϕ, ξ -1 sin ϕ) 2 = a(q, θ, ξ) cos(2ϕ) + b(q, θ, ξ) sin(2ϕ), (25) uniformly w.r.t. all θ ∈ S 1 and all q ∈ E * such that q = 1. Note that the coefficients a and b, hence also C 1 , are O(ξ -2 ). The positive part appearing in the expression p, v(ϕ) 2 + of (24) and not in (25) is not an issue, thanks to a minor technical argument presented in Lemma 2.11 below. The announced result follows from (24) and Proposition 2.9.

Lemma 2.11. Let f ∈ C 2 (R, R), and let g(x) := f (x) 2 + for all x ∈ R. Then g ∈ L ∞ loc (R), g (x) = d 2 dx 2 (f (x) 2 ) if f (x)
> 0, and g (x) = 0 for almost every x such that f (x) = 0. Proof. Clearly g is locally Lipschitz, with derivative g (x) = 2f (x)f (x) + . This expression shows that g is continuous and also locally Lipschitz, hence g is almost everywhere twice differentiable, with the announced expression. This concludes the proof.

The Dubins car model

The metric of F D of the Dubins car model is obtained by inserting the cost function

C D (κ) := 1 if |ξκ| ≤ 1, C D (κ) = +∞ otherwise, in the generic expression (3). Hence for all p = (x, θ) ∈ M and all ṗ = ( ẋ, θ) ∈ E F p ( ṗ) = ẋ if ẋ = ẋ n(θ) and |ξ θ| ≤ ẋ , (26)
and F p ( ṗ) = +∞ otherwise. The control set B D (p) is an isosceles triangle, or a half square if ξ = 1, whose apex is the origin of E.

B D (p) = {( ẋ, θ) ∈ E; ξ| θ| ≤ ẋ ≤ 1, ẋ = ẋ n(θ)} = {(an(θ), bξ -1 ); 0 ≤ |b| ≤ a ≤ 1}.
The Dubins Hamiltonian is the square of a piecewise linear function.

Lemma 2.12. For all p = (x, θ) ∈ M and all p = (x, θ) ∈ E * , one has

2H D p (p) = max{0, p, (n(θ), ξ -1 ) , p, (n(θ), -ξ -1 ) } 2 = ( x, n(θ) + ξ -1 | θ|) 2 +
Proof. The result follows from the expression (10) of the Hamiltonian, and from the observation that the linear function p, • always attains its maximum at an extreme point of the convex set B D (p), hence at one the three vertices 0, (n(θ), ξ -1 ) and (n(θ), -ξ -1 ) of this triangle.

We propose the following discretization scheme H D ε,h , with gridscale h > 0 and relaxation parameter ε ∈]0, 1]. It relies on the three dimensional instantiation of Proposition 2.2 applied to the two vectors (n(θ), ±ξ -1 ). For any U :

M h → R and any p = (x, θ) ∈ Ω h , define H D ε,h U (p) := h -2 max ė∈Z 3 ρ ε ė(n(θ), ξ -1 )(U (p) -U (p -h ė)) 2 + , ė∈Z 3 ρ ε ė(n(θ), -ξ -1 )(U (p) -U (p -h ė)) 2 + .
Proposition 2.13. The discretization scheme H D ε,h U (p) is monotone and causal. It is supported on 13 points at most, within distance C WS h/ε from p. Furthermore if U coincides with a linear function on these points then

H D p (dU (p)) ≤ H D ε,h (p) ≤ H D p (dU (p)) + ε 2 2 dU (p) 2 .
The proof, entirely similar to Proposition 2.5, is left to the reader.

Convergence analysis

This section is devoted to the convergence analysis of the discretization schemes introduced in §2, and applied to the optimal control problem (7). We prove in Theorems 3.1 and 3.2 that the resulting discrete systems of equations are solvable using the fast-marching algorithm, and that the upper and lower continuous envelopes of the discrete solutions converge to sub-and supersolutions to the HJB PDE, as the grid scale h and relaxation parameter ε tend to 0 suitably. We also establish in Proposition 3.3 some continuity properties of the value function to our optimal control problem (7), which are related to the pointwise convergence of the discrete solutions.

For that purpose, following the notations of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], we introduce a close relative û : Ω → R to the value function u defined by (7). For any q ∈ Ω, denoting by F the metric of the model û(q) := inf

p∈R 2 \Ω d F (p, q). ( 27
)
We denote by B(p, r) ⊆ M the Euclidean unit ball of center p ∈ M and radius r > 0. The discretization grid Ω, ∂Ω h , is defined in (13).

Theorem 3.1. Let Ω ⊆ M be an open and bounded domain, and let α : M →]0, ∞[ have Lipschitz regularity. Then for any h > 0 and ε ∈]0, 1] the system

∀p ∈ Ω h , H RS+ ε,h U (p) = 1 2 α(p) 2 , ∀p ∈ ∂Ω h , U (p) = 0, (28)
admits a unique solution denoted U ε,h : M h → R. This solution can be computed using the fast marching algorithm with complexity O(N h ln N h ), where U n .

N h := #(Ω h ). Let U n := U εn,
Then u ≤ u ≤ u ≤ û on Ω.
The first part of this result, on the discretized systems, is established in §3.1. The second part, on the comparison with u and û, which are respectively the smallest super-solution and the largest sub-solution to the HJB PDE (8), is addressed in §3.2. These results are of course not limited to the Reeds-Shepp forward model. Theorem 3.2. Theorem 3.1 applies to the Reeds-Shepp reversible, the Euler-Mumford elastica, and the Dubins models as well. Obviously, the discretization schemes are H RS± ε,h , H EM ε,K,h and H D ε,h , and the metrics are F RS± , F EM and F D in problem (7), respectively. In the Euler-Mumford case, the additional discretization parameter K must be taken into account as follows: the fastmarching complexity is O(KN h ln N h ), and convergence holds for

U n = U εn,Kn,hn provided ε n → 0, h n /ε n → 0 and K n → ∞ as n → ∞.
The convergence result presented in Theorem 3.1 is incomplete until one proves that u = û on a large subset of the domain Ω. Our knowledge on this topic is gathered in the following proposition, proved in §3.3. As a side product, this result establishes some continuity properties of the value function to the optimal control problem (7), which are interesting from a qualitative point of view. The interior of a set A is denoted by int(A). The most interesting aspect of Proposition 3.3 is what it does not prove. In particular, u RS+ and u EM need not be continuous at the boundary of Ω, and u D may be discontinuous in the interior of Ω, as well as u RS+ if Ω has not the shape specified in Proposition 3.3. The assumption int(Ω) = Ω, can also be formulated int T = T as in [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], where T := M \ Ω is the target set. It forbids the presence of isolated points in the target, which are inconvenient from a theoretical perspective, although they are common in applications.

In order to simplify the proofs, we introduce the dual F * of any metric F, defined for any point p ∈ M and co-vector p ∈ E * by

F * p (p) := sup ṗ =0 p, ṗ F p ( ṗ) = sup ṗ∈B(p) p, ṗ = 2H p (p). ( 29
)
where the control set B and Hamiltonian H are defined by ( 4) and ( 9), see also (10). The dual norm F * p is positively 1-homogeneous and obeys the triangular inequality:

F * p (λp) = λF * p (p) F * p (p 1 + p2 ) ≤ F * p (p 1 ) + F * p (p 2 ) (30)
for any λ ≥ 0, p ∈ M, p, p1 , p2 ∈ E * .

Existence, uniqueness and boundedness of a discrete solution

We establish the existence and uniqueness of a solution to the discretized problem (28) for the models of interest. We also prove upper bounds on this solution which are independent of the gridscale h, assuming that the relaxation parameter satisfies 0 < ε ≤ h. (This property remains valid if 0 < ε ≤ Ch, where C is an absolute constant.) Definition 3.4. Let F be a PDE discretization scheme on a finite set X, in the sense of Definition 2.1. A discrete map U ∈ R X is called a sub-(resp. strict sub-, resp. super-, resp. strict super-) solution of the scheme F iff FU ≤ 0 (resp. FU < 0, resp. FU ≥ 0, resp. FU > 0) pointwise on X. If FU = 0, then U is a solution to the scheme.

When the scheme F is obvious from context, we simply speak of sub-and super-solution. The existence, uniqueness, and computability of the solutions to PDE schemes are discussed in the next result, using the notions of monotony and causality introduced in Definition 2.1. Theorem 3.5 is not an original contribution, but gathers classical results from [Tsi95, RT92, Obe06], see also [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF] for a proof. Theorem 3.5 (Solving monotone schemes). Let F be a monotone scheme on a finite set X s.t.

(i) There exists a sub-solution U -and a super-solution U + to the scheme F.

(ii) Any super-solution to F is the pointwise limit of a sequence of strict super-solutions.

Then there exists a unique solution U ∈ R X to FU = 0, and it satisfies U -≤ U ≤ U + . If in addition the scheme is causal, then this solution can be obtained via the Dynamic-Programming algorithm, also called Dijkstra or Fast-Marching, with complexity O(M ln N ) where

N = #(X), M = #({(x, y) ∈ X × X; FU (y) depends on U (x)}). ( 31
)
The logarithmic factor ln N in the complexity comes from the cost of maintaining a priority queue. We counted as elementary, complexity-wise, a scheme dependent local update which in our case amounts to solving a univariate quadratic equation. We refer to [START_REF] Rouy | A Viscosity Solutions Approach to Shape-From-Shading[END_REF] for more details on the Fast-Marching algorithm.

General properties of our PDE schemes. In this paragraph, we prove that the PDE discretization schemes introduced in §2 are monotone, causal, admit a sub-solution, and satisfy Property (ii) of Theorem 3.5. For that purpose we express them as specializations of a generic design, described in the next proposition.

Proposition 3.6. Let X be a finite set, and let F : R X → R X be defined by

(F U (x)) 2 := max i∈I y∈X∪∂X b(i, x, y)(U (x) -U (y)) 2 + , (32)
for all U : X → R and all x ∈ X. We denoted by I and ∂X some arbitrary finite sets, and by b :

I × X × (X ∪ ∂X) → [0, ∞[ some non-negative weights. By convention, U is extended by 0 on ∂X. Then F (λU ) = λF (U ) and F (U + V ) ≤ F (U ) + F (V )
, pointwise on X, for any λ ≥ 0 and any U, V : R X → R. Let also

HU (x) := 1 2 (F U (x)) 2 FU (x) := -a(x) + HU (x), ( 33
)
where a : X →]0, ∞[ is arbitrary. Then F is a monotone and causal scheme, admits U ≡ 0 a as sub-solution, and (1 + ε)U is a strict super-solution for any super-solution U and any ε > 0.

Proof. The 1-Homogeneity of F is obvious. The triangular inequality follows from the expression

F U (x) = max i∈I b(i, x, y)(U (x) -U (y)) + y∈X∪∂X ,
and from the basic inequality (a + b)

+ ≤ a + + b + applied to a = U (x) -U (y), b = V (x) -V (y).
The scheme F is monotone since F is non-decreasing w.r.t. the differences (U (x) -U (y)) y∈X , and it is causal since F only depends on their positive part. The null map U ≡ 0 satisfies

FU (x) = -a(x) < 0, hence is a (strict) sub-solution. Finally, if U is a super-solution and if ε > 0, then F((1 + ε)U ) = (1 + ε) 2 FU + ((1 + ε) 2 -1)a ≥ 2εa > 0 pointwise, by homogeneity of F , hence (1 + ε)U is a strict super-solution.
Each of the schemes introduced in §2 can be written in the form H of Proposition 3.6. Note that a slight reformulation is required for the Reeds-Shepp models H RS+ ε,h U (p), H RS± ε,h U (p), which involve symmetric finite differences, hence expressions of the form

1≤i≤N max {0, U (p) -U (q i ), U (x) -U (r i )} 2 = 1≤i≤N max (U (p) -U (q i )) 2 + , (U (p) -U (r i )) 2 + .
This expression can be put in the form (32) using the distributivity of the "+" operator over the "max" operator, namely a + max{b, c} = max{a + b, a + c}. For consistency of (33, right) with (28), one must choose a(x) = 1 2 α(x) 2 .

Construction of a continuous super-solution. In this paragraph, we construct a supersolution to the generalized eikonal PDE (8), see §3.2 for more details on this concept. It is sampled on the discrete domain Ω h in the next paragraph, which yields a discrete super-solution bounded independently of the gridscale h > 0. In the rest of this subsection, we let for all p = (x, θ) ∈ M

u(p) := α 0 + n(θ 0 ), x + ξ 2 d S 1 (θ, θ 0 ) 2 , ( 34
)
where α 0 = C WS + max{ p ; p ∈ Ω}. (In other sections of this paper, the symbol u still stands for the value function defined by (7).) The angle θ 0 ∈ S 1 in (34) is arbitrary but fixed, and d S 1 denotes the distance function on

S 1 := R/(2πZ). For all p = (x, θ) ∈ M such that d S 1 (θ, θ 0 ) = π, one has du(p) = (n(θ 0 ), ξϕ) ∈ E * , (35)
where ϕ ∈] -π, π[ is the unique element congruent with θ -θ 0 modulo 2π. The function u is not differentiable where d S 1 (θ, θ 0 ) = π, but by convention we still define du(p) as ( 35) with ϕ := π.

As shown in the following lemma, this expression defines a super-gradient of u.

Lemma 3.7. Let u : M → R be defined by (34). Then for all p = (x, θ) ∈ M and all ṗ = ( ẋ, θ) ∈ E one has

u(p + ṗ) ≤ u(p) + du(p), ṗ + ξ 2 θ2 . Proof. Let ϕ ∈] -π, π] be congruent with θ -θ 0 modulo 2π. Then d S 1 (θ + θ, θ 0 ) 2 ≤ |ϕ + θ| 2 = ϕ 2 + 2ϕ θ + θ2 = d S 1 (θ, θ 0 ) 2 + 2 ϕ, θ + θ2
, which implies the announced result.

The next proposition lower bounds the dual metric applied to the differential of (34).

Proposition 3.8. Let F be F RS+ , F RS± , F EM , or F D . Then the dual metric (29) obeys F * p (du(p)) ≥ c 0 for any p ∈ M, where c 0 > 0 is an absolute constant.

Proof. Let p = (x, θ) ∈ Ω and let ϕ ∈] -π, π] be congruent with θ -θ 0 modulo 2π. Then

F RS+ * p (du(p)) 2 = n(θ), n(θ 0 ) 2 + + |θ -θ 0 | 2 = (cos ϕ) 2 + + ϕ 2 , F RS± * p (du(p)) 2 = n(θ), n(θ 0 ) 2 + |θ -θ 0 | 2 = (cos ϕ) 2 + ϕ 2 , 2F EM * p (du(p)) = n(θ), n(θ 0 ) + n(θ), n(θ 0 ) 2 + |θ -θ 0 | 2 = cos ϕ + cos 2 ϕ + ϕ 2 , F D * p (du(p)) = ( n(θ), n(θ 0 ) + |θ -θ 0 |) + = (cos ϕ + |ϕ|) + .
The Proposition 3.8 implies that λu is a super-solution of the HJB PDE (8), in the sense of Definition 3.14 below, where λ = α ∞ /c 0 . Indeed, this follows from the positive 1-homogeneity of F * , and the non-negativity of u on ∂Ω. We do not directly use this fact in this subsection, since we aim at constructing a discrete super-solution, but it explains the particular role played by u.

Discretization of the super-solution. We construct a discrete super-solution to the problem (28) by sampling the continuous one constructed in Proposition 3.8. This implies the existence of a bounded solution to our PDE discretization, see Corollary 3.11. For that purpose, a preliminary estimate on the discrete Hamiltonian regularity is required.

Lemma 3.9. Let F := √ 2H where H is H RS+ ε,h , H RS± ε,h , H EM ε,K,h or H D ε,h
, and where ε ≤ 1. Let p 0 ∈ Ω h , and let V, W :

M h → R obey for all ṗ ∈ M h V (p 0 + ṗ) ≥ V (p 0 ) -ṗ 2 , W (p 0 + ṗ) ≥ W (p 0 ) -ṗ .
Then for some C = C 0 (1 + ξ -1 ), where C 0 is an absolute constant, one has

F V (p 0 ) ≤ Ch/ε, F W (p 0 ) ≤ C.
Proof. For any ṅ ∈ R d one obtains using Proposition 2.2

h -2 ė∈Z d ρ ε ė( ṅ)(W (p 0 ) -W (p 0 -h ė)) 2 + ≤ h -2 ė∈Z d ρ ε ė( ṅ) h ė 2 = ṅ 2 (1 + (d -1)ε 2 ) (36)
In the case of V , observing that

V (p 0 ) -V (p 0 -h ė) ≤ h ė 2 ≤ (C WS h/ε) h
HV (p 0 ) ≤ 1 2 (Ch/ε) 2 (resp. HW (p 0 ) ≤ 1 2 C 2 ).
It follows from the above estimates since the discretized Hamiltonians of interest are sums of expressions like (36, left), with ṅ = O(1 + ξ -1 ). In the case of H EM ε,K,h one must additionally observe that K k=1 α k = 2 is bounded independently of K, see (23).

In the following two results, the notation C 1 + ξ α , where α ∈ R, means that C depends only on the parameter ξ and satisfies C ≤ C (1 + ξ α ) for all ξ ∈]0, ∞[, where C is an absolute constant.

Proposition 3.10.

Let F := √ 2H where H is H RS+ ε,h , H RS± ε,h or H D ε,h
, and where 0 < h ≤ ε ≤ 1. Define U : Z h → R by U (p) = u(p) for all p ∈ Ω h , and U (p) = 0 for all p ∈ ∂Ω h , where u is defined in (34). Then

F U (p) ≥ c 0 -C 0 h/ε (37)
for all p ∈ Ω h , where C 0 1 + ξ and c 0 is from Proposition 3.8. In the Euler-Mumford case, the constant c 0 in (37) must be replaced with c 2 0 -2C 1 K -2 , where C 1 1 + ξ -2 is the constant from Proposition 2.10.

Proof. In this proof we let p 0 = (x 0 , θ 0 ) ∈ Ω h be fixed, p = (x, θ) ∈ Ω h be an arbitrary point, and denote by θ ∈] -π, π] the unique angle congruent to θ -θ 0 . Let U : M h → R be defined by U (p) := u(p 0 ) + du(p 0 ), (x -x 0 , θ) . By consistency of the discretization, see Propositions 2.5, 2.10, 2.13, and by Proposition 3.8 we obtain, denoting by H the Hamiltonian of the model

HU (p 0 ) ≥ H p 0 (du(p 0 )) ≥ 1 2 c 2 0 , (38)
except in the Euler-Mumford case, where

HU (p 0 ) ≥ 1 2 c 2 0 -C 1 K -2 , with C 1 1 + ξ -2 . Thus F U (p 0 ) ≥ c 0 , or in the Euler-Mumford case F U (p 0 ) ≥ √ c 0 -2C 1 K -2 . For any p ∈ Ω h one has U (p) + ξ 2 θ2 ≥ U (p)
, by Lemma 3.7, with the above notation for θ ∈] -π, π]. On the other hand, if p ∈ ∂Ω h is within distance C WS of Ω h , then U (p) ≥ 0 = U (p), by choice of α 0 in the definition of u, see (34). Hence denoting V := U -U we obtain

V (p 0 + ṗ) = U (p 0 + ṗ) -U (p 0 + ṗ) ≥ -ξ 2 ṗ 2 for any ṗ ∈ hZ 3 such that ṗ ≤ C WS . Therefore F V (p 0 ) ≤ ξC 2 (1 + ξ -1 )h/ε = C 2 (1 + ξ)h/ε,
by Lemma 3.9 and the 1-homogeneity of F , where C 2 is independent of ξ. We used the fact that the expression of F V (p 0 ) only involves points within distance C WS h/ε ≤ C WS of p 0 , see Propositions 2.5, 2.10, 2.13. Using the triangular inequality F U = F (U + V ) ≤ F U + F V , pointwise on Ω h , see Proposition 3.6, we thus obtain

F U (p 0 ) ≥ F U (p 0 ) -F V (p 0 ) ≥ c 0 -C 2 (1 + ξ)h/ε.
In the Euler-Mumford case the constant c 0 in the above expression must be replaced with √ c 0 -2C 1 K -2 , see (38). The result follows.

Corollary 3.11. Assume that 0 < h ≤ ε ≤ 1 and ε ≤ K 0 h, where K 0 1 + ξ. Let H be the discretized Hamiltonian H RS+ ε,h , H RS± ε,h or H D ε,h . Then the system of equations HU (p) = α(p) 2 /2 for all p ∈ Ω h , and U (p) = 0 for all p ∈ ∂Ω h , admits an unique solution U : M h → R, which is bounded independently of h and ε.

The same holds for the Hamiltonian

H EM ε,K,h , provided K ≥ K 1 1 + ξ -1 .
Proof. Assume that the parameters ε, h and K obey the above constraints, with the constants K 0 = 2C 0 and K 1 = 4 √ C 1 , where C 0 and C 1 are from Proposition 3.10. Then the considered system of equations admits the super-solution λU , where U is defined in Proposition 3.10 and λ = 2 α L ∞ /c 0 . Applying Theorem 3.5, whose assumptions were established in Proposition 3.6 for the models of interest, except for the existence of a super-solution which we just proved, we conclude the proof.

Discontinuous viscosity solutions of eikonal PDEs

In this subsection, we establish Theorem 3.1 using the theory of discontinuous solutions to static first order Hamilton-Jacobi-Bellman PDEs. We rely on the framework of section V of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], intended for time optimal control problems without controllability, either local or global. Let F be F RS+ , F RS± , F EM or F D , and consider the PDE.

∀p ∈ Ω, F * p (du(p)) = α(p), ∀p ∈ ∂Ω, u(p) = 0 or F * p (du(p)) = α(p), (39)
where F * refers to the dual metric, defined on the co-tangent space M × E * , see (29). In comparison with our initial formulation (8), two remarks are in order: (I) we used the relation H = 1 2 (F * ) 2 relating the Hamiltonian H with the dual metric F * (29) to reformulate the HJB PDE in Ω, and (II) we emphasized that boundary conditions must be interpreted in a relaxed sense, following the notations of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], due to the possible discontinuity of the solution. The study of solutions to this system relies on one sided notions of continuity. Definition 3.12. Let (X, d) be metric space, and let u : X → R. The function u is said Lower-Semi-Continuous (resp. Upper-Semi-Continuous) iff for any converging sequence p

n → p ∈ X lim inf n→∞ u(p n ) ≥ u(p) (resp. lim sup n→∞ u(p n ) ≤ u(p).)
The acronyms used for these properties (resp. and also Boundedness) are USC and LSC (resp. BUSC and BLSC). Recall that the cartesian grid of scale h is denoted M h ⊆ M, see (40), and that B(p, r) ⊆ M is the Euclidean ball of center p ∈ M and radius r > 0.

Lemma 3.13. For each n ≥ 0 let h n > 0 and U n : M hn → R. Assume that h n → 0 as n → ∞, and that U n is uniformly bounded independently of n. Then u, u : M → R defined as follows are respectively BUSC and BLSC

u(p) := lim r→0 lim inf n→∞ inf q∈M hn q-p <r U n (q), u(p) := lim r→0 lim sup n→∞ sup q∈M hn q-p <r U n (q). ( 40
)
Proof. We focus on the case of u, since the case of u is similar. First note that u(p) is well defined and uniformly bounded w.r.t. p ∈ M. Indeed (i) for any fixed r > 0, and for sufficiently large n, the set M hn ∩ B(p, r) is non-empty since h n → 0 as n → ∞, (ii) U n is uniformly bounded, and (iii) the leftmost limit as r → 0 is monotone, namely increasing as r decreases to 0, hence well defined.

In order to establish the Lower Semi-Continuity of u, let us consider an arbitrary sequence (p n ) n≥0 converging to p ∈ M. Let also r n → 0 and ε n → 0 be vanishing sequences of positive reals. By construction, for any n ≥ 0 there exists ϕ(n) ≥ n and

q n ∈ B(p n , r n ) such that U ϕ(n) (q n ) ≤ u(p n ) + ε n . Then, as announced, since q n → p and ϕ(n) → ∞ as n → ∞, lim inf n→∞ u(p n ) ≥ lim inf n→∞ U ϕ(n) (q n ) ≥ u(p).
Following [START_REF] Bardi | On Hopf's formulas for solutions of Hamilton-Jacobi equations[END_REF], we introduce the concept of sub-and super-solutions to the system (39). Definition 3.14.

• A sub-solution of (39) is a BUSC u : Ω → R such that for any p ∈ Ω and any ϕ ∈ C 1 (Ω, R) for which u -ϕ attains a local maximum at p, one has:

p ∈ Ω ⇒ F * p (dϕ(p)) ≤ α(p), p ∈ ∂Ω ⇒ min{u(p), F * p (dϕ(p)) -α(p)} ≤ 0.
• A super-solution of (39) is a BLSC u : Ω → R such that for any p ∈ Ω and any ϕ ∈ C 1 (Ω, R) for which u -ϕ attains a local minimum at p, one has:

p ∈ Ω ⇒ F * p (dϕ(p)) ≥ α(p), p ∈ ∂Ω ⇒ max{u(p), F * p (dϕ(p)) -α(p)} ≥ 0.
It known that replacing "local maximum" with "strict global maximum" (resp. "local minimum" with "strict global minimum") in Definition 3.14 yields an equivalent definition, see [START_REF] Bardi | On Hopf's formulas for solutions of Hamilton-Jacobi equations[END_REF].

Following a classical strategy [START_REF] Bornemann | Finite-element Discretization of Static Hamilton-Jacobi Equations based on a Local Variational Principle[END_REF], we estimate in the next lemma the discretization error of our numerical scheme when applied to continuously differentiable functions, and conclude in the subsequent proposition that suitable limits of solutions to our discrete numerical schemes are sub-and super-solutions to the HJB PDE (39).

Lemma 3.15. Let F ε,h := 2H ε,h where H ε,h is H RS+ ε,h , H RS± ε,h or H D ε,h and 0 < h ≤ ε ≤ 1. Let ϕ ∈ C 1 (M,
R), and let ω be the modulus of continuity of dϕ. Then

|F ε,h ϕ(p) -F * p (dϕ(p))| ≤ C(ω(C WS h/ε) + ε dϕ(p) ). ( 41
)
In the Euler-Mumford case,

|F ε,K,h ϕ(p) -F * p (dϕ(p))| ≤ C(ω(C WS h/ε) + (ε + K -1 n ) dϕ(p) ).
Proof. In this proof, if p = (x, θ), q = (x , θ ) ∈ M = R 2 × S 1 , then qp (abusively) stands for (x -x, ϕ) ∈ E := R 2 × R where ϕ ∈] -π, π] is congruent with θ -θ. Fix the point p ∈ M, and define the tangent map Φ : q ∈ M → ϕ(p) + dϕ(p), qp . Then |ϕ(q) -Φ(q)| ≤ pq ω( pq ) for any q ∈ M. This implies F (ϕ -Φ)(p) ≤ Cω(C WS h/ε) where F := F ε,h and C 1 + ξ -1 , by Lemma 3.9, the 1-homogeneity of F , and since any point q appearing in the expression of F U (p) satisfies pq ≤ C WS h/ε. Proceeding likewise for F (Φ -ϕ)(p) and using the triangular inequality, proved for F in Proposition 3.6, we obtain

|F ϕ(p) -F Φ(p)| ≤ max{F (Φ -ϕ)(p), F (ϕ -Φ)(p)} ≤ C ω(C WS h/ε).
The second contribution to (41) comes from the estimate

|F Φ(p) -F * p (dϕ(p))| = √ 2| HΦ(p) -H p (dϕ(p))| ≤ √ 2 |HΦ(p) -H p (dϕ(p))|,
where we used the classical inequality | √ a -√ b| ≤ |a -b| for any a, b ≥ 0. Inserting the discretization error of the Hamiltonian, see Propositions 2.5, 2.10 and 2.13, we conclude the proof.

Proposition 3.16. Let u, u : M → R be defined as in Theorem 3.1 (resp. Theorem 3.2). Then u is a super-solution, and u is a sub-solution, in the sense of Definition 3.14.

Proof. We focus on the case of u, since the case of u is similar. By lemma 3.13, u is BUSC as required. Let p ∈ Ω and ϕ ∈ C 1 (Ω, R) be such that u attains a strict global maximum at p.

For each n ≥ 0, define X n := {q ∈ M hn ; d(q, Ω) ≤ C WS h n }, where d denotes the Euclidean distance, and let p n ∈ X n be a point where U n -ϕ attains its global maximum. Then U n (p n )ϕ(p n ) → u(p)-ϕ(p) as n → ∞, up to extracting a sub-sequence (because of the lim inf operator in (40)). This implies p n → p as n → ∞, by strictness of the maximum of u-ϕ at p. In addition

U n (p n ) = (U n (p n ) -ϕ(p n )) + ϕ(p n ) → (u(p) -ϕ(p)) + ϕ(p) = u(p)
as n → ∞, by choice of the sequence (p n ) n≥0 and by continuity of ϕ. In order to conclude the proof, we distinguish wether the test point p ∈ Ω lies in the interior or on the boundary of the domain.

Case where p ∈ Ω. By construction of p n , one has

U n (p n ) -U n (q) ≥ ϕ(p n ) -ϕ(q) for all q ∈ X n . Hence α(p n ) = F U n (p n ) ≥ F ϕ(q n ) ≥ F * pn (dϕ(p n )) -C(ω(C WS h n /ε n ) + ε n dϕ(p n ) ),
where the first inequality is by monotony of the discretized operator F := F ε,h , and the second one is by Lemma 3.15. Thus α(p) = lim α(p n ) ≥ lim F * pn (dϕ(p n )) = F * p (dϕ(p)) as n → ∞. Case where p ∈ ∂Ω. We distinguish two sub-cases: if p n ∈ Ω for infinitely many integers n ≥ 0, then F * p (dϕ(p)) ≤ α(p) as before, thus the (relaxed) boundary condition is satisfied as desired. Otherwise, up to extracting a subsequence, one has 0 = U n (p n ) → u(p) as n → ∞, thus u(p) = 0 and the boundary condition is satisfied.

The following result concludes the proof of Theorems 3.1 and 3.2. Proposition 3.17 (Adapted from [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]). The value function u defined by (7) is the smallest super-solution to the HJB PDE (39), and û defined by (27) is the largest sub-solution.

The fact that u is the smallest super-solution follows from Theorem 3.7 in chapter V [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], and the fact that û is the largest super-solution from Theorem 4.29 in the same chapter. To be complete, we describe below the slight reformulation of the optimal control problem (7) required to match the notations of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], and check that the assumptions used in [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] are satisfied. For that purpose, we introduce a compact and convex set A ⊆ R × R and regard its elements a = ( ẋ, θ) ∈ A as a (scalar) physical velocity, and an angular velocity. The following instantiations of A are considered

A RS+ := {( ẋ, θ) ∈ R 2 ; ẋ2 + θ2 ≤ 1, ẋ ≥ 0}, A RS± := {( ẋ, θ) ∈ R 2 ; ẋ2 + θ2 ≤ 1}, ( 42
)
A EM := {( ẋ, θ) ∈ R 2 ; ( ẋ -1/2) 2 + θ2 ≤ (1/2) 2 }, A D := {( ẋ, θ) ∈ R 2 ; 0 ≤ | θ| ≤ ẋ}. Define f : M × A → E and l : M × A → R as f ((x, θ), ( ẋ, θ)) := -( ẋn(θ), θ)/α(x, θ), l((x, θ), a) := 1, ( 43
)
where α : M →]0, ∞[ is the local cost function. Let also T := M \ Ω be the target region. Consider, following [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], the optimal control problem v(p) := inf (a,T )∈A T 0 l(q p,a (t), a(t)) e -t dt, subject to qp,a (t) = f (q p,a (t), a(t)), ∀t ∈ [0, T ],

q p,a (0) = p and q p,a (T ) ∈ T ,

where A consists of all pairs of (free) terminal times T ≥ 0 and measurable function a : [0, T ] → A. We claim that, with the choices (42), one has T 0 e -t dt = 1 -exp(-u(p)). Conversely, admissible paths for (44) can be reparametrized into admissible paths for (7), and the identity (45) follows. Similarly, one can define v by replacing T with its interior in (44), and obtain that v = 1 -exp(-û) in Ω.

v(p) = u(p) 0 e -t dt = 1 -exp(-u(p)), (45)
The following PDEs are thus equivalent, the rightmost being the one considered in [BCD08]

H p (du(p)) = 1 2 α(p) 2 , F * p (du(p)) = α(p), v(p) + F p (dv(p)) -α(p) = 0,
where p ∈ Ω is arbitrary. The main reason why [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] considers v = 1 -exp(-u) instead of u is that v remains bounded even for problems lacking global controllability, whereas u may take infinite values. This technicality is however irrelevant for the problems considered in this paper, since global controllability does hold as shown in §3.1. Finally, we check the specific assumptions of Theorems 3.7 and 4.29 in chapter V of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF]. These are (I) the compactness of the set A of controls (42), (II) the Lipschitz continuity of f and l, see (43), which follows from the Lipschitz continuity of the cost α, and its boundedness below on the compact domain Ω, and (III) the closedness of T , following from the openness of Ω.

Continuity properties of the value function

This subsection is devoted to the proof of Proposition 3.3, which describes sets where the two variants u, û of the value function coincide. The announced properties of the Reeds-Shepp reversible and the Dubins models follow from rather general arguments. Point (I) indeed follows from the locally controllability of the Reeds-Shepp reversible model, which is due to its subriemannian structure and to Chow's theorem, see [START_REF] Montgomery | A Tour of Subriemannian Geometries, Their Geodesics and Applications[END_REF] and 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]. Hence u = û and the value function u is in fact not only continuous, but 1/2-Holder continuous in this case. Point (IV) on the Dubins model follows from another general argument of [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF], involving Baire's theorem, see Lemma 3.25. In contrast the proof of points (II) and (III) on the Reeds-Shepp forward and Euler-Mumford models requires a geometrical pertubation argument, which involves lifting diffeomorphisms from R 2 to M := R 2 × S 1 , as presented below.

The argument arg( ẋ) of a non-zero vector ẋ ∈ R 2 is defined as the unique angle θ ∈ S 1 such that ẋ = ẋ n(θ). The matrix vector product is denoted by " • ". Lemma 3.18. Let ψ be a C n diffeomorphism of R 2 , where n ≥ 2. Define Ψ(x, θ) = (y, ϕ) by

y := ψ(x), ϕ := arg(dψ(x) • n(θ)), ( 46
)
for all (x, θ) ∈ R 2 × S 1 . Then Ψ is a C n-1 diffeomorphism of R 2 × S 1 . Furthermore, let ( ẋ, θ) ∈ R 2 × R and let ( ẏ, φ) := dΨ(x, θ) • ( ẋ, θ). Then K ẋ ≤ ẏ ≤ K ẋ , | φ| ≤ (K| θ| + K 2 ẋ )/K, ( 47
)
where K = (dψ(x)) -1 -1 , K := dψ(x) , and K 2 := d 2 ψ(x) .

Proof. The bijectivity of θ ∈ S 1 → ϕ := arg(dψ(x) • n(θ)), for any fixed x, follows from the invertibility of dψ(x). The estimate (47, left) follows from the definition of the operator norm A := sup ẋ =0 A ẋ / ẋ of a matrix A. The upper bound on | φ| is obtained by composing the following two Taylor expansions: the first one

arg(x + ẋ) = arg(x) + x -2 x ⊥ , ẋ + o( ẋ ),
is obtained by basic geometric reasoning, and the second one

dψ(x + ẋ) • n(θ + θ) = dψ(x) • n(θ) + (d 2 ψ(x) • ẋ) • n(θ) + θdψ(x) • n(θ) ⊥ + o( ẋ + | θ|),
by bi-linearity of the matrix-vector product.

The next lemma upper bounds the composition of the metric, of the models RS+, RS± and EM, with the tangent map to a diffeomorphism of the form (46). No similar estimate can be established for the Dubins metric, due to the hard constraint |ξ θ| ≤ ẋ appearing in (26).

Lemma 3.19. Under the assumptions of Lemma 3.18, denoting p := (x, θ), ṗ := ( ẋ, θ), q := (y, ϕ) and q := ( ẏ, φ), one has

F RS+ q ( q) ≤ K RS F RS+ p ( ṗ), F RS± q ( q) ≤ K RS F RS± p ( ṗ), F EM q ( q) ≤ K EM F EM p ( ṗ), with K 2 RS := max{K 2 +ξ 2 ε(1+ε), (1+ε)K} and K EM := max{K +ξ 2 ε(1+ε)/K, (1+ε)K 2 /K},
where ε := K 2 /K and K := K/K.

Proof. By construction of the diffeomorphism Ψ, the colinearity constraint involved in the definition (3) of the metrics is preserved: ẋ = ẋ n(θ) ⇒ ẏ = ẏ n(ϕ), and likewise for the unsigned colinearity constraint ẋ = ẋ, n(θ) n(θ) ⇒ ẏ = ẏ, n(ϕ) n(ϕ). Using the inequality (a + εb) 2 ≤ (1 + ε)(a 2 + εb 2 ), valid for any a, b ∈ R, ε > 0, we obtain as announced

ẏ 2 + ξ 2 φ2 ≤ (K ẋ ) 2 + ξ 2 (K θ + ε ẋ ) 2 ≤ K 2 + ξ 2 ε(1 + ε) ẋ 2 + (1 + ε)K 2 ξ 2 θ2 ẏ + ξ 2 φ2 ẏ ≤ K ẋ + ξ 2 (K θ2 + ε ẋ ) 2 K ẋ ≤ (K + ξ 2 ε(1 + ε)/K) ẋ + (1 + ε) K 2 K ξ 2 θ2 ẋ .
The previous lemma is next specialized to diffeomorphisms defined by the flow of a vector field.

Corollary 3.20. Let ν : R 2 → R 2 be a vector field with continuous and uniformly bounded first and second derivatives. Let (ψ t ) t≥0 be the family of C 2 diffeomorphisms of R 2 defined by

ψ 0 (x) = x, d dt ψ t (x) = ν(ψ(x)), (48)
for any x ∈ R 2 . Then for any path

γ ∈ Lip([0, 1], R 2 × S 1 ) one has length F (Ψ t • γ) ≤ (1 + Ct) length F (γ) where C = C(ν, ξ)
, where Ψ t is defined from ψ t by (46), and where the metric

F is F RS+ , F RS± or F EM .
Proof. Define for any t ≥ 0 the constants

K(t) := inf x∈R 2 (dψ t (x)) -1 -1 , K(t) := sup x∈R 2 dψ t (x) , K 2 (t) := sup x∈R 2 d 2 ψ t (x) .
The Taylor expansion ψ t (x) = x + tν(x) + o(t) can be differentiated twice w.r.t. x ∈ R 2 since the two sides are C 2 smooth. This yields dψ t (x) = Id +tdν(x) + o(t) and d

2 ψ t (x) = td 2 ν(x) + o(t). Hence K(t) = 1 + O(t), K(t) = 1 + O(t) and K 2 (t) = O(t).
Therefore the corresponding constants of Lemma 3.19 obey K RS (t) = 1 + O(t) and K EM (t) = 1 + O(t). Inserting these estimates into the path length expression (5), and using the Lipschitz regularity of the cost function α, we obtain the announced result.

The following lemma further specializes the diffeomorphisms considered, which are designed so as to offset one endpoint of a given path, and leave the other endpoint unaffected.

Lemma 3.21. Let γ : [0, 1] → R 2 × S 1 be a Lipschitz path. Denote γ(0) = (x 0 , θ 0 ), γ(1) := (x 1 , θ 1 ), and assume that x 0 = x 1 . Let F be F RS+ , F RS± or F EM and let us assume that length F (γ) < ∞. Then there exists a family of diffeomorphisms ψ ṗ, ṗ ∈ R 2 × R such that: for any sufficiently small ṗ length

F (ψ ṗ • γ) ≤ (1 + O( ṗ )) length F (γ), ψ ṗ(γ(0)) = γ(0) + ṗ, ψ ṗ(γ(1)) = γ(1).
Proof. Denote r := x 0 -x 1 > 0. Let ν 1 , ν 2 , ν 3 : R 2 → R 2 be smooth vector fields supported on B(x 0 , r/2), and defined by ν 1 (x) = (1, 0), ν 2 (x) = (0, 1), and ν 3 (x) = x 0 + (x -x 0 ) ⊥ respectively for all x ∈ B(x 0 , r/4), where B(x, r) ⊆ R 2 denotes the Euclidean unit ball centered at x ∈ R 2 and of radius r > 0. Let ψ ν t denote the diffeomorphism generated by the flow at time t of a vector field ν, as defined in (48).

By construction, one has

ψ ν 1 t 1 (x, θ) = (x+(t 1 , 0), θ), ψ ν 2 t 2 (x, θ) = (x+(0, t 2 ), θ)
, and ψ ν 3 t 3 (x, θ) = (x 0 +R t 3 (x-x 0 ), θ+t 3 ) for any sufficiently small t 1 , t 2 , t 3 ∈ R, for any x ∈ R 2 sufficiently close to x 0 , and for any θ ∈ R. Also For that purpose, we need to distinguish a particular class of degenerate paths: we say that γ : [0, 1] → R 2 × S 1 is a purely angular motion iff it has the form γ(t) = (x, θ(t)), for all t ∈ [0, 1], where x ∈ R 2 is a constant and θ : [0, 1] → S 1 is an arbitrary function.

ψ ν 1 t 1 = ψ ν 2 t 2 = ψ ν 3 t 3 = Id on a neighborhood of x 1 , for any t 1 , t 2 , t 3 ∈ R. Defining ψ ṗ := ψ ν 1 t 1 • ψ ν 2 t 2 • ψ ν 3 t 3 , where ṗ = (t 1 , t 2 , t 3 ) ∈ R 2 × R,
Proposition 3.22. Let p ∈ Ω be such that the minimal path γ : [0, 1] → Ω for u(p) is not a purely angular motion. Assume also that int(Ω) = Ω, and the model is Reeds-Shepp forward or Euler-Mumford. Then u(p) = û(p).

Proof. As noted in the introduction, a minimal path γ exists by Appendix A of [START_REF] Chen | A New Finsler Minimal Path Model With Curvature Penalization for Image Segmentation and Closed Contour Detection[END_REF] (it may be non-unique). Since γ is not a purely angular motion, there exists t * ∈]0, 1] such that denoting p 0 = (x 0 , θ 0 ) = γ(0) and p 1 = (x 1 , θ 1 ) = γ(t * ) one has x 0 = x 1 . Denote γ 0 := γ |[0,t * ] and γ 1 := γ |[t * ,1] . For any sufficiently small ṗ ∈ R 2 × R we may construct an admissible path from p 0 + ṗ to p = γ(1), by concatenation of ψ ṗ • γ 0 with γ 1 , where ψ ṗ is as defined in Lemma 3.21. If

p 0 + ṗ ∈ M \ Ω this implies û(p) ≤ length F (ψ ṗ • γ 0 ) + length F (γ 1 ) ≤ (1 + O( ṗ )) length F (γ) = (1 + O( ṗ ))u(p).
The assumption int(Ω) = Ω implies that one can find arbitrarily small ṗ ∈ E such that Proof. We only need to show that purely angular motions cannot be minimizing for the optimal control problem (7) under these assumptions, and apply Lemma 3.21. In the case where Ω = U × S 1 , the domain shape forbids that a purely angular motion has one endpoint in ∂Ω and the other in Ω, hence it is not even an admissible candidate path. In the Euler-Mumford case, any non-constant purely angular motion has infinite length, hence it cannot be minimizing.

p 0 + ṗ ∈ M \ Ω, hence û(p) ≤ u(p).
The continuity of the value function(s) u and û at the points where u = û, announced in Proposition 3.3, follows the next elementary lemma. Recall that u and û are respectively BLSC and BUSC, since they are sub-and super-solutions of (39), see Proposition 3.17. Lemma 3.24. Let (X, d) be a metric space, and let u, û : X → R be respectively LSC and USC, and obey u ≤ û on X. If p ∈ X is such that u(p) = û(p), then u and û are continuous at p.

Proof. Let (p n ) n≥0 be an arbitrary sequence converging to p ∈ X. Using successively that u is LSC, that u ≤ û, and that û is USC, one obtains

u(p) ≤ lim inf n→∞ u(p n ) ≤ lim sup n→∞ û(p n ) ≤ û(p).
If u(p) = û(p) then the sequences u(p n ) and û(p n ) must be converging to this common limit, which implies the announced continuity of u and û at p. Finally, we discuss the case of the Dubins model. By Corollary 4.30 in section V [START_REF] Bardi | Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations[END_REF] (which uses in particular the assumption int(Ω) = Ω), one has û * = u ≤ û, where v * (resp. v * ) denotes the lower (resp. upper) semi-continuous envelope of a function v. Hence u = û at each point of continuity of û. The next elementary lemma shows that these points are dense, which concludes the proof of Proposition 3.3. Recall that a residual set is an intersection of open dense sets, i.e. the complement of a meager set, and that by Baire's theorem any residual subset of a complete metric space is dense. Lemma 3.25. Let (X, d) be a metric space, and let u : X → R be BLSC or BUSC. Then the points of continuity of u form a residual set.

Proof. We focus on the BLSC case, and note that the BUSC case follows by considering -u. Let u * be the upper semi-continuous envelope of u, and for each n ≥ 1 let C n := {x ∈ X; u * (x)u(x) < 1/n}. For any x ∈ X, one has the equivalences:

u is continuous at x ⇔ u * (x) = u(x) ⇔ x ∈ ∩ n≥1 C n .
In the following, we fix n ≥ 1 and establish that C n is open and dense, hence ∩ n≥1 C n is a residual set which concludes the proof. The openness of C n follows from the upper semi-continuity of u * -u. Assume for contradiction that C n is not dense, hence that there exists an open ball B(x 0 , r) on which u * -u ≥ 1/n identically. Let also M := sup{u(x) -u(y); x, y ∈ X}, which is finite by assumption.

Let K > 2nM be an integer. We construct inductively a sequence (x k ) 0≤k≤K by choosing x k+1 such that d(x k , x k+1 ) < r/K and u(x k+1 ) > u(x k ) + 1/(2n). This is possible since

d(x k , x 0 ) < kr/K ≤ r, hence x k / ∈ C n , thus u(x k ) + 1/(2n) < u * (x k ) = lim sup y→x k u(y). Eventually we obtain u(x K ) -u(x 0 ) = K-1 i=0 u(x i+1 ) -u(x i ) > K/(2n) ≥ M ,
which contradicts the definition of M . Hence C n must be dense, and the proof is complete.

Basis reduction techniques

This section is devoted to the proof of Proposition 2.2, using techniques from lattice geometry. For that purpose, we study an optimization problem referred to as Voronoi's first reduction of quadratic forms [START_REF] Schürmann | Computational geometry of positive definite quadratic forms[END_REF], using special coordinate systems known as obtuse superbases of lattices [START_REF] Conway | Low-Dimensional Lattices. VI. Voronoi Reduction of Three-Dimensional Lattices[END_REF]. Similar techniques are used for the discretization on cartesian grids of anisotropic diffusion in [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF][START_REF] Fehrenbach | Sparse non-negative stencils for anisotropic diffusion[END_REF], of Monge-Ampere equations in [START_REF] Benamou | Monotone and Consistent discretization of the Monge-Ampere operator[END_REF], and of anisotropic eikonal equations in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF].

We first need to introduce some notation. Let d ≥ 1 be the ambient dimension, d ∈ {2, 3} for the applications intended in this paper. The canonical d-dimensional Euclidean space, and d-dimensional integer lattice are denoted

E d := R d , L d := Z d .
The

P := {D ∈ S ++ (E * d ); ∀v ∈ L * d \ {0}, v D ≥ 1}, L(M ) := inf D∈P Tr(M D). (49)
Denote by M, D := Tr(M D) the duality bracket between S(E) and S(E * ), and observe that v 2 D = D, v ⊗ v , where v ⊗ v ∈ S + (E d ) denotes the outer product of a co-vector v ∈ E * d with itself. Note also that the semi-definite constraint D 0 in (49, left) is redundant. One can therefore rephrase Voronoi's optimization problem L(M ) to make its linear structure more apparent:

minimize

M, D subject to D, v ⊗ v ≥ 1 for all v ∈ L * d \ {0}. (50)
The linear program (50) was shown by Voronoi to be feasible in arbitrary dimension d ≥ 1, in the sense that L(M ) has a non-empty and compact set of minimizers for any M ∈ S ++ (E d ).

The next lemma describes the decomposition of a symmetric matrix M using the directions vij ∈ L * d attached to a superbase. As one can suspect, it is related to the KKT relations (51), see Proposition 4.6 below.

Lemma 4.4. Let M ∈ S(E d ), and let ( ė0 , • • • , ėd ) be a superbase of Z d . Then

M = i<j λ ij vij ⊗ vij ,
where

λ ij := -ėi , M ėj (54)
Proof. For any i, j, k, l ∈ 0, d such that i = j and k = l, one has vij , ėk vij , ėl = -1 if {i, j} = {k, l}, and 0 otherwise. Hence denoting by M the r.h.s. of (54) one has ėk , M ėl = -λ kl = ėk , M ėl for all k = l. Equality also holds if k = l, using the identity ėk =i =k ėi . Since ( ė1

, • • • , ėd ) is a basis of E d , this implies M = M as announced.
The following lemma attaches to any superbase b of L d a vertex D b of the polytope P, defined in (49). Vertices of P are referred to as perfect forms, and play a central role in lattice geometry [START_REF] Schürmann | Computational geometry of positive definite quadratic forms[END_REF]. Interestingly, Voronoi proved that the polytope P has only finitely many distinct equivalence classes of vertices under the action of linear changes of coordinates preserving the lattice L d , see [START_REF] Schürmann | Computational geometry of positive definite quadratic forms[END_REF]

distinct i 0 , • • • , i k ∈ 0, d one has I⊆ 0,d δ I i 0 • • • δ I i k = #{I ⊆ 0, d ; {i 0 , • • • , i k } ⊆ I} = 2 d-k .
Hence for any pairwise distinct i, j, k, l ∈ 0, d , noting that vij , ėI = δ I i -δ I j for any I ⊆ 0, d ,

2 d vij , vik D b = I⊆ 0,d (δ I i -δ I j )(δ I i -δ I k ) = I⊆ 0,d (δ I i -δ I i δ I j -δ I i δ I k +δ I j δ I k ) = 2 d -2 d-1 -2 d-1 +2 d-1 .
This establishes vij , vik

= i<j λ ij ( D, vij ⊗ vij - D b , vij ⊗ vij ) = i<j λ ij ( vij 2 D -1) ≥ 0. The result follows.
The previous proposition leaves open the question of the existence of an M -obtuse superbase b, given M ∈ S ++ (E d ). In dimension d ∈ {2, 3}, such a b always exists, and can be constructed via Selling's algorithm [START_REF] Selling | Ueber die binären und ternären quadratischen Formen[END_REF] which is implemented "as is" in our numerical experiments §5. Note that this algorithm could in principle be accelerated by a preliminary basis reduction step [START_REF] Phong | Low-dimensional lattice basis reduction revisited[END_REF], but the advantage is only visible for extremely large condition numbers, which are irrelevant for applications to anisotropic PDEs. A closer inspection shows that Selling's algorithm is equivalent to the Simplex algorithm applied to the linear program (50). Selling's algorithm is also described in Appendix B of [START_REF] Bost | Hermitian vector bundles and extension groups on arithmetic schemes. I. Geometry of numbers[END_REF]. In dimension d = 2, Selling's algorithm is equivalent to exploring the an arithmetic structure named the Stern-Brocot tree, see [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF] for details on this approach and an application to second order HJB PDEs for stochastic control.

Radius of the decomposition

We bound in this section the Euclidean norm of the integral co-vectors involved in the matrix decomposition (51). Our estimate is sharper than the one presented in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], which has significant consequences on the convergence analysis §3, see the discussion after Theorem 4.11. The results of this subsection and of the next are new to the author's knowledge. One reason In the rest of this section, we denote by C d the constant of Proposition 4.8. Our next result describes a non-negative, lattice-adapted decomposition of symmetric tensors. It coincides with the one considered in Proposition 1.1 of [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], but the stencil radius estimate is sharper: linear in the condition number of M , instead of quadratic v ≤ C Cond(M ) 2 in dimension d = 3. This improvement has theoretical consequences, since the assumption h n /ε n → 0 in Theorem 3.1 would otherwise need to be replaced with the stronger (and unrealistic in applications) assumption h n /ε 2 n → 0.

Corollary 4.12. Let M ∈ S ++ (E d ), where d ∈ {2, 3}. Then there exists non-negative weights ρ v(M ), supported on at most It is possible to construct tensor decompositions of the form (57) by other means than Voronoi's first reduction, but the resulting stencil radius is typically larger, and provably larger [START_REF] Mirebeau | Minimal stencils for discretizations of anisotropic PDEs preserving causality or the maximum principle[END_REF] in dimension d = 2. PDE schemes of small stencil size are appreciated in numerical applications for reasons of robustness, accuracy, implementation of boundary conditions, and parallelization potential.

d(d + 1)/2 elements v ∈ L * d , such that v∈L * d ρ v(M )v ⊗ v = M. (57)
∀ ṗ ∈ E d , v∈L * d ρ ε v(n) v, ṗ 2 = n, ṗ 2 + ε 2 ( ṗ 2 -n, ṗ 2 ). ( 58
) Furthermore v ≤ 2C d ε -1 whenever ρ ε v(n) > 0, and ρ ε v(n) v 2 = 1 + ε 2 (d -1).
Proof. Apply Corollary 4.12 to M = n ⊗ n + ε 2 (Id -n ⊗ n) ∈ S ++ (E d ), whose eigenvalues are 1 with eigenvector n, and ε 2 with multiplicity d -1, thus Cond(M ) = ε -1 .

Taking positive parts of linear forms

We conclude in this section the proof of Proposition 2.2. Our main technical result, presented below, shows how to use decompositions of anisotropic symmetric positive definite tensors, as in Corollary 4.13, to build approximations of positive parts of linear forms.

Proposition 4.14. Let n ∈ E * d with n = 1, and let

ε > 0. Let v1 , • • • , vK ∈ E * d and ρ 1 , • • • , ρ K ≥ 0 be such that for all ṗ ∈ E d 1≤k≤K ρ k vk , ṗ 2 = n, ṗ 2 + ε 2 ( ṗ 2 -n, ṗ 2 ). ( 59
)
Assume that5 vk , n ≥ 0 for each 1 ≤ k ≤ K. Then for all ṗ ∈ E d one has

n, ṗ 2 + ≤ 1≤k≤K ρ k vk , ṗ 2 + ≤ n, ṗ 2 + + ε 2 ( ṗ 2 -n, ṗ 2 ). ( 60
)
Proof. In this proof, we identify E d with its dual E * d thanks to the Euclidean structure, and thus drop the "dot" and "hat" superscripts on vectors and co-vectors for clarity. We assume that

ρ k = 1, for all 1 ≤ k ≤ K, up to replacing v k with √ ρ k v k . Denote by v * k := v k -v k , n n the orthogonal projection of v k on the hyperplane orthogonal to n, where 1 ≤ k ≤ K. Then by (59) 1≤k≤K n, v k 2 = 1, 1≤k≤K v * k ⊗ v * k = ε 2 (Id -n ⊗ n).
The proof of (60) is split into two parts, depending on the sign of n, p . If n, p ≤ 0, then

v k , p = v k , n n, p + v * k , p ≤ v * k , p for all 1 ≤ k ≤ K, thus as announced 1≤k≤K v k , p 2 + ≤ 1≤k≤K v * k , p 2 + ≤ 1≤k≤K v * k , p 2 = ε 2 ( p 2 -n, p 2 ).
In contrary if n, p ≥ 0, then the second inequality of (60) is immediate. In addition

v k , p 2 + ≥ v k , p 2 -v * k , p 2 for any 1 ≤ k ≤ K. (Indeed, if v k , p ≥ 0 then v k , p 2 + = v k , p 2 ≥ v k , p 2 -v * k , p 2 , and in contrary if v k , p ≤ 0 then 0 ≥ v k , p = v k , n n, p + v * k , p ≥ v * k , p , thus v k , p 2 + = 0 ≥ v k , p 2 -v * k , p 2 .) Hence, we conclude 1≤k≤K v k , p 2 + ≥ 1≤k≤K v k , p 2 -v * k , p 2 = n, ṗ 2 + ε 2 p * 2 -ε 2 p * 2 = n, p 2 ,
where we denoted p * := pp, n n, so that p * 2 = ṗ 2 -n, ṗ 2 .

Combining Corollary 4.13 with Proposition 4.14 we conclude the proof of Proposition 2.2, thanks to the following two remarks. (I) The assumption vk , n ≥ 0 for all 0 ≤ k ≤ K in Proposition 4.14 is not restrictive, since one can always replace vk with its opposite -v k without incidence on (59). (II) The roles of E d and of its dual E * d are exchanged in Proposition 2.2.

Numerical experiments

This section is devoted to a numerical validation of the proposed PDE schemes, and to an illustration of their potential applications. For validation, compare in §5.1 the geodesics obtained by solving Hamilton's ODE of geodesics against the results of our PDE discretization. A second validation is presented in §A, and based on numerical counterparts of the control sets of Figure 2. We investigate applications to motion planning in §5.2, and to vessel tracking in simulated medical data in §5.3. All the test cases presented in this paper are based on synthetic data. Experiments closer to applications and involving real data will be the subject of future work.

Free and open source numerical codes for reproducing (most of) the numerical experiments, as well as additional examples, are available on the author's webpage6 .

Comparison with geodesic shooting

dp dt = - ∂H ∂ p , dp dt = ∂H ∂p , (61)
where p(t) and p(t) are respectively the state and the co-state of the geodesic at time t ∈ [0, 1]. Given an initial point p 0 , and a target point p 1 , one can try to adjust the initial co-state p0 so that the solution to (61) starting from (p(0), p(0)) = (p 0 , p0 ) satisfies p(1) = p 1 . We tried this procedure, referred to as geodesic shooting, using a fourth order Runge-Kutta method for solving (61), and a Newton method for adjusting p0 . Another approach to the computation of geodesics is to numerically solve the PDE (8), and then extract the minimal geodesics as streamlines of the geodesic flow (11). We use a second order Euler scheme for that latter ODE, together with the following upwind estimate of the geodesic flow direction dH p (du(p)), where p ∈ M. Assume that the local Hamiltonian H p is approximated in the following form

H(p) = 1 2 1≤k≤K ρ k p, ėk 2 + . (62)
Differentiating we obtain the first order approximation

dH(p) = 1≤k≤K ρ k p, ėk + ėk , thus dH p (du(p)) ≈ 1≤k≤K ρ k u(p) -u(p -h ėk ) h + ėk .
The slightly more general form (15) is handled similarly. Let us mention that our numerical codes also implement a second backtracking method, inspired by the diffuse numerical geodesics described in [START_REF] Benmansour | Derivatives with respect to metrics and applications: subgradient marching algorithm[END_REF], and which yields similar results. However we observed that many alternative backtracking methods did (often) fail, in particular with the Dubins model, due to the discontinuity of value function u.

In favorable cases, geodesic shooting is more precise than the PDE and backtracking approach. However it is also much less general and robust, for the following reasons. (I) Geodesic shooting cannot address the Reeds-Shepp forward and Dubins models, due to their non twice A comparison of geodesic shooting with the proposed approach is presented on Figure 4. A similar experiment was presented in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF] for the Reeds-Shepp reversible model, but not for the Euler-Mumford model.

Motion planning

Motion planning is a natural field of application for minimal path methods, see for instance [START_REF] Kimmel | Multivalued distance maps for motion planning on surfaces with moving obstacles[END_REF], and the different models considered in this paper may account for the maneuverability constraints of a number of vehicles. For instance, the reversible Reeds-Shepp model describes wheelchairs or wheelchair-like robots. The forward only variant is appropriate for vehicles of the same type, but which cannot see behind themselves. Minimal paths w.r.t. the Euler-Mumford elastica model have the advantage of being smooth, and their curvature penalization is physically meaningful. The hard upper bound constraint of the Dubins model on the curvature is relevant for the numerous vehicles subject to a minimum turning radius. See [START_REF] Mirebeau | Automatic differentiation of nonholonomic fast marching for computing most threatening trajectories under sensors surveillance[END_REF] for more applications to motion planning of the numerical methods presented in this paper, including two player games where an opponent choses the obstacles.

Our experiments, presented on Figures 7,8, 9 and 10, show that the numerical schemes introduced in this paper can be used to solve complex motion planning problems, on domains involving numerous walls, in CPU time below one second on a standard laptop8 . Within that time, the PDE (8) is numerically solved on the full domain, here discretized on a 90 2 × 60 grid, which yields a complete strategy to reach the given target. Larger grid scales yield more accurate paths, at the price of longer computation times (but still quasi-linear in the number of pixels). Memory usage is dominated by the storage of the value function, namely one floating point number per grid point. Once the PDE is solved, extracting a minimal path from an arbitrary point in the domain has an almost negligible cost.

The chosen domain for this experiment is [0, 1] 2 × S 1 , the relaxation parameter is ε = 0.1, and various values are considered for ξ, which we recall is homogeneous to a radius of curvature.

The different qualitative properties of the four models appear clearly, in particular the cusps of the Reeds-Shepp resersible model minimal paths, the in-place rotations of the Reeds-Shepp forward paths, the smoothness of the Euler-Mumford paths, and the lower-bounded radius of curvature of the Dubins paths.

On the rightmost sub-figures of Figure 10, illustrating the Dubins model, the author acknowledges that two paths fail, in some places, to obey the prescribed lower-bound on the radius of curvature. This is due to the approximate nature of the PDE discretization, and to the fact that no path obeying this lower bound seems to exist between these endpoints. These artifacts can easily be eliminated by post-processing.

Tubular structure segmentation

We present a numerical experiment, based on synthetic data, designed to illustrate the qualitative differences between the four implemented minimal path models in applications to tubular structure segmentation, see the Discussion in §1.3. Our experiment is illustrated on Figures 11, 12 and 13. We choose a local cost function α = α(x) taking small values iff x is close to a vessel, then we extract a minimal path between prescribed endpoints w.r.t. some curvature cost and parameter ξ, see the full description below. Some model/parameter combinations succeed, in the sense that the extracted minimal path goes along the desired vessel, and other fail, as discussed below. We would however like to emphasize that our results could be considerably improved by choosing a cost function α = α(x, θ) also depending on the orientation θ, see [BC10, SBD + 15]. Our "inefficient" choice of cost function α is intended to exacerbate the differences between the different models, by making the problem harder.

As noted by a reviewer, real medical data always contains noise, which is absent from our synthetic test case, see Figure 12. Due to noise, real applications do not directly define the speed function α as the input image, as we do here. Instead, the data-driven speed function is obtained as the result of a complex pre-processing step, aimed at reducing noise and locally detecting the image features, see the above references on e.g. tubular structure segmentation. We choose to avoid this complex machinery (thus do not consider noise) so as to keep the experiment as simple, reproducible, and interpretable as possible. Indeed, our objective is not to address a difficult test case, but to emphasize the qualitative differences between the four implemented path models.

The orientation lifted image domain Ω = [0, 1] 2 × S 1 , is discretized on a 73 2 × 60 grid, so that computation times remain within ≈ 1s. Note that the O(N ln N ) time complexity and O(N ) space complexity of the proposed numerical method makes it suitable for larger problems as well. The cost function α = α(x), shown in grayscale on Figures 11, 12 and 13, is set to α ≈ 1/4 in the neighborhood of the S shaped tubular structure, α = 1/3 close to the straight one, and α = 1 elsewhere. The HJB PDE (8) is numerically solved, with u = 0 at two seed points shown in blue on Figure 11, and with outflow boundary conditions on ∂Ω. Note that by point we do mean an element of M = R 2 × S 1 , displayed as an arrow on the figures. Minimal paths from two points, shown in red, are backtracked. The experiment is called a success if these paths do follow the S tubular structure, and end in the leftmost blue point.

The Reeds-Shepp reversible model fails this test, because the minimal paths shift into reverse gear, and thus end at the incorrect blue point. The Reeds-Shepp forward model performs slightly better, ending at the correct endpoint, and correctly extracting the bottom part of the S structure for the larger value of ξ. However the top part of the S structure is not correctly extracted, because the minimal path intended to do so takes a shortcut through the second straight structure, due to its ability to perform in-place rotations. (Again, let us emphasize that the Reeds-Shepp models can be perfectly suitable for tubular structure segmentation, when contrary to this experiment an orientation dependent cost function α = α(x, θ) is provided, see [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 model is able to extract the S shaped tubular structure for a large range of parameters ξ, from 0.55 to 1.25. Excessively small values of ξ allow tight turns and therefore shortcuts through the intersecting structure (recall that ξ should be interpreted as a radius of curvature). Excessively large values of ξ make it too costly to follow the S structure, hence the path takes wider turns in the background of the image. We regard the Euler-Mumford model as the best choice among those considered here for image segmentation tasks, in view of our numerical experiments and of the literature [START_REF] Mumford | Elastica and computer vision[END_REF][START_REF] Chen | Global Minimum for a Finsler Elastica Minimal Path Approach[END_REF].

The Dubins model is also able to extract the S shaped tubular structure, but for a narrower range of parameters ξ than in the Euler-Mumford case, from 0.25 to 0.38. Excessively small or excessively large choices of ξ lead to symptoms, respectively shortcuts and excursions in the image background, similar to those observed in the Euler-Mumford case. We expect the Dubins model to be less efficient than the Euler-Mumford model in practical image segmentation tasks, since it requires a finer tuning of the parameter ξ, matching the curvature of the extracted structures.

Conclusion

In this paper, we introduced numerical PDE methods for solving a family of non-holonomic optimal control problems, associated to the Reeds-Shepp (reversible or forward only), Euler-Mumford and Dubins models. The design and analysis of these methods uses tools from different branches of mathematics, including (I) an original reformulation of the Euler-Mumford Hamiltonian, (II) a convergence analysis in the setting of discontinuous viscosity solutions to HJB PDEs, and (III) a finite differences scheme based on lattice geometry. The discretized PDEs are solved in a single pass via the dynamic programming principle, which guarantees fast computation times. Synthetic test cases illustrate the potential of our methods in motion planning and image segmentation, and can be reproduced using free and open source codes.

Future works include (i) computing three dimensional9 Euler-Mumford and Dubins minimal paths, by solving PDEs on R 3 × S 2 , and addressing other instances of non-holonomic optimal control problems, (ii) implementing GPU accelerations of the solver, in the spirit of [WDB + 08], (iii) designing mesh-based, instead of cartesian grid based, discretizations of the HJB PDEs considered in this paper, and (iv) developing practical applications of our globally optimal, curvature penalized minimal paths.

A Local validation of the numerical scheme, via the control sets

We present a local validation of our PDE discretization procedure, by comparing the model control sets (4), see also Figure 2, with some numerical counterparts.

Consider a compact and convex set B ⊆ E = R 3 , containing the origin, and the corresponding Hamiltonian H. Contrary to the rest of this paper, the possible dependency of B and H on some underlying base point p ∈ M is not considered, since the discussion is purely local. Let also H be an approximation of H, for instance of the form (15). Consider the set

B := { ṗ ∈ E; ∀p ∈ E * , H(p) ≤ 1/2 ⇒ p, ṗ ≤ 1}, (63)
and note that B = B if H = H. We regard the closeness of B and B, inspected visually, as a good witness of the closeness of H and H. On Figure 5, left, is illustrated the case where B is an ellipsoid, with principal axes of length (1, 0.1, 0.1). In that case H is a quadratic function, and the discrete representation using Corollary 4.13 is exact. One therefore has H = H, thus B = B as can be observed. This particular case is at the foundation of [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF].

On figure 5, right, is illustrated the case where B = [0, ṅ] is a segment, where ṅ ∈ E is a unit vector, and therefore H(p) = 1 2 p, ṅ 2 + . The discretization is performed using Proposition 2.2 with ε = 0.1, via the basis reduction techniques presented in §4. As can be observed, the vectors ėk , 1 ≤ k ≤ K, are almost aligned with ṅ, and the set B is close to the segment B in the Hausdorff distance, although slightly fatter.

Figure 6 is devoted to the models of interest in this paper. The parameters are ξ = 0.2 (appearing in the curvature cost (2)), θ = π/3 (the current orientation), and ε = 0.1 (tolerance in Proposition 2.2). The offsets ėk , 1 ≤ k ≤ K, are illustrated on Figure 3 page 9. Comparing with Figure 2, we can confirm that the sets B are close to the corresponding ellipse 13 , half ellipse, non-centered ellipse, and triangle B respectively. However the true control sets B are flat, with Haussdorff dimension 2, whereas their counterparts B are slightly fatter.

B Convexity of the metric

We prove in this appendix that the metrics F RS+ p , F EM p , F D p : E → [0, ∞] are convex, for any fixed p ∈ Ω, due their construction (3) and to the following two results. 13 The choice ξ = 1 in Figure 2, versus ξ = 0.2 in Figure 6, yields a round disk instead of an ellipse Proof. Lower semi-continuity, homogeneity and positivity are obvious. In addition for any (n, t), (n , t ) ∈]0, ∞[×R one obtains

f (n + n , t + t ) n + n = C t n n n + n + t n n n + n ≤ n n + n C( t n ) + n n + n C( t n ) = f (n, t) + f (n , t ) n + n .
Note that 1-positive homogeneity and the triangular inequality together imply convexity. We recall the notation E := R 2 × R, used in the next result.

Corollary B.2. Let C : R → [1, ∞] be convex, lower semi-continuous, and such that l(ε) := lim C(εt)/t as t → ∞ exists and belongs to ]0, ∞], for each ε ∈ {-1, 1}. Let F : E → [0, ∞] be defined as F ( ẋ, θ) := ẋ C( θ/ ẋ ) for each ( ẋ, θ) ∈ E such that ẋ = ẋ n and ẋ > 0, where n ∈ S d-1 is a fixed direction. Let also F (0, 0) = 0, F (0, θ) = | θ|l( θ/| θ|) if θ = 0, and F ( ẋ, θ) := +∞ otherwise. Then F is lower semi-continuous and obeys

• (1-positive homogeneity) F (λ ṗ) = λF ( ṗ), for all λ > 0 and all ṗ ∈ E.

• (Separation) F ( ṗ) = 0 iff ṗ = 0, for all ṗ ∈ E.

• (Triangular inequality) F ( ṗ + q) ≤ F ( ṗ) + F ( q), for all ṗ, q ∈ E.

Proof. By the previous lemma, F obeys the announced properties on the convex set (]0, ∞[n)×R. These properties are also satisfied on the closure ([0, ∞[n) × R since, clearly, F is extended to it by its lower continuous envelope, and since the limits l(1) and l(-1) are positive (for separation). Finally the announced properties hold for the trivial extension of F to R 2 × R by +∞ since the subset ([0, ∞[n) × R is closed and convex.

Figure 1 :

 1 Figure 1: Globally minimal paths for the Reeds-Shepp reversible model (ξ = 0.3), Reeds-Shepp forward model (ξ = 0.3), Euler-Mumford elastica model (ξ = 0.2), and Dubins model (ξ = 0.2), with uniform cost α ≡ 1. Seed point (1/2, 1/2, 0) ∈ R 2 × S 1 , tip points (k/3, l/3, π/4), k, l ∈ {0, • • • , 3}. (Results obtained using our numerical schemes, see §2, with relaxation parameter ε = 0.1, and angular resolution 2π/60.)

FFigure 2 :

 2 Figure 2: Control sets of the Reeds-Shepp reversible, Reeds-Shepp forward, Euler-Mumford and Dubins models, see (4). All have empty interior, reflecting the non-holonomy of the models.

1h

  (u(p) -u(p -h ė)) + , where h > 0 is the discretization grid scale.

Figure 3 :

 3 Figure 3: Discretization stencils used for the Reeds-Shepp reversible, Reeds-Shepp forward, Euler-Mumford, and Dubins models. Note the sparseness and anisotropy of the stencils. Model parameters: θ = π/3, ξ = 0.2. Discretization parameters: ε = 0.1, and for Euler-Mumford K = 5.

  hn , where ε n → 0 and h n /ε n → 0 as n → ∞. Define for all p ∈ Ω u(p) := lim r→0 lim inf n→∞ inf M h ∩B(p,r) U n , u(p) := lim r→0 lim sup n→∞ sup M h ∩B(p,r)

  Proposition 3.3. Under the assumptions of Theorem 3.1, and in addition int(Ω) = Ω. The value functions u, û : Ω → R are equal in the following cases: • (Reeds-Shepp reversible model) u = û on Ω. • (Reeds-Shepp forward model) u = û on Ω, if this domain has the form Ω = Ω 0 × S 1 . • (Euler-Mumford model) u = û on Ω. • (Dubins model) u = û on a dense subset of Ω. Furthermore, in each case, u and û are continuous at each point p ∈ Ω where u(p) = û(p).

  right-hand sides are continuous and non-vanishing functions of ϕ ∈ [-π, π]. More precisely, one easily finds the following lower bounds: for the two Reeds-Shepp models and the Dubins model c 0 = 1, attained for ϕ = 0; for the Euler-Mumford elastica model c 0 = π/4, attained for ϕ = π/2.

  for any p ∈ Ω. Indeed, consider a Lipschitz path γ : [0, 1] → Ω from ∂Ω to p ∈ Ω and such that u(p) = length F (γ) < ∞. Introduce a time reparametrization η : [0, T ] → Ω of γ, from ∂Ω to p, at unit speed w.r.t. the reversed metric in the sense that F η(t) (-η(t)) = 1 for all t ∈ [0, T ]. Then clearly T = length F (γ) = u(p), and one can uniquely define controls a : [0, T ] → A by f (η(t), a(t)) = -η(t) for all t ∈ [0, T ]. Introducing an additional time-reversal reparametrization, t ∈ [0, T ] → T -t, required since (44) considers contrary to us paths from p to ∂Ω, we obtain v(p) ≤

  and using Corollary 3.20 one obtains the announced result. Points (II) and (III) of Proposition 3.3 are established in the next proposition and corollary.

  Finally û(p) = u(p) since by construction û(p) ≥ u(p). Corollary 3.23. Under the assumptions of Proposition 3.3. The Reeds-Shepp forward model satisfies u = û on Ω if this domain has the shape Ω = U × S 1 . The Euler-Mumford model satisfies u = û on Ω.

  dual space and dual lattice are denoted E * d and L * d . Thanks to the Euclidean structure, there is a canonical identification E d ∼ = E * d and L d ∼ = L * d , but the distinction is kept for clarity. Let S(E d ) ⊆ L(E d , E * d ) denote the set of symmetric linear maps, and let S ++ (E d ) ⊆ S + (E d ) ⊆ S(E d ) be the subsets of positive definite and semi-definite ones. To each M ∈ S ++ (E d ) is associated a scalar product •, • M and a norm • M , defined for all ė, ḟ ∈ E d by ė, ḟ M := M ė, ḟ , ė M := ė, ė M . Let us point out that it is equivalent, up to a linear change of coordinates, to study the geometry of an arbitrary lattice of E d w.r.t. the canonical Euclidean norm, or to study the geometry of the canonical lattice L d w.r.t. the norm • M induced by an arbitrary M ∈ S ++ (E d ). Voronoi's first reduction [Sch09] consists of a convex set P ⊆ S ++ (E * d ), and for each M ∈ S ++ (E d ) of a linear programming problem L(M ):

  . The perfect forms {D b ; b superbase of L d } are such a class, and in fact are the only one 4 in dimension d ≤ 3. Lemma 4.5. Let b = ( ė0 , • • • , ėd ) be a superbase of L d , and let D b Then vij D b = 1 for all 0 ≤ i < j ≤ d and v D b ≥ 1 for any v ∈ L * d \ {0}. In particular D b is a vertex of the polytope P ⊆ S(E d ), at the intersection of the facets defined by {D ∈ P; D, vij ⊗ vij = 1}, i, j ∈ 0, d . Proof. Let δ I i = 1 if i ∈ I and δ I i = 0 otherwise, for any i ∈ 0, d and I ⊆ 0, d . Clearly (δ I i ) 2 = δ I i , and for any pairwise

  D b = 1/2. Likewise vij , vkl D b = 0, and vij 2 D b = vij , vij D b = 1 as announced. Since (v 12 , • • • , vd0 ) is a basis of L * d , see Proposition 4.3, we obtain by bilinearity that v 2 D b ∈ Z for all v ∈ L * d . Since D b is positive definite, one has v D b > 0 for all v ∈ E * d \{0}, hence v D b ≥ 1 for all v ∈ L * d \ {0} as announced. By definition D b ∈ P, and D b obeys the d(d + 1)/2 = dim S(E d ) linear equalities D b , vij ⊗ vij = 1, for all 0 ≤ i < j ≤ d, which are linearly independent by Lemma 4.4. Thus D b is a vertex of P. The two previous lemmas, combined, yield a complete solution to the optimization problem (49), when an M -obtuse superbase is known. Proposition 4.6. Let M ∈ S ++ (E d ), and let us assume that there exists an M -obtuse superbase b of L d . Then D b is a minimizer of L(M ), and applying Lemma 4.4 to the basis b yields the explicit value L(M ) = i<j λ ij , as well as an optimality certificate: the Karush-Kuhn-Tucker conditions are (54, left). Proof. Clearly D b ∈ P, and for any D ∈ P one has M, D -D b

  Proposition 4.7 (Selling's algorithm). Let M ∈ S ++ (E d ), where d ∈ {2, 3}, and let b = ( ė0 , • • • , ėd ) be a superbase of L d . Define a second superbase b of L d by Case d = 2 : b := (-ė0 , ė1 , ė0 -ė1 ), Case d = 3 : b := (-ė0 , ė1 , ė0 + ė2 , ė0 + ė3 ). Then Tr(M D b ) -Tr(M D b ) = 2 2-d ė0 , ė1 M . Selling's algorithm consists in iteratively, and until b is an M -obtuse superbase: (a) reordering the superbase b so that ė0 , ė1 M > 0, and (b) applying the transformation b ← b . This algorithm terminates, and in particular there exists an M -obtuse superbase of L d . Proof. Applying Definition 4.1 we find that b is indeed a superbase of L d . The expression of Tr(M D b ) -Tr(M D b ) follows from a direct computation. Denoting by b n be the superbase obtained at the n-th step of Selling's algorithm, one observes that Tr(M D bn ) strictly decreases as n increases by construction. In view of 2 d Tr(M D b ) = there are only finitely many superbases b = ( ė0 , • • • , ėd ) ∈ L d d such that Tr(M D b ) is below any given bound, hence Selling's algorithm must terminate, and the result follows.

  Theorem 4.11. Let ( ė0 , • • • , ėd ) be an M -obtuse superbase of Z d , where M ∈ S ++ (E d ) and d ∈ {2, 3}. Then vij ≤ 2C d Cond(M ) for all i = j. Proof. Case d = 2. By Propositions 4.8 and 4.9, one has vij ≤ C d Cond(D) = C d Cond(M ) as announced (the factor 2 is here useless). Case d = 3. By Propositions 4.8 and 4.10, assuming w.l.o.g. that b0 is D-obtuse, the co-vectors v01 , v12 , v23 , v30 have their norm bounded by C d Cond(D) = C d Cond(M ). The remaining two co-vectors v02 and v13 can be expressed as: v13 = ė2 × ė0 = ė2 × (-ė1 -ė2 -ė3 ) = -ė2 × ė1 -ė2 × ė3 = v30 + v01 , and likewise v02 = v01 + v12 . Thus v13 ≤ v30 + v01 ≤ 2C d Cond(M ), and likewise v02 ≤ 2C d Cond(M ). The announced result follows.

Furthermore, v ≤

 v 2C d Cond(M ) whenever ρ v(M ) > 0, and v∈L * d ρ ε v(n) v 2 = Tr(M ).Proof. By Proposition 4.7 there exists an M -obtuse superbase, by Lemma 4.4 it yields a decomposition of the announced form (58), and by Theorem 4.11 the norms of the support co-vectors are bounded as announced. The last identity follows by taking the trace of (57).

Corollary 4. 13 .

 13 Let d ∈ {2, 3}, let 0 < ε ≤ 1, and let n ∈ E * d such that n = 1. Then there exists non-negative weights ρ ε v(n), supported on at most d(d + 1)/2 elements v ∈ L * d , such that

  Consider a model with a smooth Hamiltonian H, such as the reversible Reeds-Shepp model or the Euler-Mumford model. Paths of minimal energy are known to obey 7 the Hamilton equations of geodesics, which read

Figure 4 :

 4 Figure 4: Comparison of Minimal geodesics obtained by PDE resolution and backtracing (black), or geodesic shooting (blue), for Euler-Mumford (left) and Reeds-Shepp reversible models in 2D (center) or 3D (right).

Figure 5 :

 5 Figure 5: Illustration of §A. Left: Discretization stencil for a quadratic anisotropic Hamiltonian. The Hamiltonian and control set representation are exact. Right: Discretization stencil of Proposition 2.2, for a Hamiltonian H(p) = 1 2 p, ṅ 2 + . The approximate control set (63) is close to the segment [0, ṅ], but slightly fatter.

Figure 6 :

 6 Figure 6: Approximate control set (63) for the Reeds-Shepp reversible, Reeds-Shepp forward, Euler-Mumford and Dubins models. The shapes of the original control sets, see Figure 2 are recognizable, elongated due to the model parameter choice ξ = 0.2, and slightly fatter due to the effect of discretization, with relaxation parameter ε = 0.1. Orientation θ = π/3. Seed Figure 3 for the corresponding stencils.

Figure 11 :

 11 Figure 11: Tubular segmentation test with the Reeds-Shepp models, reversible (left), and forward (right). First image ξ = 0.2, second ξ = 0.8. This segmentation test is mostly "failed", on purpose, for different reasons, see §5.3. CPU time ≈ 0.3s.

Figure 12 :

 12 Figure 12: Tubular segmentation using Euler-Mumford elasticas. Curvature penalization is left: insufficient (ξ = 0.2), middle: adequate (ξ = 0.6), right: exaggerate (ξ = 1.5). CPU time ≈ 1.2s.

Figure 13 :

 13 Figure 13: Tubular segmentation using the Dubins model. Curvature penalization is left: insufficient (ξ = 0.2), middle: adequate (ξ = 0.35), right: exaggerate (ξ = 0.5). CPU time ≈ 0.5s.

  ė and reasoning similarly, we obtain the upper bound (36, right) multiplied by (C WS h/ε) 2 .The announced result is equivalent to

The Reeds-Shepp models are in this paper mostly discussed for comparison with the Euler-Mumford and Dubins models, since our numerical results in the Reeds-Shepp case are actually quite similar to[START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF], despite the distinct discretization.

This specific discretization requires the metric tensors condition number to remain below ( √

+ 1)/2, see[START_REF] Mirebeau | Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction[END_REF] for an explanation of this bound and an unconditional method.

Using a semi-lagrangian discretization for three dimensional the Reeds-Shepp reversible model, and the present discretization for the three dimensional forward model.

Indeed, as shown below in Proposition 4.6, the minimum (49) of the linear program L(M ) on K is always attained at a matrix D b attached to a superbase b.

The scalar product vk , n makes sense thanks to the Euclidean structure on E d and E * d .

github.com/Mirebeau/HamiltonFastMarching

Except perhaps abnormal geodesics, which we do not discuss here, see[START_REF] Montgomery | A Tour of Subriemannian Geometries, Their Geodesics and Applications[END_REF].

Laptop processor: 2.7GHz Intel R Core i7 using a single core

The three dimensional Reeds-Shepp model is already adressed here and in[START_REF] Duits | Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis[END_REF].

Post-Doctoral researcher at University Paris-Dauphine.

PhD student under the direction of R. Duits at TU/e University, Eindhoven.

Post-Doctoral researcher at TU/e University, Eindhoven.

Acknowledgement. The author thanks Da-Chen 10 , Jorg Portegies 11 and Erik Bekkers 12 , for careful testing, bug-fixing, and feed-back on the numerical codes.

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

The Karush-Kuhn-Tucker optimality conditions for this problem read as follows: there exists non-negative weights λ k ≥ 0 and integral co-vectors vk ∈ L * d , where 1

This formula is reminiscent of the eigenvector-eigenvalue decomposition, but the co-vectors vk have integer coordinates, and the number of terms is larger: d(d + 1)/2 instead of d. It is used to design monotone finite differences PDE schemes in [START_REF] Bonnans | A fast algorithm for the two dimensional HJB equation of stochastic control[END_REF][START_REF] Fehrenbach | Sparse non-negative stencils for anisotropic diffusion[END_REF][START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using Voronoi's first reduction of quadratic forms[END_REF], implicitly in the first two references. In §4.1 we describe Selling's algorithm which is a constructive, simple and efficient method for solving (50) in dimension d ≤ 3, used as-is in our numerical experiments. We estimate in §4.2 the largest norm of the co-vectors appearing in (51). We finally conclude in §4.3 the proof of Proposition 2.2.

Obtuse superbases and Selling's algorithm

We introduce in this section the concept of obtuse superbase of lattice, a preferred coordinate system which provides, in particular, a complete solution to Voronoi's first reduction (49), see Proposition 4.6. An obtuse superbase exists for all lattices of dimension two and three, and can be constructed using Selling's algorithm, see Proposition 4.7. The results of this subsection are mostly reformulations of [CS92, BK10], but they are prerequisites for the original results of the next subsections. Proofs are provided for completeness.

We attach to each superbase a family of d(d + 1) co-vectors. Definition 4.2. Let ( ė0 , • • • , ėd ) be a superbase of of L d . For all distinct i, j ∈ 0, d , we define a co-vector vij ∈ L * d by the linear equalities

Note that vij = -v ji , and that the definition (52) of the d-dimensional co-vector vij using d + 1 linear constraints makes sense thanks to the (single) redundancy relation

In dimension d ∈ {2, 3} the construction of Definition 4.2 has a simple geometrical interpretation

where (i, j, k) is a permutation of {0, 1, 2}, (resp. (i, j, k, l) of {0, 1, 2, 3}), and the ± sign refers to its signature.

for this is that most studies on lattice geometry consider a single norm, and not the interaction of two norms as here, namely the Euclidean norm • and an anisotropic norm • M . Our first step is to upper bound the norm of the elements of an M -obtuse superbase. Throughout this subsection, we denote by vij ∈ L * d the co-vectors associated by Definition 4.2 to a superbase ( ė0 , • • • , ėd ). Observe also that Id, D b ∈ P, and that D b is optimal for (50). Thus, as announced

In order to proceed, we recall some identities relating the scalar products associated to M and D := M -1 , where M ∈ S ++ (E d ). For any ė0 , ė1 , ė2

where × denotes the cross product of three dimensional vectors. These identities are easily checked when M = Id, and otherwise follow from a change of variables by M Recalling that vij = -v ji for any i = j one observes the following: if α 1 ≥ 0 ≥ α 2 then b0 is D-obtuse, if α 2 ≥ 0 ≥ α 0 then b1 is D-obtuse, and if α 0 ≥ 0 ≥ α 1 then b2 is D-obtuse. On the other hand, expressing (α i ) 2 i=0 in terms of cross products yields

and inserting the identity ė0 := -( ė1 + ė2 + ė3 ) in the above expressions one obtains that α 0 + α 1 + α 2 = 0. Therefore these three scalars cannot have the same sign strictly. Thus in the sequence α 0 , α 1 , α 2 , α 0 there is at least one non-negative scalar α i followed with a non-positive scalar α i+1 , where indices are understood modulo 3. By the previous argument, bi-1 is D-obtuse, and the announced result follows.