Gabriel Peyré
email: gabriel.peyre@polytechnique.fr

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Heuristically Driven Front Propagation for Geodesic Paths Extraction

Keywords:

In this paper we present a simple modification of the Fast Marching algorithm to speed up the computation using a heuristic. This modification leads to an algorithm that is similar in spirit to the A * algorithm used in artificial intelligence. Using a heuristic allows to extract geodesics from a single source to a single goal very quickly and with a low memory requirement. Any application that needs to compute a lot of geodesic paths can gain benefits from our algorithm. The computational saving is even more important for 3D medical images with tubular structures and for higher dimensional data. 80% 60% 40%

Shortest Path: Continuous and Discrete Algorithms

A large class of problems can be formulated as the extraction of a shortest paths for a given discrete or continuous metric. The applications of shortest paths thus range from tubular structures extraction in 3D medical images [START_REF] Deschamps | Fast Extraction of Minimal Paths in 3D Images and Applications to Virtual Endoscopy[END_REF] to path finding in video games [START_REF] Stout | Smart moves: Intelligent path-finding[END_REF]. The ability to quickly extract a geodesic path is highly desirable. In this paper we propose a modification of front propagation methods. We introduce a heuristic that drives the computations, and greatly reduces the computation time the Fast Marching algorithm. In this section we briefly review some basic facts about the Fast Marching algorithm and geodesics extraction and we also present a common framework that will allow to introduce our new approach.

Fast Marching Algorithm

The classical Fast Marching algorithm is presented in [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], and a similar algorithm was also proposed in [START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF]. This algorithm performs a front propagation that can be used to extract geodesic paths. The minimal length properties of geodesic 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]. We will illustrate various applications of geodesic paths in section 3.

In the continuous setting, a geodesic curve minimizes the weighted length of the curve for a given metric. In R d , we are given a potential function P (x) > 0, and the weighted geodesic distance between two points x 0 , x 1 ∈ R d , is defined as

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

To compute the distance function U (x)

def.

= d(x 0 , x) with an accurate and fast algorithm, this minimization can be reformulated as follows. The level set curve C t def.

= {x \ U (x) = t} propagates following the evolution equation ∂Ct ∂t (x) =

1 P (x)
-→ n x and the function U satisfies the nonlinear Eikonal equation:

||∇U (x)|| = P (x).
(

) 2
The function F = 1/P > 0 can be interpreted as the propagation speed of the front C t . The Fast Marching algorithm on an orthogonal grid makes use of an upwind finite difference scheme to compute the value u of U at a given point x i,j of a grid (the equation is written in R 2 for simplicity), solving:

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 P (x i,j ) 2 .
(

) 3
This is a second order equation that is solved as detailed for example in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF]. An optimal ordering of the grid points is chosen so that the whole computation only takes O(N log(N )), where N is the number of points.

Front Propagation Methods for Shortest Path

We now give a common framework for front propagation algorithms, which includes the Fast Marching procedure [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], the Dijkstra algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF] to compute shortest paths on graphs, and our heuristical front propagation procedure. These methods can be described using the following formalism:

• Alive is the set of grid 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 implementation-dependant sub-functions • A way to update the value U (y) at a given Trial point y. In the Fast Marching and our Heuristical front propagation, U (y) is computed by solving equation 3.

• 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, P(x) is chosen to minimize the number of visited points. We will explain in section 2 how to actually construct a function P that makes use of an heuristic. -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 = x1 is reached, stop the algorithm.

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

Geodesic Extraction

For the applications we have in mind, the computation of the geodesic distance U to a point x 0 is only a tool needed before the extraction of a geodesic that links this point to another point x 1 . This curve can be computed by extracting the parametric curve C(t) that solves the back propagation equation:

dC dt = - --→ ∇U with C(0) = x 0 .
This gradient descent is a very 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.

Heuristic 0% 40% 100% 60% 80%

1
x 0 x Fig. 3. Path planning using a multiresolution heuristic.

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. We also give some insights about how to choose the parameters of our method.

Fast Marching 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 will try to advance the front toward the goal point x 1 , and not isotropically. In order to do so, assume that together with the the current weighted distance to the start point U (x) = d(x 0 , x), we are able to have an estimate of the weighted distance that remains to be marched V (x) ≈ d(x 1 , x). Then our heuristical front propagation algorithm follows the implementation of the table 1 with a priority map

P(x) = U (x) + λ V (x). ( 4
)
This introduces the first parameter of our algorithm: the weight λ we use for the heuristic, which typically range from 0% to 100%. The rational 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]. On figure 2, one can see the effect of various choices for the parameter λ. The value λ = 0 corresponds to the classical Fast Marching propagation, which result in a very large region of Alive points (shaded in gray). However, as we increase the value of λ toward 100%, the explored region shrinks around the geodesic path that links x 0 to x 1 . There is 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 section 2.4 that although these numerical errors can accumulate during the propagation, the Hausdorff error on the extracted geodesic remains low both in synthetic case studies and in real applications. • We do not have an immediate access to the remaining distance d(x, x 1 ), since it would involve performing another front propagation from x 1 . We explain in the next section how to overcome this problem. Our heuristic ordering strategy is a reminiscence of the A * algorithm [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF], which was first introduced to solve problems in artificial intelligence. The use of heuristics and fast exploration strategies (such as IDA * [START_REF] Korf | Depth-first iterative-deepening: an optimal admissible tree search[END_REF]) is crucial in the artificial intelligence field to avoid the complexity of exploring the whole configuration space. These technics are also heavily used for path finding, for example in video games [START_REF] Stout | Smart moves: Intelligent path-finding[END_REF]. We show in section 3 that a large class of applications in computer vision can also gain benefit from this paradigm.

Multiresolution Heuristic

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 P is of size n × n, the query of P(y) thus requires:

• The pre-computation of a coarse potential map P R of size (Rn) × (Rn). This is done by first a pre-filtering of P (to avoid aliasing of high frequencies) and then a cubic spline re-interpolation on a coarser grid. • The pre-computation of the approximate distance map V of size (Rn) × (Rn). This is done by performing a full Fast Marching on a coarse grid, using potential P R , and starting from point x 1 . • During the heuristical front propagation starting from point x 0 , when P(y) is queried, we interpolate with cubic splines the value of V on the coarse grid to retrieve a value on the original grid. There is clearly a tradeoff between choosing a low R to reduce the computation time, and a high R so that V (x) approximates well d(x, x 1 ). We show in section 2.4 some insight about the correct value for this parameter R.

The new algorithm we propose allows to use multiresolution computation for the extraction of geodesic curve. 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 [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 purpose. However the problem of segmentation is in some sense dual to the geodesic extraction, in which we are not interested in the level sets of the distance function.

Reducing Memory Usage

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 chose 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: 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. This pointer-based representation of the neighboring relation is very convenient to extract the geodesic with a gradient descent. There is some memory overhead due to the fact that we explicitly store pointer to neighbors, but the fact that our scheme can explore significantly less cells that the classical fast marching allows to save much more memory, as shown in next section. The computing time overhead due to the use of a hash table is about 40% in all our tests.

cell NULL (0,0) (1,0) (1,

Numerical Validation

A Matlab implementation of our algorithm, together with the scripts needed to reproduce the figures of this paper, is freely available on Matlab Central [START_REF] Peyré | Fast marching toolbox, available on Matlab Central[END_REF].

In order to estimate the precision of the results, we use the Hausdorff error between the found paths and the paths obtained by fast marching without heuristic. On 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:

• 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 a lot the resolution without too much impact on the accuracy of the Heuristic. In contrast, one can see that the satellite image suffers from too huge variation when the resolution parameter R becomes smaller than 60% and then again for 90%. This is due to strong topological change in the path, as depicted in figure 5. 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.

Applications

In this section, we show the versatility of our method by proposing various applications where the extraction of geodesics is a central issue. We explain why our algorithm can bring a speed improvement and allows us to use simpler methods.

Volumetric Geodesics Extraction

3D geodesic extraction is very useful in medical volumetric data analysis. It can be applied to perform tubular structures 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]. On figure 7 one can see the extraction 3D geodesics on synthetic data (top and middle rows) and on real medical data (bottom row) for R = 20%. The red surface shows the boundary of the explored regions of alive cells. The computation time gain (C.t. gain) is also indicated.

2D Path Planning

Direct application of AI techniques has been intensively used for path finding, and A * is the method of choice for video games [START_REF] Stout | Smart moves: Intelligent path-finding[END_REF]. For the Euclidean case, faster and more complex graph-based methods have been proposed, such as visibility graphs [START_REF] Rohnert | Shortest paths in the plane with convex polygonal obstacles[END_REF].

The Fast Marching method can be used to produce continuous paths with a sub pixel precision [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], even when a non Euclidean metric is used. This is very powerful since it allows the modeling of smooth obstacles (such as forest, enemy repulsion, etc). In [START_REF] Melchior | Consideration of obstacle danger level in path planning using a* and fast-marching optimisation: comparative study[END_REF], the authors compare the Fast Marching and the A * algorithms for path planning. However they do not go one step further and connect these two powerful methods. On figure 3 one can see various paths extracted from a 2D map for a resolution R = 20% and a varying λ. Although the computed path can deviate from the real geodesic, our algorithm always produces locally acceptable path (smooth and avoiding obstacles).

Constrained Path Planning

Geodesics can be used to compute the path of a robot with various shape and motion constraints [START_REF] Kimmel | Optimal algorithm for shape from shading and path planning[END_REF]. Basically, each additional degree of freedom add 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 2.3 is crucial to scale to complex problems. The resolution R of the heuristic should be chosen carefully as a function of the typical width of the corridors and rooms the robot must pass through. In our tests we set R = 20% and we use a varying heuristic strength λ. Note however that although our algorithm can produce a wrong path for agressive heuristic usage, it never produces non-admissible moves.

On figure 8, 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 9 shows the influence of the heuristic strength λ on the cells explored by the front propagation.

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 about 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 and then the 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

P (x) = 1 ||C -x|| 1 1 + ||∇I(x)|| 2 + ε, (5)
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-andbound paradigm, which is a dichotomic search that avoids computing the closed geodesic for each point S on the segment CD. However, with our 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 10, we have shown a globally optimal circular geodesic, computed with various heuristic strength λ.

Heuristic 0% 40% 100% 60% 80% Fig. 10. 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 the 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 use a method derived from harmonic mesh parameterization [START_REF] Floater | Parameterization of Manifold Triangulations[END_REF]. This involves the resolution of a sparse linear system that search 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 once for all 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 disk-like charts that covers the region of alive vertices. On figure 1 and 11, one can see the algorithm in action on various meshes, and for various values of the parameter λ.

Coarse mesh Heuristic 0% 50% 100% 90% Fig. 11. Heuristically driven front propagation on 3D meshes.

Conclusion

In this paper we have presented a simple modification of the Fast Marching to speed up the computation time and the memory requirement of the algorithm. This modification is straightforward to implement, since it only involves the computation of a multiresolution heuristic to propagate the front toward the correct direction. Numerical tests on synthetic and real data show that this modification does not result in large distortion of the extracted curve. We examine some potential applications of this algorithm and show that it fit nicely into various existing computational frameworks.

Fig. 1 .

 1 Fig. 1. Heuristically-driven front propagation on a 3D mesh shown with increasing values of heuristic proportion. The various colors indicate the level sets of the distance function. Only the colored region is explored by our front propagation algorithm.

1 xFig. 2 .

 12 Fig. 2. An example of 2D path planning. The set of alive points according to increasing heuristic is shown in gray.

•

  Alive set: the starting point x0; • Trial set: the neighbors of x0; • 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:

  far, open or close cell* neighbors[START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF]; }; // pointers to the 4 neighbors

Fig. 4 .Fig. 5 .

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

•Fig. 6 .

 6 Fig. 6. Graphical display of extracted geodesics for various heuristic strengths.

Fig. 7 .

 7 Fig. 7. Extraction of geodesics in 3D.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Examples of constrained path planning.