Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance

07/11/2017
Publication GSI2017
OAI : oai:www.see.asso.fr:17410:22600
DOI :
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit
 

Résumé

We consider a two player game, where a rst player has to install a surveillance system within an admissible region. The second player needs to enter the monitored area, visit a target region, and then leave the area, while minimizing his overall probability of detection. Both players know the target region, and the second player knows the surveillance installation details. Optimal trajectories for the second player are computed using a recently developed variant of the fast marching algorithm, which takes into account curvature constraints modeling the second player vehicle maneuverability. The surveillance system optimization leverages a reverse-mode semi-automatic di erentiation procedure, estimating the gradient of the value function related to the sensor location in time O(N lnN).

Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance

Collection

application/pdf Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance Jean-Marie Mirebeau, Johann Dreo
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance
application/pdf Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance (slides)

Média

Voir la vidéo

Métriques

0
0
1.28 Mo
 application/pdf
bitcache://d374bc58c01b7a5a80a3e4f1a85398b66df5d19f

Licence

Creative Commons Aucune (Tous droits réservés)

Sponsors

Sponsors Platine

alanturinginstitutelogo.png
logothales.jpg

Sponsors Bronze

logo_enac-bleuok.jpg
imag150x185_couleur_rvb.jpg

Sponsors scientifique

logo_smf_cmjn.gif

Sponsors

smai.png
logo_gdr-mia.png
gdr_geosto_logo.png
gdr-isis.png
logo-minesparistech.jpg
logo_x.jpeg
springer-logo.png
logo-psl.png

Organisateurs

logo_see.gif
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xmlns="http://datacite.org/schema/kernel-4"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/17410/22600</identifier><creators><creator><creatorName>Jean-Marie Mirebeau</creatorName></creator><creator><creatorName>Johann Dreo</creatorName></creator></creators><titles>
            <title>Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2018</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><subjects><subject>Anisotropic fast-marching</subject><subject>motion planning</subject><subject>sensors placement</subject><subject>game theory</subject><subject>optimization</subject></subjects><dates>
	    <date dateType="Created">Thu 8 Mar 2018</date>
	    <date dateType="Updated">Thu 8 Mar 2018</date>
            <date dateType="Submitted">Mon 15 Oct 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">d374bc58c01b7a5a80a3e4f1a85398b66df5d19f</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>37325</version>
        <descriptions>
            <description descriptionType="Abstract">We consider a two player game, where a rst player has to install a surveillance system within an admissible region. The second player needs to enter the monitored area, visit a target region, and then leave the area, while minimizing his overall probability of detection. Both players know the target region, and the second player knows the surveillance installation details. Optimal trajectories for the second player are computed using a recently developed variant of the fast marching algorithm, which takes into account curvature constraints modeling the second player vehicle maneuverability. The surveillance system optimization leverages a reverse-mode semi-automatic di erentiation procedure, estimating the gradient of the value function related to the sensor location in time O(N lnN).
</description>
        </descriptions>
    </resource>
.

Automatic differentiation of non-holonomic fast marching for computing most threatening trajectories under sensors surveillance Jean-Marie Mirebeau1 and Johann Dreo2 1 University Paris-Sud, CNRS, University Paris-Saclay, jean-marie.mirebeau@math.u-psud.fr 2 THALES Research & Technology, johann.dreo@thalesgroup.com Abstract. We consider a two player game, where a first player has to install a surveillance system within an admissible region. The second player needs to enter the monitored area, visit a target region, and then leave the area, while minimizing his overall probability of detection. Both players know the target region, and the second player knows the surveil- lance installation details. Optimal trajectories for the second player are computed using a recently developed variant of the fast marching algo- rithm, which takes into account curvature constraints modeling the sec- ond player vehicle maneuverability. The surveillance system optimization leverages a reverse-mode semi-automatic differentiation procedure, esti- mating the gradient of the value function related to the sensor location in time O(N ln N). Keywords: Anisotropic fast-marching, motion planning, sensors place- ment, game theory, optimization 1 Introduction This paper presents a proof of concept numerical implementation of a motion planning algorithm related to a two player game. A first player selects, within an admissible class Ξ, an integral cost function on paths, which takes into account their position, orientation, and possibly curvature. The second player selects a path, within an admissible class Γ, with prescribed endpoints and an interme- diate keypoint. The players objective is respectively to maximize and minimize the path cost C(Ξ, Γ) := sup ξ∈Ξ inf γ∈Γ C(ξ, γ), where C(ξ, γ) := Z T (γ) 0 Cξ(γ(t), γ0 (t), γ00 (t)) dt, (1) where the path γ is parametrized at unit Euclidean speed, and the final time T(γ) is free. From a game theoretic point of view, this is a non-cooperative zero- sum game, where player Ξ has no information and player Γ has full information over the opponent’s strategy. The game (1) typically models a surveillance problem [20], and exp(−C(Ξ, Γ)) is the probability for player Γ to visit a prescribed keypoint without being de- tected by player Ξ. For instance player Ξ is responsible for the installation of radar [1] or sonar detection systems [20], and would like to prevent vehicles sent by player Γ from spying on some objectives without being detected. The dependence of the cost Cξ w.r.t. the path tangent γ0 (t) models the varia- tion of a measure of how detectable the target is (radar cross section, directivity index, etc.) w.r.t. the relative positions and orientations of the target and sensor. The dependence of Cξ on the path curvature γ00 (t) models the airplane maneu- verability constraints, such as the need to slow down in tight turns [9], or even a hard bound on the path curvature [8]. Strode [20] has shown the interplay of motion planning and game theory in a similar setting, on a multistatic sonar network use case, but using isotropic graph-based path planning. The same year, Barbaresco [2] used fast-marching for computing threatening paths toward a single radar, but without taking into account curvature constraints and without considering a game setting. The main contributions of this paper are as follows: 1. Anisotropy and curvature penalization: Strategy optimization for player Γ is an optimal motion planning problem, with a known cost function. This is addressed by numerically solving a generalized eikonal PDE posed on a two or three dimensional domain, and which is strongly anisotropic in the presence of a curvature penalty and a detection measurement that depends on orientation. A Fast-Marching algorithm, relying on recent adaptive sten- cils constructions, based on tools from lattice geometry, is used for that purpose [9, 14, 15]. In contrast, the classical fast marching method [17] used in [5] is limited to cost functions Cξ(γ(t)) independent of the path orienta- tion γ0 (t) and curvature γ00 (t). 2. Gradient computation for sensors placement: Strategy optimization for player Ξ is typically a non-convex problem, to which various strategies can be ap- plied, yet gradient information w.r.t. the variable ξ ∈ Ξ is usually of help. For that purpose, we implement efficient differentiation algorithms, forward and reverse, for estimating the gradient of the value function of player Ξ ∇ξC(ξ, Γ), where C(ξ, Γ) := inf γ∈Γ C(ξ, γ). (2) Reverse mode differentiation reduced the computation cost of ∇ξC(ξ, Γ) from O(N2 ), as used in [5], to O(N ln N), where N denotes the number of discretization points of the domain. As a result, we can reproduce exam- ples from [5] with computation times reduced by several orders of magnitude, and address complex three dimensional problems. Due to space constraints, this paper is focused on problem modeling and numerical experiments, rather than on mathematical aspects of wellposedness and convergence analysis. Free and open source codes for reproducing (part of) the presented numerical experiments are available on the first author’s webpage3 . 3 github.com/Mirebeau/HamiltonianFastMarching 2 Mathematical background of trajectory optimization We describe in this section the PDE formalism, based on generalized eikonal equations, used to compute the value function minγ∈Γ C(ξ, γ) of the second player, where ξ is known and fixed. Their discretization is discussed in §3. We dis- tinguish two cases, depending on whether the path local cost function Cξ(x, ẋ, ẍ) appearing in (1) depends on the last entry ẍ, i.e. on path curvature. 2.1 Curvature independent cost Let Ω ⊆ E := R2 be a bounded domain, and let the source set Υ and target set Θ be disjoint subsets of Ω. For each x ∈ Ω, let Γx denote the set of all paths γ ∈ C1 ([0, T], Ω), where T = T(γ) is free, such that γ(0) ∈ Υ, γ(T) = x and ∀t ∈ [0, T], kγ0 (t)k = 1. The problem description states that the first player needs to go from Υ to Θ and back, hence its set of strategies is Γ = S x∈Θ Γ+ x × Γ− x , where Γ+ x = Γ− x := Γx, and C(ξ, Γ) = inf x∈Θ u+ ξ (x) + u− ξ (x), where u± ξ (x) := inf γ∈Γ ± x C± (ξ, γ). (3) Here and below, the symbol “±” must be successively replaced with “+” and then “−”. We denoted by C± the path cost defined in terms of the local cost Cξ(x, ±ẋ). In practice though, we only consider symmetric local costs, obeying Cξ(x, ẋ) = Cξ(x, −ẋ), hence the forward and return paths are identical and we denote uξ := u+ ξ = u− ξ . Define the 1-homogenous metric Fξ : Ω × E → [0, ∞], the Lagrangian Lξ and the Hamiltonian Hξ by Fξ(x, ẋ) := kẋkCξ(x, ẋ/kẋk), Lξ := 1 2 F2 ξ , Hξ(x, x̂) := sup ẋ∈E hx̂, ẋi − Lξ(x, ẋ). Here and below, symbols denoting tangent vectors are distinguished with a “dot”, e.g. ẋ, and co-vectors with a “hat”, e.g. x̂. Under mild assumptions [3], the function uξ : Ω → R is the unique viscosity solution to a generalized eikonal equation ∀x ∈ Ω \ Υ, Hξ(x, ∇xuξ(x)) = 1/2, ∀x ∈ Υ, uξ(x) = 0, with outflow boundary conditions on ∂Ω. The discretization of this PDE is discussed in §3. We limit in practice our attention to Isotropic costs Cξ(x), and Riemannian costs Cξ(x, ẋ) = p hẋ, Mξ(x)ẋi where Mξ(x) is symmetric positive definite, for which efficient numerical strategies have been developed [12, 14]. 2.2 Curvature dependent cost Let Ω ⊆ R2 × S1 be a bounded domain, within the three dimensional space of all positions and orientations. As before, let Υ, Θ ⊆ Ω. For all x ∈ Ω let Γ± x be the collection of all γ ∈ C2 ([0, T], Ω), such that η := (γ, ±γ0 ) satisfies η(0) ∈ Υ, η(T) = x and ∀t ∈ [0, T], kγ0 (t)k = 1. Since the first player needs to go from Υ to Θ and back, its set of strategies is Γ = S x∈Θ Γ+ x × Γ− x . Equation (3) holds, where C± denotes the path cost defined in terms of the local cost Cξ(p, ±ṗ, p̈). Consider the 1-homogeneous metric F± ξ : TΩ → [0, ∞], defined on the tan- gent bundle to Ω ⊆ R2 × S1 by F± ξ ((p, n), (ṗ, ṅ)) := ( +∞ if ṗ 6= kṗkn, kṗkCξ(p, n, ±ṅ/kṗk) else, where p ∈ R2 , n ∈ S1 is a unit vector, and the tangent vector satisfies ṗ ∈ R2 , ṅ ⊥ n. This choice is motivated by the fact that R T 0 F± ξ (η(t), η0 (t))dt is finite iff η : [0, T] → Ω is of the form (γ, ±γ0 ), and then it equals R T 0 Cξ(γ(t), ±γ0 (t), γ00 (t))dt. Introducing the Lagrangian L± ξ = 1 2 (F± ξ )2 on TΩ, and its Legendre-Fenchel dual the Hamiltonian H± ξ , one can again under mild assumptions characterize u± ξ as the unique viscosity solution to the generalized eikonal PDE H± ξ (x, ∇u± ξ (x)) = 1/2 with appropriate boundary conditions [3]. In practice, we choose cost func- tions of the form Cξ(p, ṗ, p̈) = C◦ ξ (p, ṗ)C?(|p̈|), where C? is the Reeds-Shepp car or Dubins car [8] curvature penalty, with respective labels ? = RS and D, namely CRS(κ) := p 1 + ρ2κ2, CD(κ) := ( 1 if |ρκ| ≤ 1, +∞ otherwise, where ρ > 0 is a parameter which has the dimension of a curvature radius. The Dubins car can only follow paths which curvature radius is ≤ ρ, whereas the Reeds-Shepp car (in the sense of [9] and without reverse gear), can ro- tate into place if needed. The Hamiltonian then has the explicit expression H((p, n), (p̂, n̂)) = 1 2 C0 ξ (p, n)−2 H∗(n, (p̂, n̂)) where HRS = 1 2 (hp̂, ni2 + + kn̂/ρk2 ) and HD = 1 2 max{0, hp̂, ni + kn̂/ρk}2 . 3 Discretization of generalized eikonal equations We construct a discrete domain X by intersecting the computational domain with an orthogonal grid of scale h > 0: X = Ω ∩ (hZ)d , where d = 2 for the curvature independent models, and d = 3 for the other models which are posed on R2 ×(R/2πZ) —using the angular parametrization S1 ∼ = R/2πZ (in the latter periodic case, 2π/h must be an integer). We design weights cξ(x, y), x, y ∈ X such that for any tangent vector ẋ at x ∈ Ω one has Hξ(x, ẋ) ≈ h−2 X y∈X c2 ξ(x, y)hx − y, ẋi2 +, (4) where a+ := max{0, a} (expression (4) is typical, although some models require a slight generalization). The weights cξ(x, y) are non-zero for only few (x, y) ∈ X at distance kx − yk = O(h). Their construction exploits the additive structure of the discretization grid X and relies on techniques from lattice geometry [19], see [9, 14, 15] for details. The generalized eikonal PDE Hξ(x, ∇xuξ(x)) = 1/2, which solution uξ(x) should be regarded as a distance map, is discretized as X y∈X c2 ξ(x, y)(Uξ(x) − Uξ(y))2 + = h2 /2, (5) with adequate boundary conditions. The solution Uε : X → R to this system of equations is computed in a single pass with O(N ln N) complexity [17], using a variant of the Fast-Marching algorithm. This is possible since the l.h.s. of (5) is a non-decreasing function of the positive parts of the finite differences (Uξ(x) − Uξ(y))y∈X. Note that the eikonal PDE discretization (5), based on upwind finite differences, differs from the semi-Lagrangian approach [18], which can also be solved in a single pass but is usually less efficient due to the large cardinality and radius of its stencils. Image segmentation techniques relying on the numerical solutions to anisotropic eikonal PDEs were proposed in [6] using Riemannian metrics, and in [4, 7] based on the reversible Reeds-Shepp car and Euler elastica curvature penalized models respectively. However these early works rely on non- causal discretizations, which have super-linear complexity O(N1+1/d ) where the unspecified constant is large for strongly anisotropic and non-uniform metrics. This alternative approach yields (much) longer solve times, incompatible our application - where strongly anisotropic three dimensional eikonal PDEs are solved as part of an inner loop of an optimization procedure. To be able to use the gradient to solve the problem (1), we need to differ- entiate the cost C(ξ, Γ) w.r.t. the first player strategy ξ ∈ Ξ. In view of (3), this only requires the sensitivity of the discrete solution values Uξ(x∗) at the few points x∗ ∈ X ∩ Θ, w.r.t to variations in the weights cξ(x, y), x, y ∈ X. For that purpose we differentiate (5) w.r.t. ξ at an arbitrary point x ∈ X \ Υ, and obtain X y∈X ωξ(x, y) (dUξ(x) − dUξ(y) + (Uξ(x) − Uξ(y)) d ln cξ(x, y)) = 0, where ωξ(x, y) := c2 ξ(x, y)(Uξ(x) − Uξ(y))+. Therefore dUξ(x) = X y∈X αξ(x, y)dUξ(y) + X y∈X βξ(x, y)dcξ(x, y), (6) where αξ(x, y) := ωξ(x, y)/ P y ωξ(x, y), and βξ(x, y) := αξ(x, y)/cξ(x, y). We first choose x = x∗ in (6), and then recursively eliminate the terms dUξ(y) by applying the same formula at these points, except for points in the source set y ∈ Υ for which one uses the explicit expression dUξ(y) = 0 (since Uξ(y) = 0 is in this case independent of ξ). This procedure terminates: indeed, whenever dUξ(x) depends on dUξ(y) in (6), one has αξ(x, y) > 0, thus ω(x, y) > 0, hence Uξ(x) > Uξ(y). It is closely related to automatic differentiation by reverse accu- mulation [11], and has the modest complexity O(N). 4 Numerical results The chosen physical domain R is the rectangle [0, 2]×[0, 1] minus some obstacles, as illustrated on Figure 1. Source point is (0.2, 0.5) and target keypoint (1.8, 0.5). The computational domain is thus Ω = R for curvature independent models and Ω = R × S1 for curvature dependent models, which is discretized on a 180 × 89 or 180 × 89 × 60 grid. No intervention from the first player. The cost function is Cξ(p, ṗ, p̈) = C∗(|p̈|), where C∗(κ) is respectively 1, p 1 + ρ2κ2 and (1 if ρκ ≤ 1, otherwise +∞), with ρ := 0.3. The differences between the three models are apparent: the curvature independent model uses the same path forward and back; the Reeds-Shepp car spreads some curvature along the way but still makes an angle at the target point; the Dubins car maintains the radius of curvature below the bound ρ, and its trajectory is a succession of straight and circular segments. A referee notes that following an optimal trajectory for the Dubins model is dangerous in practice, since any small deviation is typically impossible to correct locally, and may drive into an obstacle; these trajectories are also easier to detect due to the large circular arc motions. Curvature independent Reeds-Shepp car, forward only Dubins car Fig. 1. Shortest path from the blue point (left) to the red keypoint (right) and back. Next we study three games where player one aims to detect player two along its way from the source set Υ to the target Θ and back, using different means. If the first player does not intervene, see Figure 1, or if its strategy is not optimized, see Figure 3, then there is typically a unique optimal path (optimal loop in our games) for player two. In contrast, an interesting qualitative property of the optimal strategy ξ ∈ Ξ for the first player is that it has a large number of optimal responses from player two, see Figure 4, in some cases even a continuum, see Figure 2 (bottom) and [5]. This is typical of two player games. Fresh paint based detection. In this toy model, see Figure 2, the first player spreads some fresh paint over the domain, and the second player is regarded as detected if he comes back covered in it from his visit to the keypoint. The cost function is Cξ(p, ṗ, p̈) = ξ(p)C∗(|p̈|), where ξ : R → R+ is the fresh paint density, decided by the first player, and C∗(κ) is as above. For wellposedness, we impose upper and lower bounds on the paint density, namely 0.1 ≤ ξ(p) ≤ 1, and subtract the paint supply cost R R ξ(p)dp to (1). The main interest of this specific game, also considered in [5], is that C(ξ, Γ) is concave w.r.t. ξ ∈ Ξ. The observed optimal strategy for player Ξ is in the curvature independent case to make some “fences” of paint between close obstacles, and in the curvature penalized models to deposit paint at the edges of obstacles, as well as along specific circular arcs for the Dubins model. Visual detection. The first player places some cameras, e.g. with 360-degree field of view and mounted at the ceiling, which efficiency at detecting the second Curvature independent Reeds-Shepp car, forward only Dubins car Fig. 2. Top: Optimal distribution of paint, to mark a path from the blue point (left) to the red keypoint (right) and back. Bottom: Geodesic density at the optimal paint distribution. player decreases with distance and is blocked by obstacles, see Figure 3. The cost function is Cξ(p, ṗ, p̈) = C∗(κ) X q∈ξ [p,q]⊆R 1 kq − pk2 , (7) where ξ ∈ Ξ is a subset of R with prescribed cardinality, two in our experiments. The green arrows on Fig 3 originate from the current (non optimal) camera position, and point in the direction of greatest growth ∇C(ξ, Γ) for the first player objective function. Curvature independent Reeds-Shepp forward Dubins car Curvature independent Reeds-Shepp forward Dubins car Fig. 3. Field of view of the cameras (black gradients), optimal furtive paths (red lines), local direction of improvement of the camera position (green arrows). Radar based detection. The first player places some radars on the domain R = [0, 2] × [0, 1], here devoid of obstacles, and the second player has to fly by unde- tected. The cost function is Cξ(p, ṗ, p̈) = C∗(|p̈|) v u u t X q∈ξ hṗ, npqi2 + δ2hṗ, n⊥ pqi kp − qk4 (8) where npq := (q−p)/kq−pk. The first player strategy ξ contains the positions of three radars, constrained to lie in the subdomain [0.4, 1.6]×[0, 1]. The parameter δ is set to 1 for an isotropic radar cross section (RCS), or to 0.2 for an anisotropic RCS. In the latter case a plane showing its side to radar is five times less likely to be detected than a plane showing its nose or back, at the same position. Green arrows on Figure 4 point from the original position to the (locally) optimized position for player Ξ. At this position, several paths are optimal for player Γ, shown in red on Fig 4. Curvature independent Reeds-Shepp forward Dubins car Curvature independent Reeds-Shepp forward Dubins car Fig. 4. Optimal radar placement with an isotropic (top) or anisotropic (bottom) radar cross section. Computational cost On a standard Laptop computer (2.7Ghz, 16GB ram), op- timizing the second player objective, by solving a generalized eikonal equation, takes ≈ 1s in the curvature dependent case, and ≈ 60 times less in the curvature independent case thanks to the absence of angular discretization of the domain. Optimizing the first player objective takes ≈ 100 L-BFGS iterations, each one taking at most 8s. For the stability of the minimization procedure, the problems considered were slightly regularized by the use of soft-minimum functions and by “blurring” the target keypoint over the 3 × 3 box of adjacent pixels. 5 Conclusion We have modeled a motion planning problem that minimize an anisotropic prob- ability of detection, taking into account navigation constraints while computing the gradient of the value function related to the sensors location. This model is thus useful for surveillance applications modeled as a two-player zero-sum game involving a target that tries to avoid detection. References 1. F. Barbaresco and B. Monnier, Minimal geodesics bundles by active contours: Radar application for computation of most threathening trajectories areas & corridors, 10th European Signal Processing Conference, Tampere, 2000, pp. 1-4. 2. F. Barbaresco, Computation of most threatening radar trajectories areas and cor- ridors based on fast-marching & Level Sets, IEEE Symposium on Computational Intelligence for Security and Defense Applications (CISDA), Paris, 2011, pp. 51-58. 3. M. Bardi and I. Capuzzo-Dolcetta, Optimal control and viscosity solutions of Hamilton-Jacobi-Bellman equations, Bikhauser, 1997. 4. E. Bekkers, R. Duits, A. Mashtakov, and G. Sanguinetti, Data-Driven Sub- Riemannian Geodesics in SE(2). In Proc. SSVM 2015, pages 613-625, 2015. 5. F. Benmansour, G. Carlier, G. Peyré, and F. Santambrogio, Derivatives with respect to metrics and applications: subgradient marching algorithm, Numerische Mathe- matik, vol. 116, no. 3, pp. 357-381, May 2010. 6. F. Benmansour and L. Cohen, Tubular Structure Segmentation Based on Minimal Path Method and Anisotropic Enhancement. International Journal of Computer Vi- sion, 92(2), pages 192-210, 2011. 7. D. Chen, J.-M. Mirebeau and L. D. Cohen, A New Finsler Minimal Path Model with Curvature Penalization for Image Segmentation and Closed Contour Detection, in Proc. CVPR 2016, Las Vegas, USA. 8. L. E. Dubins, On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents, Amer. J. Math., vol. 79, pp. 497-516, 1957. 9. R. Duits, S.P.L. Meesters, J.-M. Mirebeau, J. M. Portegies, Optimal Paths for Variants of the 2D and 3D Reeds-Shepp Car with Applications in Image Analysis, Preprint available on arXiv 10. J. Fehrenbach and J.-M. Mirebeau, Sparse Non-negative Stencils for Anisotropic Diffusion, Journal of Mathematical Imaging and Vision, pp. 1-25, 2013. 11. A. Griewank, A. Walther, Evaluating derivatives: principles and techniques of al- gorithmic differentiation. Society for Industrial and Applied Mathematics, 2008. 12. J.-M. Mirebeau, Anisotropic Fast-Marching on cartesian grids using Lattice Basis Reduction, SIAM J. Numer. Anal., vol. 52, no. 4, pp. 1573-1599, Jan. 2014. 13. J.-M. Mirebeau, Minimal stencils for discretizations of anisotropic PDEs preserv- ing causality or the maximum principle, SIAM J. Numer. Anal., vol. 54, no. 3, pp. 1582-1611, 2016. 14. J.-M. Mirebeau, Anisotropic fast-marching on cartesian grids using Voronois first reduction of quadratic forms, preprint available on HAL 15. J.-M. Mirebeau, Fast Marching methods for Curvature Penalized Shortest Paths, preprint available on HAL 16. D. Mumford, Elastica and Computer Vision, no. 31, New York, NY: Springer New York, pp. 491-506, 1994. 17. E. Rouy and A. Tourin, A Viscosity Solutions Approach to Shape-From-Shading, SIAM J. Numer. Anal., vol. 29, no. 3, pp. 867-884, Jul. 1992. 18. J. Sethian and A. Vladimirsky, Ordered upwind methods for static Hamilton-Jacobi equations. Proceedings of the National Academy of Sciences, 98(20), 2001. 19. A. Schürmann, Computational geometry of positive definite quadratic forms, Uni- versity Lecture Series, 2009. 20. C. Strode, Optimising multistatic sensor locations using path planning and game theory, IEEE Symposium on Computational Intelligence for Security and Defense Applications (CISDA), Paris, 2011, pp. 9–16.