Jami Torki

Ali Sari

Charles

Charles Joubert

Symposium on Materials for Energy Storage and Conversion Akyaka, Mugla

come

Modeling the distribution of ions in a electrolytic capacitor

The increase of the density of electrical energy stored in energy storage systems is a crucial element in reducing the volume of on-board systems. Capacitors are made up of two conductive armatures, called electrodes, separated by a dielectric allowing opposing electric charges to be stored on each electrode. This increase can be visualized through an increase in the surface area of the electrodes (engraved, porous [START_REF] Simon | Materials for electrochemical capacitors[END_REF] or based on nanotubes interdigitated principle [START_REF] Chiou | Cycle stability of the electrochemical capacitors patterned with vertically aligned carbon nanotubes in an LiPF6-based electrolyte[END_REF]). This work aims to develop better generation of energy storage system using concepts based on aluminum electrolytic capacitors (AEC) principle (Figure 1 [START_REF] Hayek | Analyse de vieillissement, estimation de la durée de vie et méthode de surveillance de l'état de santé des condensateurs électrolytiques[END_REF]). The modeling will allow us to understand the physics involved in an energy storage mechanism: an electrolytic capacitor. One of the current highlights of my research is to study the density distribution of ionic species in the electrolyte in an electric field. This understanding will allow us to achieve an optimization of macroscopic parameters such as the capacitance, the Equivalent Series Resistance (ESR) or the temporal response of the system which can be deduce from microscopic analysis. To introduce my modeling, AEC consist of two aluminum electrodes (anode and related cathode), an oxide layer acting as a dielectric and an electrolyte usually boric acid or benzoic dissolved in ethylene glycol or glycerin and water [START_REF] Hayek | Analyse de vieillissement, estimation de la durée de vie et méthode de surveillance de l'état de santé des condensateurs électrolytiques[END_REF]. These chemical mixtures allow the electrolyte to be a very good conductor, over a wide temperature range, in order to bring the negative potential of the cathode (real cathode) closer to the dielectric via ionic transport in the electrolyte. The sheets of electrolyte impregnated paper serve as physical separators limiting the risk of short circuits between the two electrodes. The thin layer of Al2O3 (alumina) is caused by natural oxidation of aluminum which limits corrosion.

My current research consists of studying the distribution of ion concentration, potential and electrical field in the electrolyte for different models at steady state first and then in transient regime. Indeed, the positively or negatively charged ions will be attracted by the electric potential of the electrodes. In addition, they will in turn influence the electric field. This is the coupling between the spatial charge density with the electric potential and this is why the equations for describing the interaction between the ions and the electric field is nonlinear [START_REF] Chien | Modeling of Ion Diffusion Induced Current through a Nanochannel[END_REF]. My modeling is based on microscopic systems with simplified geometries in order to facilitate the resolution of this type of non linear equation. The first model will be a theoretical 2D model with face-to-face plane electrodes in order to compare my model with analytical solutions and check the consistency of the simulation results. The second case will be to study the impact of an etching similar to that of an anode of an AEC (Figure 1). And finally, we will consider a particular geometry allowing, like the etching, an increase in the surface of the electrode. Therefore, we can compare the results and highlight the advantages / disadvantages of each structure. The analysis of the transient regime in these three models makes it possible to study the behavior of the ions subjected to the strong electric field during the charge of the AEC. The temporal evolution of the relaxation resulting from the diffusion of the ions during the discharge of the electrolytic capacitor can also be simulated in order to characterize the electrolytic capacitors operating.

Figure 1 .

 1 Figure 1. Aluminum electrolytic capacitor.


Gabriel Peyré
email: peyre@ceremade.dauphine.fr

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

Heuristically Driven Front Propagation for Fast Geodesic Extraction

Keywords: geodesics, Fast Marching, shortest paths, heuristic, A * algorithm

This paper presents a new method to quickly extract geodesic paths on images and 3D meshes. We use a heuristic to drive the front propagation procedure of the classical Fast Marching. This results in a modification of the Fast Marching algorithm that is similar to the A * algorithm used in artificial intelligence. In order to find very quickly geodesic paths between any given pair of points, two methods are proposed to devise an heuristic that restrict the front propagation. The multiresolition heuristic computes the heuristic using a propagation on a coarse map. For applications where pre-computation is acceptable, the landmark-based heuristic pre-computes distance maps to a sparse set of landmark points. We introduce various distortion metrics in order to quantify the errors introduced by the heuristically driven propagations. We show that both heuristic approaches bring a large speed-up for large scale applications that require the extraction of geodesics on images and 3D meshes.

In order to compute the geodesic path between two points on the surface, computation is needed on a large part of the surface when using classical fast marching (left). When the distance map to a set of landmarks is pre-computed, and the propagation is heuristically driven by these maps, only the colored region is explored, and it becomes smaller and smaller as the number of landmarks increases.

Previous works

The extraction of shortest paths is a building block for a large class of applications ranging from graph theory to computer vision. The computation can be carried in a totally discrete setting (such as a graph) or over a discretized domain (such as an image or a 3D mesh).

Graphs and discrete computation. The canonical method to compute shortest paths on graphs is the Dijkstra algorithm (see for instance [START_REF] Cormen | Introduction to Algorithms[END_REF]). Fast exploration strategies have been used to speed-up the computation. The A * algorithm [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF] makes use of a heuristic to reduce the search space. Other exploration strategies have been proposed in the field of artificial intelligence, such as IDA * [START_REF] Korf | Depth-first iterative-deepening: an optimal admissible tree search[END_REF].

Images and continuous setting. In order to extract geodesics for a continuous metric, the Fast Marching algorithm [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF] uses a front propagation to solve non-iteratively a finite difference approximation of the Eikonal equation. A similar algorithm was also proposed in [START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF]. The minimal length properties of geodesics has been applied in computer vision, for example to solve global minimization problems for deformable models [START_REF] Cohen | Global Minimum for Active Contour models: A Minimal Path Approach[END_REF]. The continuous nature of this method is particularly attractive for image processing, for instance to extract tubular structures and centerlines in 3D medical data [START_REF] Deschamps | Fast Extraction of Minimal Paths in 3D Images and Applications to Virtual Endoscopy[END_REF].

Path planning: from discrete to continuous. Discrete computation on graphs gives rise to numerous schemes to perform motion planning and the A * algorithm is extensively used for path-finding in video games [START_REF] Stout | Smart moves: Intelligent path-finding[END_REF]. For the case of Euclidean metrics, faster algorithms have been devised that exploit specific data structures such as visibility graphs [START_REF] Rohnert | Shortest paths in the plane with convex polygonal obstacles[END_REF]. The Fast Marching algorithm can be used to extract paths with a non-Euclidean metric [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. The authors of [START_REF] Melchior | Consideration of obstacle danger level in path planning using A* and fast-marching optimisation: comparative study[END_REF] compare the Fast Marching and the A * algorithms to perform motion planning, but they do not propose a heuristic modification of the Fast Marching.

Shortest Path: Continuous and Discrete Algorithms

In this section we develop a unifying framework which includes the Fast Marching [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], Dijkstra [START_REF] Cormen | Introduction to Algorithms[END_REF] and A * [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF] algorithms together with our new heuristically driven propagation.

Front Propagation Methods for Shortest Path

We work on a discrete set of points and for each point x we have access to its neighbors y, defining the relation y ∼ x. These points can be embedded in discretete grid (for Fast Marching and for our method) or can be the vertices of a graph (for Dijkstra and A * ). Our goal is to compute the distance function U (x) = d(x 0 , x) to some starting point x 0 .

Initialization:

• Alive set: the starting point x 0 ; • Trial set: the neighbors of x 0 ; • Far: the set of all other grid points.

Loop:

• Let x be the Trial point with the smallest priority P(x); • Move it from the Trial to the Alive set;

• For each neighbor y of the current point x:

-if y is Far, then add it to Alive and compute a new value for U (y), -if y is Alive, recompute the value U (y), and update it if the new value is smaller, -recompute the priority P(y). • If the end point x = x 1 is reached, stop the algorithm.

Table 1: Pseudo-code for the common framework for front propagation.

The propagation methods label the points during the computation according to: • Alive is the set of points at which the distance value U has been computed and will not change; • Trial is the set of next grid points to be examined and for which an estimate of U has been computed; • Far is the set of all other grid points, for which there is not yet an estimate for U . Table 1 shows the main steps of the algorithms. Each algorithm must implement the following subfunctions

• A way to update the value U (y) at a given Trial point y. This computation uses the values of U at the adjacent locations. This computation depends on the metric used, which can be defined on a graph (for discrete methods) or on the whole space (for continuous methods). We explain in sections 2.2 and 2.3 specific instantiations for the Dijkstra and Fast Marching algorithms. • A priority map P orders the set of Alive points according to some computational criterion. In the Fast Marching and Dijkstra algorithm, P(x) = U (x) is the current distance to the starting point. In our heuristical front propagation as in the A * algorithm, P(x) is chosen to minimize the number of visited points. We explain in sections 3 to 5 how to actually construct a priority function P that makes use of a heuristic.

Discrete Case: Dijkstra

In the discrete setting, a symmetric weight g(x, y) is used to define the distance between two adjacent points x ∼ y of the graph. The length of a path

v = (v 1 ∼ . . . ∼ v m ], of m adjacent points v i is L(v) = g(v i , v i+1
), and we define the distance between two vertices

d(x 0 , x 1 ) = min v {L(v) \ v 1 = x 0 , v m = x 1 } .
The distance U (x) at a vertex x in the alive set is updated during the propagation according to

U (x) = min y∼x U (y) + g(x, y) .
The shortest path v from x 0 to x 1 is tracked backward using

v 0 = x 1 and v i+1 = argmin y∼vi U (y).

Continuous Case: Fast Marching

In R d , we are given a potential function g(x) > 0, and the weighted geodesic distance between two points x 0 , x 1 ∈ R d , is defined as

d(x 0 , x 1 ) = min γ 1 0 ||γ ′ (t)|| g(γ(t))dt , (1)
where γ is a piecewise regular curve with γ(0) = x 0 and γ(1) = x 1 . When g = 1, the integral in (1) corresponds to the length of the curve γ and d is the classical Euclidean distance.

The Fast Marching method uses the fact that the function U satisfies the nonlinear Eikonal equation:

||∇U (x)|| = g(x). (2)
The distance U (x) = u at a point x = x i,j in the trial set is updated during the propagation according to the solution of

max(u -U (x i-1,j ), u -U (x i+1,j ), 0) 2 + max(u -U (x i,j-1 ), u -U (x i,j+1 ), 0) 2 = h 2 g(x i,j ) 2 .
This is a second order equation (the equation is written in R 2 for simplicity) and it can be solved as detailed for example in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF].

The geodesic curve γ from x 0 to x 1 can be computed by extracting the parametric curve C(t) that solves the back propagation equation:

dC dt = - --→ ∇U with C(0) = x 1 .
This gradient descent is a local computation, and it only uses the value of U for a small fraction of the visited grid points. Note that these grid points are those located in the Alive set at the end of the front propagation procedure.

Heuristically Driven Front Propagation

In this section we explain our algorithm in the 2D setting, and show some numerical results that illustrate the main features of this method.

No heuristic, 3D view of the potential

Propagation with a Heuristic

In order to minimize the number of Alive points at the end of the front propagation procedure, one should use a priority function P that tries to advance the front toward the goal point x 1 . In order to do so, we assume that, together with the current weighted distance to the start point U (x) = d(x 0 , x), we have an estimate of the remaining weighted distance V (x) ≈ d(x 1 , x). Our heuristical front propagation algorithm follows the implementation of table 1 with a priority map

P(x) = U (x) + V (x). (3)
The rationale behind the definition of P is that d(x 0 , x) + d(x 1 , x) is minimal and constant along the geodesic path joining x 0 and x 1 , see [START_REF] Kimmel | Finding shortest paths on surfaces using level sets propagation[END_REF].

A * Algorithm. For discrete graphs, it has been shown that if the heuristic satisfies V (x) d(x 1 , x), then the extracted geodesic is a path of minimum length. This leads to the A * algorithm of [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF]. Various strategies have been proposed to devise admissible heuristics, see [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] and the references therein.

Heuristically driven Fast Marching

In figure 2, one can see a front propagation, where we have used the oracle heuristic V (x) = λ d(x 1 , x), with a parameter λ ∈ [0, 1]. The value λ = 0 corresponds to the classical Fast Marching propagation, which results in a large region of Alive points (colored in red). However, as we increase the value of λ toward 1, the explored region shrinks around the geodesic path that links x 0 to x 1 .

There are however two important issues with this ordering of the Trial set: • This ordering can break the monotone condition that is required by the Fast Marching algorithm to produce a valid approximation of the continuous underlying distance function. We show in the numerical results presented in sub-section 5.2 that the geometric error on the extracted geodesic remains low. • We do not have immediate access to the remaining distance d(x, x 1 ), since it would involve performing another front propagation from x 1 . We explain in sections 4 and 5 two methods to overcome this problem.

Evaluation of a Heuristic

We cast the problem of finding a good heuristic into the problem of approximating the distance function d(x, y) between two points (x, y) by some function d. We define the heuristic using

V (x) = d(x 1 , x) ≈ d(x 1 , x).
This approximated distance d should satisfy d d in order not to perturb the propagation from the true shortest path (computed with V = 0). It must be fast to evaluate and can use a reasonable amount of pre-computed data.

Section 4 uses a propagation on a coarse map in order to compute such an approximate d. This multiresolution heuristic is fast to compute and does not require the storage of additional data. Section 5 uses pre-computed distances to landmark points to compute d. This landmark-based heuristic requires additional data but leads to highly accurate estimation of the real distance d while enforcing d d.

To evaluate the quality of the heuristic, we propose two metrics • The approximation metric

E 1 (d, d) = |d(x, y) -d(x, y)| 2 dx dy,
is approximated during evaluation using a finite number of precomputed distance maps d(p i , x) to a set of points {p i } m i=1 chosen at random

E 1 (d, d) = 1 m m i=1 |d(p i , y) -d(p i , y)| 2 dy.
• The computation gain metric measures the usefulness of the heuristic for extracting geodesics between a set of pairs of points {(x i 0 , x i 1 )} i . For each front propagation from x i 0 to x i 1 , we measure the area A i ( d) of the Alive set at the end of the propagation. We also measure the area A i (0) of the Alive set for a propagation without the heuristic. We define the computation gain metric to be

E 2 (d, d) = 1 m m i=1 A i ( d)/A i (0).
When the geodesics queries are random, these two metrics tend to be the same. However, the computationgain metric is able to better adapt to typical applications where queries can be highly non-uniform (for instance in road or tubular structure extraction).

4 Multiresolution Heuristic

Coarse-to-fine Geodesic Computations

In order to compute the remaining distance V (x) ≈ d(x, x 1 ) with a fast algorithm, we perform a Fast Marching front propagation starting from the point x 1 , but on a coarser grid. We thus have introduced a second parameter for our heuristical front propagation: the resolution R ∈ (0, 1) we use for the coarse grid.

If the original potential map g is of size n × n, the query of P(y) thus requires:

• The pre-computation of a coarse potential map g R of size (Rn) × (Rn). This is done by first pre-filtering g (to avoid aliasing of high frequencies) and then applying a cubic spline re-interpolation on a coarser grid. • The pre-computation of the approximate distance map V of size (Rn) × (Rn) is done by performing a full Fast Marching on a coarse grid, using potential g R , and starting from point x 1 . • During the heuristical front propagation starting from point x 0 , when P(y) is queried, the value of V is interpolated with cubic splines on the coarse grid to retrieve a value on the original grid. In figure 3 one can see the coarse map g R , and high frequency details such as small roads disappear as the resolution gets too low.

There is clearly a tradeoff between choosing a low R to reduce the computation time, and a high R so that V (x) approximates d(x, x 1 ) well. The new algorithm we propose allows us to use multiresolution computation for the extraction of geodesic curves. Using a multiresolution framework for solving the point-to-point geodesic problem is not so easy because it is a boundary problem, and for instance, multigrid methods are not suitable. Adaptive mesh [START_REF] Droske | An adaptive level set method for interactive segmentation of intracranial tumors[END_REF] and multigrid methods [START_REF] Papandreou | A fast multigrid implicit algorithm for the evolution of geodesic active contours[END_REF] have been used in conjunction with geodesic active contours for segmentation. In order to estimate the precision of the results, we use the Hausdorff error between the paths obtained by fast marching with and without the multiresolution heuristic. In figure 6 one can see the geodesics extracted for different values of λ. Figure 4 shows the result of our algorithm for various settings on (a) a synthetic map and (b) a satellite image. We have depicted:

=100% R =30% R =15% R

Numerical Validation of Multiresolution Heuristic

• The 2D map: The red curves indicate the boundary of the visited region. One can see that these curves shrink toward the geodesic (central blue curve) as one increase the strength of the heuristic from 0% to 100%. • Hausdorff error vs. heuristic strength λ: We have set the heuristic resolution R to 50%. One can see that the error is higher for the synthetic map (a). This is due to the fact that this map contains large flat areas, where a small error in the computed geodesic distance leads to deviation of the extracted geodesic.

In contrast, the geodesic in the satellite image (b) contains very anisotropic areas, which stabilize the extracted geodesic. • Hausdorff error vs. heuristic resolution R: We have set the heuristic strength λ to 50%. One can see that the synthetic map (a) is nearly insensitive to the resolution of the coarse map used to compute the heuristic. This is because the underlying function is very smooth, so one can reduce the resolution a lot without too much impact on the accuracy of the heuristic. In contrast, one can see that the satellite image suffers from large variation when the resolution parameter R becomes smaller than 60% and then again for 90%. This is due to strong topological changes in the path, as depicted in figure 5. For R = 20%, this overhead is balanced by the heuristic saving as soon as λ 5%. These tests clearly show that our algorithm can bring a large computational speed-up, but the parameters should be finely tuned to adapt to the characteristic of each map. For instance, these experiments show that the user must have some prior knowledge about the typical width of the tubular structures he wants to extract, and set the resolution R so that the coarse map V R still contains these structures.

Applications of a Multiresolution Heuristic

Volumetric Geodesic Extraction 3D geodesic extraction is very useful in medical volumetric data analysis. It can be applied to perform tubular structure extraction, and it is extended to virtual endoscopy in [START_REF] Deschamps | Fast Extraction of Minimal Paths in 3D Images and Applications to Virtual Endoscopy[END_REF]. In figure 7 one can see the extraction of 3D geodesics from synthetic data (top and middle rows) and from real medical data (bottom row) for R = 20%. The red surface shows the boundary of the explored regions of alive cells. The computational time gain (Comp. gain) is also indicated. Globally Optimal Geodesic Active Contours The concept of circular geodesics was first introduced in [START_REF] Sun | Circular shortest path in images[END_REF]. The authors of [START_REF] Appleton | Globally optimal geodesic active contours[END_REF] proposed a simple way to compute circular geodesics around a point, in order to compute a globally optimal geodesic, with an application to object segmentation. The user simply select a point C inside the object to segment (see figure 8), and then the algorithm virtually "cuts" the image along a horizontal line that links C to the boundary of the image. This way, one can force a geodesic path to go around C by running a classical Fast Marching from a point S to itself, but forbidding the front to pass through the segment CD. For an underlying image I, the globally optimal geodesic around C is defined as the closed geodesic curve with minimum length, where the metric is defined as

g(x) = 1 ||C -x|| 1 1 + ||∇I(x)|| 2 + ε, (4)
where ||C -x|| is the distance from the curve point x to the center C. The authors of [START_REF] Appleton | Globally optimal geodesic active contours[END_REF] proposed a powerful algorithm based on the branch-and-bound paradigm, which is a binary search that avoids computing the closed geodesic for each point S on the segment CD. However, with our heuristic front propagation, we have tested a simpler algorithm that works well in practice. We simply compute the circular geodesics that pass though a given fixed number of points along the cut segment CD. These extractions can be performed quickly using our heuristically driven front propagation, with the restriction that the front should not pass though the cut segment.

In figure 9, we have shown a globally optimal circular geodesic, computed with various heuristic strengths λ.

Heuristic 0%

40% 100% 60% 80%

Figure 9: Globally optimal circular path extraction with increasing heuristic.

Geodesic Extraction on 3D Meshes The Fast Marching algorithm has been extended to 3D meshes in [START_REF] Sethian | Computing Geodesic Paths on Manifolds[END_REF]. Our heuristic algorithm also extends to 3D meshes, with the following modifications with respect to the Euclidean setting:

• We must construct a coarse mesh approximation of the original 3D mesh. Mesh simplification is a large topic, and several greedy methods exist, see for example [START_REF] Hoppe | Progressive meshes[END_REF]. In our tests, we use the farthest point strategy proposed in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF] for remeshing, since it uses Fast Marching as a building block. • Once the heuristic function has been computed on the coarse mesh, it must be interpolated on the original dense mesh. Several methods for data interpolation on 3D meshes exist, and we have used a method derived from harmonic mesh parameterization [START_REF] Floater | Parameterization of Manifold Triangulations[END_REF]. This involves the solution of a sparse linear system that describes a harmonic function that fits the values computed on the coarse mesh. These two steps are quite computationally intensive, but note that:

• The coarse mesh can be pre-computed, and can be re-used for multiple geodesic extraction.

• To avoid the computational overhead of computing the interpolation on the whole mesh, we use the local parameterization strategy of [START_REF] Surazhsky | Isotropic remeshing of surfaces: a local parameterization approach[END_REF]. We compute the interpolation only on a small set of overlapping disklike charts that cover the region of alive vertices. In figure 1 and 10, one can see the algorithm in action on various meshes, and for various values of the parameter λ.

Coarse mesh Heuristic 0% 50% 100% 90%

Figure 10: Heuristically driven front propagation on 3D meshes with a multiresolution heuristic.

Landmark-based Front Propagation

In this section, we describe an alternate way to implement the heuristic. This heuristic is computed using an approximated distance d evaluated with a set of pre-computed distances to landmark points.

Landmark-based heuristic

A new method for distance evaluation on a graph has been introduced in [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] as an admissible heuristic for the A * algorithm. We explain why this method can be useful for our heuristically driven Fast Marching and give a quantitative numerical study.

This method exploits the triangle inequality, since we have, for every pair of points x and y, d(x, y) = sup

In the following, we drop the z 1 , . . . , z n dependencies and call the approximated distance d. We note that this approximation always satisfies the condition d d.

Figure 11 gives an intuitive explanation of the efficiency of this approximation. In the ideal case, the geodesic γ k joining a landmark z k to x also passes through y. In this case, we have d(x, y) = d(x, y) and there is no approximation. But in most cases, this is not true, but a geodesic γ k passes close to y and d is indeed a good approximation of the real distance.

γ γ k γ k x x z k y y y k Ideal case Real life case d ( x , y ) d ( y , z k ) - d ( x , z k ) = d ( x , y ) d ( y , z k ) - d ( x , z k ) >
Figure 11: Justification of the approximation properties of d.

Landmark-based Path Extraction

This landmark-based propagation can be used to speed-up the computation of most path planing applications, including path finding in video games or robotic path planing.

No heuristic 1 landmark 3 landmarks 12 landmarks 20 landmarks In figure 12, one can see the active region explored by the propagation algorithm for an increasing number of landmark points. These points are chosen at random on the 2D image. The explored region progressively shrinks toward the central extracted curve when we add more points, since the heuristic is becoming more accurate.

Hausdforff error

Hausdforff error In figure 13, one can see a numerical evaluation of the precision of the extracted geodesic. It reports the distortion ||γ -γ|| H between the original geodesic γ (computed without the heuristic) and γ (computed with the heuristic). We use the symmetric mean-square Hausdorff metric

||γ -γ|| 2 H = 1 2 γ min y∈γ ||x -y|| 2 2 dx + γ min x∈γ ||x -y|| 2 2 dx ,
since this captures the geometric distortion caused by our partial propagation. On 2D maps containing important curvilinear features (such as the road in the example on the left), the distortion caused by the heuristic is nearly unnoticeable. On the contrary, for relatively flat maps (such as the one depicted on the right) the distortion can be relatively high (about few pixels for a map of 256 × 256 pixels). This is because the salient features of this map catch the geodesic and avoid large lateral moves when the front propagation is modified.

In order to improve the quality of the heuristic and thus the resulting speed-up, one has to carefully seed the landmarks across the image. In the next section we setup a complete framework for the evaluation of this sampling, together with strategies to get a high-quality seeding of the base points.

Seeding strategies for the landmarks

The quality of a sampling can be measured using the previously introduced metrics

E i (z 1 , . . . , z n ) = E i (d, dz1...zn ) for i = 1, 2.
In order to find sampling locations {z 1 , . . . , z n } for the landmarks, we use various strategies, among which: (a) A manual sampling that exploits specific knowledge of the 2D map g. This semi-automatic method is not studied in this paper. (b) A random sampling in the image. (c) A uniform sampling according to the distance d. This can be accomplished using the farthest point sampling procedure proposed in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF], where one chooses

z n+1 = argmax z ( argmin
1 k n d(z k , z)).
(d) A uniform seeding on the boundary of the domain.

(e) A greedy strategy minimizing the E 1 metric. We define the locations recursively

z n+1 = argmin z∈R E 1 (z 1 , . . . , z n , z),
where R is some small set of random candidates. (f) Same as (e) but using E 2 instead of E 1 .

Distance approximation evaluation.

A good seeding strategy should be able to put a landmark before any couple of points in the domain, near the geodesic path connecting these points. For a constant (i.e. Euclidean) metric, (d) is thus the optimal sampling scheme. However, for a complex metric, this is no more the case.

In figure 14, one can see the approximated distance d(x 1 , x) to a target point x 1 , computed for n = 32 landmarks using various seeding strategies. In figure 15, one can see a plot of the approximation error E 1 (x 1 , . . . , x n ). The error-driven seeding strategy (e) clearly performs the best as expected by its definition. For a smooth potential g (left) the seeding strategy on the boundary (d) performs well, but it fails to capture the topology of complex 2D maps (right).

Computational saving evaluation. In order to evaluate our algorithm in a real setting, we have set-up a complete framework using 500 typical queries {(x i 0 , x i 1 )} i for various 2D maps. Some resulting shortest paths queries are shown on the left of figure 16, they are relevant for many standard applications such as path finding in video games or robot path planing. The right of figure 16 shows how the computational metric E 2 (x 1 , . . . , x n ) decreases for several sampling strategies. For greedy sampling strategy (f), we have used another set of 100 typical queries, in order not to bias the computation of the error E 2 .

One can see that the computation-driven sampling strategy (f) clearly outperforms the other strategies. In particular, the approximation-based strategy (e) does not give good results, mainly because typical paths are not correlated to areas where the approximation error ||d -d|| is large.

Reducing Memory Usage

Cells Representation Classical methods, such as using an octree data structure, can be used to reduce the memory usage of level set algorithms, for example in order to perform image segmentation [START_REF] Droske | An adaptive level set method for interactive segmentation of intracranial tumors[END_REF].

We choose to implement a simple data structure to reduce the memory usage by allocating the grid cell on the fly during the propagation. A typical cell data structure, for 2D computation, is: struct cell { // current geodesic distance double distance; // either alive, trial or far char state; // pointers to the 4 neighbors cell * neighbors [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF]; };

To be able to retrieve a given cell in constant time, we also store the list of allocated cells in a hash table.

This is important because when a new cell is allocated, we need to connect it to the existing cells. Figure 17 shows a graphical display of the overall data structures. This pointer-based representation of the neighboring relation is convenient to extract the geodesic with a gradient descent. There is some memory overhead due to the explicit storage of pointers to neighbors, but the fact that our scheme explores significantly fewer cells than classical Fast Marching allows to save much more memory, as shown in the next section. The computing time overhead due to the use of a hash table is about 40% in all our tests. • Be asymetric, since we do not care about the compression time, but we need fast access to distance values.

• Give low decompression error, since we need an accurate heuristic and we need to satisfy the condition d d as much as possible.

We use a two-step procedure:

• Differential representation: A 1 (x) = d(x 1 , x) and

A k+1 (x) = d(x k+1 , x) -d x1...x k (x k+1 , x),
where we use the approximated distance dx1...x k to give an estimate of d(x k+1 , x) using the first available k landmarks. This new representation removes the redundancy that exists between the different distance maps. • We code each map A i using a vector quantization scheme [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]. In our tests, we use codeblocks of size 3 × 3, and we quantize the resulting vectors of R 9 using a codebook of size 256. This results in a memory gain of 36 : 1 with respect to storing single precision floats. The constraints we have on the computation from compressed data are similar to the process of rendering from compressed textures [START_REF] Beers | Rendering from compressed textures[END_REF]. We use the differential coding strategy to cope with the particular redundancy of distances functions. The compression produces neither noticeable artifacts nor distortion of the extracted geodesics. The computational time overhead due to decompression is about 15% in CPU time.

Applications of Landmark-based Heuristics

Our landmark-based propagation can be applied verbatim to higher dimensional propagation and to triangulated surfaces in order to speed-up shortest path queries.

Constrained Path Planning Geodesics can be used to compute the path of a robot with various shape and motion constraints [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. Basically, each additional degree of freedom adds a new dimension to the domain in which the front propagation should be performed. Solving such high dimensional problems is time and memory consuming, so the use of a heuristic is highly desirable. In our experiment, the most important issue is the memory used by the full-grid classical Fast Marching, and the memory management strategy exposed in subsection 5.4 is crucial to scale to complex problems. In figure 18, one can see two examples of path extractions in 2D with one rotational additional degree of freedom. This results in 3D front propagation, and the corresponding speed function is depicted on the left. Figure 19 shows the influence of the heuristic strength λ on the cells explored by the front propagation. Geodesic Extraction on 3D Meshes The Fast Marching algorithm has been extended to 3D meshes in [START_REF] Sethian | Computing Geodesic Paths on Manifolds[END_REF]. Our heuristically driven front propagation extends to 3D surfaces in a straightforward manner. We use a constant metric g = 1 and as the surfaces considered do not have boundaries, the uniform seeding strategy (c) gives the best results for the E 1 metric. In figure 20, one can see the approximated distance d(x 1 , x) for various numbers of landmarks. In figure 1, one can see the front for various numbers of landmark points.

Conclusion

In this paper we have presented a simple modification of Fast Marching to speed-up the extraction of geodesics on images, higher dimensional data and triangulated surfaces. This modification uses a heuristic that drives the propagation and reduces the number of visited cells. This heuristic takes into account the remaining distance to the end point and two approaches are proposed to estimate this geodesic distance.

The first approach performs a backward propagation on a coarse grid and uses an interpolation scheme to retrieve the remaining distance on the full resolution grid. This method provides an important speed-up and numerical evidences show its efficiently as long as the coarse grid can resolve fine scale details of the metric. The other approach uses pre-computed distances to a set of landmark points that are combined during the propagation according to the triangle inequality. This method provides estimations with arbitrary accuracy when the number of landmarks is increased.

Figure 1 :

 1 Figure 1: In order to compute the geodesic path between two points on the surface, computation is needed on a large part of the surface when using classical fast marching (left). When the distance map to a set of landmarks is pre-computed, and the propagation is heuristically driven by these maps, only the colored region is explored, and it becomes smaller and smaller as the number of landmarks increases.

Figure 2 :

 2 Figure 2: An example of 2D path planning. The set of alive points according to increasing heuristic is shown in red.

Figure 3 :

 3 Figure 3: Resulting potential map g R for various values of R.

1Figure 4 :Figure 5 :

 45 Figure 4: Influence of heuristic strength and resolution on number of visited cells, Hausdorff error and computation time reduction.

Figure 6 :

 6 Figure 6: Graphical display of extracted geodesics for various heuristic strengths.• Computational time saving vs. heuristic strength λ: The saving is computed relative to the time spent by classical Fast Marching. In 2D the computation times decrease roughly linearly with the strength of the heuristic.• Computational time saving vs. heuristic resolution R (not shown):There is a constant overhead due to the coarse resolution computation (which results in an offset between the curves for R = 50% and R = 20%).For R = 20%, this overhead is balanced by the heuristic saving as soon as λ 5%. These tests clearly show that our algorithm can bring a large computational speed-up, but the parameters should be finely tuned to adapt to the characteristic of each map. For instance, these experiments show that the user must have some prior knowledge about the typical width of the tubular structures he wants to extract, and set the resolution R so that the coarse map V R still contains these structures.

Figure 7 :

 7 Figure 7: Extraction of geodesics in 3D.

Figure 8 :

 8 Figure 8: Cutting the square domain to compute circular paths.

z

  |d(x, z) -d(z, y)| . This equality is still valid in the continuous framework of Fast Marching, and in order to derive a useful heuristic, one can chose a small set z 1 , . . . , z n of Landmark points. The set of distance maps d k (x) = d(z k , x) is pre-computed, and we define the approximation dz1...zn (x, y) = sup k=1...n |d k (x) -d k (y)| .

Figure 12 :

 12 Figure 12: Example of propagation using a landmark-based heuristic with random seeding of base points.

Figure 13 :

 13 Figure 13: Distortion caused by our heuristic propagation on the extracted geodesics. The error is reported using the Hausdorff distance expressed in pixels (the size of the image is 256 × 256).

1 xFigure 14 :

 114 Figure 14: Graphical display of the approximated distance d(x 1 , x) for various seeding strategies.

1 Figure 15 :

 115 Figure 15: Decreasing of the approximation error E 1 (z 1 , . . . , z n ) for various seeding strategies.

Figure 16 :

 16 Figure 16: Computation speed-up E 2 (z 1 , . . . , z n ) obtained with a various number of landmarks. The areas colored in red correspond to the explored region without heuristic.

Figure 17 :

 17 Figure 17: Data structures used for the propagation.

Figure 18 :

 18 Figure 18: Examples of constrained path planning.

Figure 19 :

 19 Figure 19: Explored area for constrained path planning.

Figure 20 :

 20 Figure 20: Approximation of the distance function d(x 1 , x) using d(x 1 , x) on a 3D mesh, with an increasing number of landmarks.

Acknowledgment. We would like to thank Anthony Yezzi for his comments and corrections on the manuscript.