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

Laurent Cohen
email: cohen@ceremade.dauphine.fr

Geodesic Methods for Shape and Surface Processing

This paper reviews both the theory and practice of the numerical computation of geodesic distances on Riemannian manifolds. The notion of Riemannian manifold allows to define a local metric (a symmetric positive tensor field) that encodes the information about the problem one wishes to solve. This takes into account a local isotropic cost (whether some point should be avoided or not) and a local anisotropy (which direction should be preferred). Using this local tensor field, the geodesic distance is used to solve many problems of practical interest such as segmentation using geodesic balls and Voronoi regions, sampling points at regular geodesic distance or meshing a domain with geodesic Delaunay triangles. The shortest path for this Riemannian distance, the so-called geodesics, are also important because they follow salient curvilinear structures in the domain. We show several applications of the numerical computation of geodesic distances and shortest paths to problems in surface and shape processing, in particular segmentation, sampling, meshing and comparison of shapes.

Manifold Geometry of Surfaces

In [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF], it was shown that finding the weighted distance and geodesic paths to a point leads to fast algorithms for image segmentation. In this chapter, we give a more general framework that is illustrated by different important applications.

This section introduces some basic definitions about local metric (a tensor field) on a Riemannian manifold and the associated notion of geodesic distance and minimal paths. The important point is that the geodesic distance to a set of starting points satisfies a non-linear differential equation, the Eikonal equation, which is solved to compute numerically the geodesic distance.

Riemannian Manifold

Parametric surface. A parameterized surface embedded in Euclidean space

M ⊂ R k is a mapping u ∈ D ⊂ R 2 → ϕ(u) ∈ M.
This definition can be extended to include surfaces not topologically equivalent to a disk, by considering a set of charts {D i } i that overlap in a smooth manner.

A curve is defined in parameter domain as a 1D mapping t ∈ [0, 1] → γ(t) ∈ D. This curve can be traced over the surface and its geometric realization is γ(t) def.

= ϕ(γ(t)) ∈ M. The computation of the length of γ in ambient k-dimensional space R k follows the usual definition, but to do the computation over the parametric domain, one needs to use a local metric (the first fundamental form) defined as follow.

Definition 1 (First fundamental form). For a parametric surface ϕ, one defines

I ϕ = ∂ϕ ∂u i , ∂ϕ ∂u j i,j=1,2 .
This local metric I ϕ defines at each point the infinitesimal length of a curve

L(γ) def. = 1 0 ||γ ′ (t)||dt = 1 0 γ ′ (t) T I ϕ (γ(t))γ ′ (t)dt.
This fundamental form is an intrinsic invariant that does not depend on how the surface is isometrically embedded in space (since the lengths depend only on this tensor field I ϕ ). In contrast, higher order differential quantities such as curvature might depend on the bending of the surface and are thus usually not intrinsic (with the notable exception of invariants such as the gaussian curvature).

Riemannian manifold. A parameterized surface is embedded into some Euclidean domain R k , which allows to define a local metric thanks to the first fundamental form I ϕ . It is however possible to consider directly a field of positive definite tensors on a parametric domain D = R s (in practice here s = 2 for surfaces or s = 3 for volumes). With a slight abuse in notations, we assimilate the resulting abstract surface M with D. Once again, we consider only surfaces globally parameterized by some Euclidean domain D and handling generic surfaces requires to split the manifold into overlapping charts.

Definition 2 (Riemannian manifold). A Riemannian manifold is an abstract parametric space M ⊂ R s equipped with a metric x ∈ M → H(x) ∈ R s×s positive definite.

Using the Riemannian metric, one can compute the length of a piecewise smooth curve γ :

[0, 1] → M L(γ) def. = 1 0 γ ′ (t)
T H(γ(t))γ ′ (t)dt.

At each location x, the Riemannian tensor can be diagonalized as follow H(x) = λ 1 (x)e 1 (x)e 1 (x) T + λ 2 (x)e 2 (x)e 2 (x)

T with 0 λ 1 λ 2 , (1)
and e 1 , e 2 are two orthogonal eigenvector fields. In fact, e i should be understood as direction (un-oriented) field since both e i and -e i are eigenvectors of the tensor. A curve γ passing at location γ(t) = x with speed γ ′ (t) has a shorter local length if γ ′ (t) is colinear to e 1 (x) rather than any another direction. Hence shortest paths (to be defined in the next section) tend to be tangent to the direction field e 1 .

In practice, the Riemannian metric H is given by the problem one wishes to solve. In image processing, the manifold is the image domain M = [0, 1] 2 equipped with a metric derived from the image (for instance its gradient). Figure 1 shows some some frequently used geodesic metric spaces:

Euclidean space: M = R s and H(x) = Id s . 2D shape: M ⊂ R 2 and H(x) = Id 2 .
Isotropic metric: H(x) = W (x)Id s , W (x) > 0 being some weight function. Parametric surface: H(x) = I ϕ (x) is the first fundamental form. Image processing: given an image I : [0, 1] 2 → R, one can use an edge-stopping weight W (x) = (ε+||∇ x I||) -1 . This way, geodesic curves can be used to perform segmentation since they will not cross boundaries of the objects. DTI imaging: M = [0, 1] 3 , and H(x) is a field of diffusion tensors acquired during a scanning experiment. For DTI imaging, the direction field e 1 indicates the direction of elongated fibers of the white matter (see [START_REF] Prados | Control theory and fast marching methods for brain connectivity mapping[END_REF]). The anisotropy of a metric H(x) is defined as

α(x) = λ 1 -λ 2 λ 1 + λ 2 = 2 √ ab -c 2 a + b ∈ [0, 1], for H(x) = a c c b . (2)
A metric with α(x) close to 1 is highly directional near x, whereas a metric with α(x) = 1 is locally isotropic near x.

Geodesic Distances

The local Riemannian metric H(x) allows to define a global metric on the space M using shortest paths. This corresponds to the notion of geodesic curves.

Definition 3 (Geodesic distance). Given some Riemannian space (M, H) with M ⊂ R s , the geodesic distance is defined as

∀ (x, y) ∈ M 2 , d M (x, y) def.

= min

γ∈P(x,y)

L(γ)

where P(x, y) denotes the set of piecewise smooth curves joining x and y

P(x, y) def. = {γ \ γ(0) = x and γ(1) = y} .
The shortest path between two points according to the Riemannian metric is called a geodesic. If the metric H is well chosen, then geodesic curves can be used to follow salient features on images and surfaces.

Definition 4 (Geodesic curve). A geodesic curve γ

∈ P(x, y) is such such that L(γ) = d M (x, y).
A geodesic curve between two points might not be unique, think for instance about two anti-podal points on a sphere. In order to perform the numerical computation of geodesic distances, we fix a set of starting points S = (x k ) k ⊂ M and consider only distance and geodesic curves from this set of points.

Definition 5 (Distance map). The distance map to a set of starting points

S = (x k ) k ⊂ M is defined as ∀ x ∈ M, U S (x) def. = min k d(x, x k ).
The main theorem that characterizes the geodesic distance is the following, that replaces the optimization problem of finding the minimum distance by a non-linear partial differential equation.

Theorem 1 (Eikonal equation).

If the metric H is continuous, then for any S ⊂ M, the map U S is the unique viscosity solution of the Hamilton-Jacobi equation

||∇ x U S || H(x) -1 = 1 with ∀ k, U S (x k ) = 0, ( 3
)
where ||v|| A = √ v T Av.
It is important to notice that, even if the metric x → H(x) is a smooth function, the distance function U S might not be smooth (it exhibit gradient discontinuities). This is why the machinery of viscosity solution is needed to give a sense to the solution of the Hamilton-Jacobi equation. See for instance [START_REF] Crandall | User's guide to viscosity solutions of second order partial differential equations[END_REF] for an introduction to viscosity solutions.

Once the distance map U S has been computed by solving the Eikonal equation (3), one can extract a geodesic joining any point x to its closest point x k ∈ S using a gradient descent on the function U S .

Theorem 2 (Gradient descent). The geodesic curve γ between x and its closest point in S solves

γ ′ (t) = - H(γ(t)) -1 ∇ γ(t) U S ||H(γ(t)) -1 ∇ γ(t) U S || with γ(0) = x.
The geodesic curve γ extracted using this gradient descent is parameterized with unit speed since ||γ ′ || = 1, so that γ : [0, T ] → M where T = d M (x, x k ).

Figure 2 shows examples of geodesic curves computed from a single starting point S = {x 1 } in the center of the image M = [0, 1] 2 and a set of points on the boundary of M. The geodesics are computed for a metric H(x) whose anisotropy α(x) (defined in equation ( 2)) is decreasing, thus making the Riemannian space progressively closer to the Euclidean space.

Image f α = 1 α = .5 α = .1 α = 0 Fig. 2. Examples of geodesics for a tensor metric with an decreasing anisotropy α (see equation [START_REF] Prados | Control theory and fast marching methods for brain connectivity mapping[END_REF] for a definition of this parameter). The tensor field H(x) is computed from the structure tensor of f as defined in equation [START_REF] Bronstein | Numerical Geometry of Non-Rigid Shapes[END_REF], its eigenvalues fields λi(x) are then modified to impose the anisotropy α.

For the particular case of an isotropic metric H(x) = W (x) 2 Id 2 , the geodesic distance and the shortest path satisfies

||∇ x U S || = W (x) and γ ′ (t) = - ∇ x U S ||∇ x U S || . ( 4
)
This corresponds to the Eikonal equation, that has been used to compute minimal paths weighted by W [START_REF] Cohen | Global Minimum for Active Contour models: A Minimal Path Approach[END_REF].

Numerical Computations of Geodesic Distances

In order to make all the previous definitions effective in practical situations, one needs a fast algorithm to compute the geodesic distance map U S . This section details Fast Marching algorithms based on front propagation that enable to compute the distance map by propagating the distance information from the starting points in S.

The basic Fast Marching algorithm and several extensions are exposed in the book on Fast Marching methods [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF]. For other applications to computer graphics and image processing one can see [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] and [START_REF] Kimmel | Numerical Geometry of Images: Theory, Algorithms, and Applications[END_REF]. The recent book [START_REF] Bronstein | Numerical Geometry of Non-Rigid Shapes[END_REF] treats all the details of the geometry of non-rigid surfaces, including geodesic distance computation and shape comparison. One can also see the two books [START_REF] Osher | Geometric Level Set Methods in Imaging, Vision, and Graphics[END_REF][START_REF] Paragios | Handbook of Mathematical Models in Computer Vision[END_REF] that contain review articles with some applications of Fast Marching and geodesic methods and in particular [START_REF] Cohen | Minimal paths and fast marching methods for image analysis[END_REF].

Front Propagation Algorithms

Depending on the properties of the metric, one needs to consider several algorithms, that all rely on the idea of front propagation. This family of algorithms allows to sort the computations in such a way that each point of the discretization grid is visited only once. This ordering is feasible for distance computation because the distance value of a grid point only depends (and can be computed) from a small number of points having only smaller distances. If one can sort the grid points with increasing distance, then one gets a coherent ordering of the computations. Of course, this is not that easy since this distance ordering would require the knowledge of the solution of the problem (the distance itself). But depending on the application, it is possible to devise a selection rule that actually select at each step the correct grid point.

A front propagation labels the points of the grid according to a state S(x) ∈ {Computed, F ront, F ar}.

During the iterations of the algorithm, a point can change of label according to

F ar → F ront → Computed.
Computed points S(x) = Computed are those that the algorithm will not consider any more (the computation of U S (x) is done for these points). Front points S(x) = F ront are the points being processed (the value of U (x) ≈ U S (x) is well defined but might change in future iterations). Far points S(x) = F ar are points that have not been processed yet.

In practice, a front propagation algorithm requires three key ingredients: Given a point x in the grid, a local set of neighbors N eigh(x) connected to x. A priority P(x) among points x in the front, that allows to select the point to process at a given iteration. In most application, this priority is computed as the current value of the distance P(x) def.

= U (x). Section 3.2 shows how to change this priority in order to speed up computations. A procedure x → U pdate(x) ∈ R that computes the distance value U (x) approximating U S (x) knowing the value U (x) for computed point and an approximate value for points in the front. This procedure usually solves some kind of equation that discretizes the Eikonal equation (3) one wishes to solve.

Listing 1 gives the details of the front propagation algorithm that computes a distance map U approximating U S (x) on a discrete grid. The following section details for actual implementations of the U pdate procedure for different metrics.

The numerical complexity of this scheme is O(n log(n)) for a discrete set of n points. This is because all the points are visited (tagged Computed ) once, and the selection of min P from the front points takes at most log(n) operations with a special heap data structure (although in practice it takes much less and the algorithm is nearly linear in time). -If S(y) = F ront, then U (y) ← min(U (y), U pdate(y)).

1. Initialization: ∀ x ∈ S, U (S) ← 0, S(x) ← F ront, ∀ y / ∈ S, S(y) ← F ar. 2. Select point: x ←-argmin S(z)=F ront P(z).
-Recompute the priority P(y).

Stop:

If x = x1, go back to 2.
Table 1: Front propagation algorithm.

Eikonal Equation Discretization

On a square grid. The classical Fast Marching algorithm, introduced by Sethian [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF], is a fast procedure to solve the Eikonal equation ( 3) for an isotropic metric H(x) = W (x) 2 Id s for a uniform regular grid that discretizes [0, 1] s . We recall this procedure for a planar domain s = 2 although it can be extended to any dimension.

In order to capture the viscosity solution of an Hamilton Jabobi equation, one cannot use standard finite differences because of the apparition of shocks and singularities in the solution of the equation. One needs to choose, at each grid point, the optimal finite difference scheme (differentiation on the left or on the right to approximate d/dx for instance). This optimal differentiation should be chosen in the direction where the solution of the equation decreases. This is called an upwind finite difference scheme, and on a 2D grid with spacing h it leads to find u = U pdate(x) at a grid point

x = x i,j that is the smallest 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 W (x i,j ) 2 . ( 5
)
The smallest solution of this equation leads to a stable and convergent scheme that can be used in the front propagation algorithm listing 1.

On a triangulation. The classical Fast Marching algorithm is restricted to isotropic metrics on a regular grid. This setting is useful for image and volumetric data processing, but in order to deal with arbitrary Riemmanian surfaces embedded in R k , one needs to modify equation [START_REF] Sethian | Level Sets Methods and Fast Marching Methods[END_REF].

Kimmel and Sethian [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF] have developed a version of the Fast Marching algorithm for a surface M ⊂ R k with metric W (x) for x in embedding space R k . In the continuous setting, a parametric surface (M, ϕ) embedded with a metric W (x) in ambient space corresponds to a Riemannian manifold with a metric

I ϕ (x)W (ϕ(x)) in parameter space x ∈ R 2 .
The algorithm of Kimmel and Sethian works on a triangulated mesh and treats the triangles of this mesh as locally flat and equipped with an isotropic metric W (x) 2 . The same algorithm can be used to process an anisotropic metric H(x) ∈ R 2×2 defined on a square lattice (an image), by locally connecting a pixel x to its 4 direct neighbors in order to create 4 adjacent triangles (that are flat). In order to describe the algorithm for these two settings (curved triangulated surface embedded in R k and Riemannian manifold with arbitrary metric H(x) for x ∈ R 2 ), we solve the Eikonal equation

||∇ x U || H -1 (x) = W (x)
locally on the triangle faces f ∈ F x adjacent to x in order to compute U pdate(x) ≈ U (x).

In order to compute the update value at a given vertex x, the algorithm computes an update value U pdate f (x) for each triangle f ∈ F x in the face 1-ring around x, F x = {f 1 , . . . , f k }. The resulting Fast Marching update step is defined as

U pdate(x) = min f ∈Fx U pdate f (x).
In order to derive the expression for U pdate f (x), one considers a planar triangle f = (x, x 1 , x 2 ) and denotes X = (x 1 -x, x 2 -x) ∈ R 2×2 . The known distances are u = (U (x 1 ), U (x 2 ))

T ∈ R 2 and one wishes to solve for U pdate f (x) = p = U (x).

The linear interpolation of U S can be written at the point x 1 and x 2 as

for i ∈ {1, 2}, U S (x i ) ≈ g, x i -x + p where g ≈ ∇ x U S .
With this approximation, equation (3) leads to a quadratic equation

U = X T g + pI ||g|| 2 H -1 (x) = W (x) 2 =⇒ I T QIp 2 + 2(I T Qu)p + (u T Qu -W (x) 2 ) = 0.
where

I = (1, 1) T ∈ R 2 and Q = (XH(x) -1 X T ) -1 ∈ R 2×2 .
The only admissible solution to this problem is

U pdate f (x) = p = I T Qu + (I T Qd) 2 + I T QI(u T Qu -W (x) 2 ) I T QI
There is some technical difficulties with this scheme on triangulations that contain obtuse angles or with metric H(x) with a large anisotropy, because the update procedure might not be monotone anymore. More accurate monotone schemes have been developped, see for instance [START_REF] Spira | An efficient solution to the eikonal equation[END_REF][START_REF] Bronstein | Weighted distance maps computation on parametric three-dimensional manifolds[END_REF][START_REF] Prados | Control theory and fast marching methods for brain connectivity mapping[END_REF]. We shall ignore these difficulties here and focus on the application of the numerical computation of geodesic distances.

Examples of propagations

2D isotropic propagation on a square grid. Figure 3 shows some examples of front propagation with the Fast Marching, for an isotropic

H(x) = W (x) 2 Id 2 .
The colored area shows, at some given step of the algorithm, the set of computed points (its boundary being the set of front points). During the iterations, the front propagates outwards until all the grid points are visited. The numerical complexity of this scheme is O(n log(n)) for a grid of n points.

Fig. 3. Examples of isotropic front propagation. The colormap indicates the values of the distance functions at a given iteration of the algorithm.

On rows 1 and 2, the potential W is computed using W (x) = f (x) so that geodesics tend to follow bright regions. On row 3, the potential W is computed using [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF] where c is chosen to match the intensity of the road to extract.

Figure 4 shows examples of distance functions to a starting point x 0 with the corresponding geodesics γ(t) extracted from some ending point x 1 . The front propagation is stopped when S(x 1 ) = Computed to avoid performing useless computations. The idea of using geodesics in order to extract salient curves in images as been introduced in [START_REF] Kimmel | Finding shortest paths on surfaces using level sets propagation[END_REF]. In practice, the difficult task is to design a metric W in order have meaningful geodesics. Here are some examples of possible choices, for image processing with an input image f :

Pixel value based potential: in many applications, one simply wishes to extract curves with a constant value c. In this case, one can use a potential like

W (x) = 1 ε + |f (x) -c| . (6)
Figure 4, left and middle, shows examples of such curves extractions. Also in many applications related to segmentation of tubular shapes, like vessels, we are looking for curves that are located in brighter or darker regions. In this case the potential can be chosen respectively as

W (x) = f (x) or W (x) = -f (x) (7)
and W should be rescaled to fill the range [ε, 1].

Gradient-based potential: for application such as edge detection one would like the geodesics to follow regions with high gradients. One can choose a potential such as

W (x) = ε + G σ * ||∇ x f ||,
where G σ is a smoothing kernel.

3D isotropic propagation on a square grid. The Fast Marching works the same way in any spacial dimension k and in particular can be used to extract shortest paths in 3D volumetric medical data. Such a volume is a discretization of a mapping f : [0, 1] 3 → R. Figure 5 shows a 3D display with a semi-transparent mapping that removes more or less parts of the data. The transparency at point (x, y, z) is defined as ρ(f (x, y, z)) where ρ : [f min , f max ] → [0, 1] is the α-mapping.

Figure 6 shows in red the front of the Fast Marching propagation, displayed as an isosurface of U S .

(x) = (|f (x) -f (x 0 )| + ε) -1
where x 0 is a point given by the user and supposed to lie inside some vessel. A geodesic follows nicely a vessel since its density is constant and thus the value of f is approximately equal to f (x 0 ) inside the vessel. Figure 7 shows other application of shortest path 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]. Isotropic propagation on a triangulated mesh. Figure 8 shows an example of propagation on a triangulated surface. The colored region corresponds to the points that are computed (its boundary being the front).

The propagation can be started from several starting points S = (x k ) k in order to compute the geodesic distance map U S . Figure 9 shows examples of such distances to several points together with geodesics. A geodesic γ links a point x to its closest point in S.

Anisotropic propagation on a square grid. In order to better follow the salient structures of an image f , one can replace the isotropic metric H(x) = W (x) 2 Id s (examples are given here in s = 2 dimensions) by a fully anisotropic metric H(x) ∈ R 2×2 which is a symmetric tensor field. This field might be given by the physical problem, such as the tensor field of DTI imaging [START_REF] Prados | Control theory and fast marching methods for brain connectivity mapping[END_REF]. Another option is to infer this field from some input image f .

The local orientation of a feature around a pixel x is given by the vector orthogonal to the gradient v(x) = (∇ x f ) ⊥ , which is computed numerically with finite differences (using maybe some little smoothing to cancel noise). This local direction information can be stored in a rank-1 tensor T 0 (x) = v(x)v(x)

T . In order to evaluate the local anisotropy of the image, one needs to average this tensor where the 4 entries of the tensor are smoothed against a gaussian kernel G σ of width σ > 0. The metric H corresponds to the so-called structure tensor, see for instance [START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF]. This local tensor T is able to extract both the local direction of edges and the local direction of textural patterns (see figure 11, left). Another option, that we do not pursue here, is to use the square of the Hessian matrix of f instead of the structure tensor.

T (x) = T 0 * G σ (x) (8)
In order to turn the structure tensor into a Riemannian metric, one can apply non-linear mapping to the eigenvalues, where ψ i is a decreasing function, for instance ψ i (x) = (ε + |x|) -1 for a small value of ε.

Fig. 10. Examples of anisotropic front propagation (from 9 starting points). The colormap indicates the values of the distance functions at a given iteration of the algorithm.

The metric is computed using the structure tensor, equation (8), of the texture f shown in the background.

Figure 10 shows an example of Fast Marching propagation using an anisotropic metric H(x). The front propagates faster in the direction of the main eigenvector field e 1 (x). Figure 11 shows distance map for a tensor field H(x) whose anisotropy α is progressively decreased, so that the geodesic distance becomes progressively Euclidean. 3 Applications and Extensions of Geodesic Distances

Shape Analysis

In order to analyze the shape of planar objects, one can consider the metric space obtained by restricting the plane to the inside of a planar domain.

Definition 6 (2D shape). A 2D shape S is a connected, closed compact set S ⊂ R 2 , with a piecewise-smooth boundary ∂S.

The geodesic distance inside such a shape is obtained by constraining the curve to lie inside S.

Definition 7 (Geodesic distance in S). The geodesic distance in S for the uniform metric is

d S (x, y) def. = min γ∈P(x,y) L(γ)
where

L(γ) def. = 1 0 |γ ′ (t)|dt.
where P(x, y) ⊂ S are the paths with starting point x and ending point y.

Figure 12 shows examples of shapes together with the geodesic distance to a starting point. The geodesic curve is the union of segments inside S and pieces of the boundary ∂S. The geodesic distance can be used to define several functions on the 2D shape. This section studies the eccentricity of a shape, as introduced by [START_REF] Ion | Shape matching using the geodesic eccentricity transform[END_REF] to perform shape recognition. These eccentric points are in fact located along the boundary.

Definition 8 (Eccentricity). The eccentricity E

S : M → R is E S (x)

Theorem 3 (Location of eccentric points). One has E(S) ⊂ ∂S.

A more general definition of eccentricity allows to replace the maximum by a weighted average of geodesic distances.

Definition 10 (α-eccentricity). The α eccentricity of some shape S is defined as

E α S (x) def. = S d S (x, y) α dy 1/α .
This eccentricity allows to generalize the notion of gravity center to the geodesic setting.

Definition 11 (Euclidean gravity center). The Euclidean gravity center is

argmin x S ||x -y|| 2 dy. The α-eccentric center is argmin x E α S (x).
Remark 1. For α = 2, the eccentric center is called geodesic gravity center (and equivalent to the Euclidean center in the case of an uniform metric).

Having defined a function such as E S inside a shape S, one can collect information about the shape using the histogram of that function.

Definition 12 (Descriptors). The eccentricity histogram descriptor

h S ∈ R m of a shape is ∀ i = 1, . . . , m, h S (i) = 1 |S| # x ∈ S \ i -1 m E S (x) -min(E S ) max(E S ) -min(E S ) < i m .
In particular, one can compare shapes by measuring the distance between the histograms

δ(h, h) 2 def. = m i=1 (h(i) -h(i)) 2 .
These histograms are invariant if one modifies a shape isometrically. In the plane, geodesic isometry of shapes are not interesting since they are rotations and translations. One can however consider approximate isometries such as articulations, that are useful to model deformations of planar shapes, as defined in [START_REF] Ling | Using the inner-distance for classification of articulated shapes[END_REF].

Definition 13 (ε-articulated object). An articulated object S can be split as

S = m i=1 S i i =j J ij , (a disjoint union) with diam(J ij ) ε.
Definition 14 (Articulation). An articulation is a mapping between two articulated shapes S, S ′ such that

f : S → S ′ = m i=1 S ′ i i =j J ′ ij is rigid on S i → S ′ i .
The eccentricity is approximately invariant for shapes that are modified by articulation. Starting from a shape library {S 1 , . . . , S p }, one can use the shape signature h S to do shape retrieval using for instance a nearest neighbor classifier, as shown in table 2. Figure 15 shows examples of typical shape retrievals. More complex signatures can be constructed out of geodesic distances and un-supervised recognition can also be considered. We refer to [START_REF] Ion | Shape matching using the geodesic eccentricity transform[END_REF] for a detailed study of the performance of shape recognition with eccentricity histograms. In a similar way, the eccentricity can be used to perform 3D surface retrieval, using the histograms displayed in figure 14.

Theorem 4 (Articulation and isometry). If f is an articulation, then

Heuristically Driven Propagation

The various implementations of the front propagation algorithm, pseudocode 1, use a simple priority P(x) = U x0 (x), where U (x) ≈ d M (x 0 , x) is the current value of the distance to the starting point. This strategy leads to an isotropic grow of the front which enforces the exploration of a large area of the computational grid. The advantage of using this priority is it does not favor any points and thus produces provable valid approximations of geodesic distance (both on a graph with Dijkstra and on a square/triangular grid with Fast Marching).

In order to reduce the computational burden, one could think about using more aggressive ordering of the front that favors some specific direction in the front. The hope is that the front would advance faster in the direction of the goal x 1 one wishes to reach. Ultimately, one would like the front to explore only points along the geodesic γ ∈ P(x 0 , x 1 ) joining the starting point to the ending point.

If one has an oracle: V (x) ≈ d(x 1 , x) that estimates the remaining geodesic distance from the current point x to the end x 1 , one can use as priority map

P(x) = U (x) + V (x).
The map V is called a heuristic since the exact distance d(x 1 , x) is not available in practice. The value of a good heuristic close to the real distance is revealed by the following theorem.

Theorem 5 (Geodesic segment). The function ψ(x) = d(x 0 , x) + d(x 1 , x) is minimal and constant ψ(x) = d(x 0 , x 1 ) along the geodesic path joining x 0 and x 1 .

In the setting of graph theory, the Dijkstra algorithm can be replaced by the A * (A-star), [START_REF] Nilsson | Problem-solving Methods in Artificial Intelligence[END_REF], which uses a heuristic to speed up computations. The following theorem proves the validity of this approach.

Theorem 6 (A * validity). If the heuristic satisfies V (x) d(x 1 , x), then the curve γ ∈ P(x 0 , x 1 ) extracted from the front propagation, algorithm 1, is a geodesic between x 0 and x 1 .

Over a continuous domain, one can invoke a similar (but weaker) theorem.

Theorem 7 (Explored area). If the heuristic satisfies V (x) d(x 1 , x), then the geodesic γ ∈ P 1 (x 0 , x 1 ) between x 0 and x 1 satisfies

{γ(t) \ t ∈ [0, 1]} ⊂ {x \ P(x) = U (x) + V (x) P(x 1 )} .
This theorem shows why it is important to estimate the geodesic distance by below, since otherwise the region explored by the algorithm might not contain the true geodesic. Figure 16 shows examples of heuristics that approximate the true remaining distance by bellow. One can see how the explored area of the propagation progressively shrinks while containing the true geodesic. Such a heuristic is however impossible to use in practice since one does not have direct access to the distance during the propagation. Many strategies can be used to estimate a heuristic. For instance, on a Riemannian metric (M, H(x)), one could use

V (x) = ρ||x -x 1 || where ρ = min x =0,v =0 ||v|| H(x) .
In this case, ρ is the minimum eigenvalue of all the tensors H(x). This heuristic estimates the geodesic distance with a Euclidean distance and satisfies V (x)

d(x 1 , x).
For a propagation on a graph (A * algorithm) that is embedded in Euclidean space according to i ∈ V → x i ∈ R k , one could also define

∀ i ∈ V, V (i) = ||x i1 -x i ||,
where i 1 is the index of the ending point. This heuristic also satisfies V (i) d(i 1 , i).

These Euclidean heuristics performs poorly on spaces that are not relatively flat. In order to compute more accurate heuristic, we use an expression of the geodesic distance as a minimization. This kind of approximation has been used first in graph theory [START_REF] Goldberg | Computing the shortest path: A* search meets graph theory[END_REF] and it is defined in a continuous setting in [START_REF] Peyré | Heuristically driven front propagation for fast geodesic extraction[END_REF]. This leads to a heuristic V (x) = d(x, x 1 ) that has the following properties. In a numerical application that requires the extraction of many geodesics in real time over a large domain, one can pre-compute (off-line) the set of distance maps to the landmarks {d(x, z i )} m i=1 . At run time, this set of distances is used to compute the heuristic and speed up the propagation. Figure 17 shows how the quality of the heuristic increases with the number of landmarks. Figure 18 shows an application to geodesic extraction on 3D meshes.

Surface Sampling

In order to acquire discrete samples from a continuous surface, or to reduce the number of samples of an already acquired mesh, it is important to be able to seed evenly a set of points on a surface. This is relevant in numerical analysis to have a good accuracy in computational simulations, or in computer graphics to display 3D models with a low number of polygons. In practice, one typically wants to enforce that the samples are approximately at the same distance from each other. The numerical computation of geodesic distances is thus a central tool, that we are going to use both to produce the sampling and to estimate the connectivity of a triangular mesh.

Farthest Point Sampling

A sampling of a Riemannian surface M is a set of points {x 1 , . . . , x n } ⊂ M. If the surface is parameterized by ϕ : [0, 1] 2 → M, the easiest way to compute a sampling is to seed points regularly over the parametric domain

∀ (i, j) ∈ {1, . . . , √ n} 2 , x i,j = ϕ(i/ √ n, j/ √ n).
This strategy performs poorly if the mapping ϕ introduces heavy geodesic distortion and the sampling might not be regular any more for the geodesic metric on the surface. In order to ensure the quality of a sampling, one can use the notion of a well separated covering.

Definition 15 (ε-covering). A sampling {x 1 , . . . , x n } ⊂ M is an ε-covering if i B ε (x i ) = M where B ε (x) def.
= {y \ d M (x, y) ε} .

Definition 16 (ε-separated). A sampling {x 1 , . . . , x n } ⊂ M is ε-separated if max(d M (x i , x j )) ε.
The farthest point sampling algorithm is a simple greedy strategy able to produce quickly a good sampling. This algorithm has been introduced in image processing to perform image approximation [START_REF] Eldar | The farthest point strategy for progressive image sampling[END_REF]. It is used in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF] together with geodesic Delaunay triangulation (to be defined in the next section) to do surface remeshing. The detection of saddle points (local maxima of the geodesic distance) is used in [START_REF] Cohen | Multiple Contour Finding and Perceptual Grouping Using Minimal Paths[END_REF] to perform perceptual grouping.

3 gives the details of this iterative algorithm. In particular, note that the update of the distance d(x) to the set of already seeded points goes faster at each iteration since the domain of update is smaller when the number of points increases.

, . . . , x n } is an ε-covering that is ε-separated for ε = max i=1,...,n min j=1,...,n d M (x i , x j ).
Note however that there is no simple control on the actual number of samples n required to achive a given accuracy ε. We refer to [START_REF] Clarkson | Building triangulations using epsilon-nets[END_REF] for an in-depth study of the approximation power of this greedy sampling scheme.

Figure 19 shows examples of farthest point sampling with a uniform (top row) and a spatially varying isotropic metric W (x) (bottom row). One can see that this scheme seeds more points in areas where the metric W is large. One can thus control the sampling density by modifying the metric W .

Triangulations

Having computed, for instance with farthest points, a sampling {x i } i∈V ⊂ M, the next step is to compute some connectivity between the samples in order to build a graph, or even better, a triangulation. The problem of surface remeshing has been studied extensively in computer graphics, see the survey [START_REF] Alliez | Recent advances in remeshing of surfaces[END_REF]. This section explains a solution based on the geodesic Delaunay triangulation.

The following definition generalizes the notion of an Euclidean Voronoi diagram, to an arbitrary surface.

Definition 17 (Voronoi segmentation). The Voronoi segmentation of a sampling {x

i } i∈V ⊂ M is M = i V i with (V i ) m i=1 def. = Voronoi M ({x i } i ) where V i def. = {x \ ∀ j = i, d M (x, x i ) d M (x, x j ).}
Each Voronoi cell V i is thus composed of points that are closer to x i than to any other sampling point. The boundary between two adjacent cells V i and V j is thus a piece of curve at equal distance between x i and x j . One can then compute the graph dual to a given partition, which joins together pair of adjacent cells. This leads to the notion of Delaunay graph.

Definition 18 (Geodesic Delaunay graph). The Delaunay graph

(V, E) of a sampling {x i } i∈V ⊂ M is defined for V = {1, . . . , n} as E = {(i, j) ∈ V \ ∂V i ∩ ∂V j = 0} .
The main interest of this Delaunay graph is that, if the number of points is large enough to capture the topology of the surface (for instance at least 4 points are needed on a sphere), then one gets a valid triangulation. This theorem means that one can find a set of faces F such that (V, E, F ) is a triangulated mesh. One can see [START_REF] Leibon | Delaunay triangulations and voronoi diagrams for riemannian manifolds[END_REF] for a theoretical study of geodesic Delaunay triangulations.

Examples of Meshing and Remeshing

This Delaunay triangulation can thus be used to perform a geodesic meshing or re-meshing of any Riemannian surface, as explained in [START_REF] Peyré | Geodesic remeshing using front propagation[END_REF].

Figure 20 shows examples of Voronoi segmentations on the plane for various isotropic Riemannian metrics W (x). The Delaunay graph allows to define a planar mesh of points evenly sampled according to the metric. Figure 21 shows examples of Voronoi cells on a surface embedded in R 3 .

Instead of using a constant or an isotropic metric W (x), one can use a fully anisotropic metric H(x) ∈ R 2×2 . The local dominant eigenvector e 1 (x) given in the decomposition (1) of the tensor gives the local preferred direction of the triangles and the anisotropy λ 1 (x)/λ 2 (x) describe how much the triangles should be stretched in this direction. Figure 22 shows an example of meshing with a metric of decreasing anisotropy. Figure 23 shows an anisotropic farthest point meshing with an increasing number of sampling points.

In order to mesh the interior of a planar shape S ⊂ R 2 , one can use the Euclidean metric inside the shape and compute a geodesic Delaunay triangulation. Some care should be made during the algorithm so that the boundary of the domain is included in the delaunay triangulation. This requires splitting boundary edges if they disappear from the Delaunay graph during the algorithm. Figure 24 shows some examples of shape meshing with this uniform metric. This triangulation is however very close to the usual definition of a planar Euclidean Delaunay triangulation. In contrast, one can use a non-uniform metric W (x) and compute a sampling inside the shape that conforms itself to this density. Figure 24 shows a sampling and meshing that uses a metric W (x) = (ε + d(x, ∂S)) -1 that tends to seed more points on the boundary of the shape S. Figure 25 shows an example of uniform remeshing of a 3D surface acquired from medical imaging with an increasing number of points. Figure 26 shows how one can adapt the density by defining a non-constant isotropic metric on the surface.

An option to compute this metric is to use a texture mapped on the surface. Starting from some parametric surface: ϕ : D ⊂ [0, 1] 2 → M, a texture T is a mapping T : [0, 1] 2 → R. It allows to define an isotropic metric using for instance an edge adaptive function

∀ x ∈ D, H(x) = ψ T (x)Id 2 .
where the edge-based stopping function is ψ T (x) = (||∇ x T || + ε) -1 . Figure 27 shows examples of remeshing with a texture-adapted metric with a decreasing value of ε (increasing adaptivity).

Conclusion

This chapter has reviewed several applications of Riemannian metrics in computer vision and graphics. In particular, the use of geodesic distances and shortest paths is useful in many areas of these fields. The design of adapted isotropic tensor field 41 points inserted 75 points inserted 150 points inserted or anisotropic metrics allows to solve efficiently segmentation, sampling, meshing and recognition problems with fast algorithms.

Semi-

Fig. 27. Adaptive remeshing with a density given by a texture.

Fig. 1 .

 1 Fig. 1. Examples of Riemannian metrics (top row) and geodesic distances and curves (bottom row). The blue/red colormap indicates the geodesic distance to the starting point. From left to right: euclidean (H(x) = Id2 restricted to M = [0, 1] 2 ), planar domain (H(x) = Id2 restricted to M = [0, 1] 2 ), isotropic (H(x) = W (x) 2 Id2 with W computed from the image using (6)), Riemannian manifold metric (H(x) is the structure tensor of the image, see equation (8)) and 3D surface (H(x) corresponds to the first fundamental form).

3 .

 3 Tag: S(x) ← Computed. 4. Update neighbors: for all y ∈ N eigh(x), -If S(y) = F ar, then S(y) ← F ront and U (y) ← U pdate(y).

Fig. 4 .

 4 Fig. 4. Example of distance functions (top row) and geodesics (bottom row).

Fig. 5 .

 5 Fig. 5. Example of semi-transparent display of volumetric data..

Figure 6

 6 Figure6shows some examples of geodesic extraction on a medical image that represents tubular structures (blood vessels) around the heart. The potential W (x) is chosen as W (x) = (|f (x) -f (x 0 )| + ε) -1 where x 0 is a point given by the user and supposed to lie inside some vessel. A geodesic follows nicely a vessel since its density is constant and thus the value of f is approximately equal to f (x 0 ) inside the vessel. Figure7shows other application of shortest path 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].

Fig. 6 .

 6 Fig. 6. Example of volumetric Fast Marching evolution (top row) and geodesic extractions (bottom row).

Fig. 7 .

 7 Fig. 7. Left: vessel extraction. Right: tubular structure extraction.

Fig. 8 .

 8 Fig. 8. Example of Fast Marching propagation on a triangulated mesh.

Fig. 9 .

 9 Fig. 9. Examples of geodesic extraction on a mesh with an increasing number of starting points.

T (x) = µ 1 e 1 e 1 T + µ 2 e 2 e 2 T

 12 =⇒ H(x) = ψ 1 (µ 1 )e 1 e 2 T + ψ 2 (µ 2 )e 2 e 2 T . (9)

Fig. 11 .

 11 Fig. 11. Left image: example of texture together with the structure tensor field, computed using equation (9). Right: examples of anisotropic distances (top row) and Voronoi diagrams (bottom row) with a decreasing anisotropy α (see equation (2) for a definition of this parameter).

Fig. 12 .

 12 Fig. 12. Geodesics inside a 2D shape.

  d S (x, y) = max y∈∂S d S (x, y).

Fig. 13 .

 13 Fig. 13. Example of eccentricity ES and corresponding histograms hS.

Figure 13 (

 13 Figure 13 (top row) shows several examples of eccentricity. The colormap indicates in blue points with small eccentricity. The points for which the minimum in the definition of E S is obtained are called eccentric. The set of eccentric point is denoted as E(S). Definition 9 (Eccentric points). An eccentric point x ∈ E(S) satisfies ∃y ∈ S, E S (y) = d(x, y).These eccentric points define regions of influence which perform a segmentation of the shape as follow S =

Fig. 14 .

 14 Fig. 14. Example of eccentricity and corresponding histograms for 3D surfaces.

  |d S (x, y) -d S ′ (x, y)| mε and |E S (x) -E S ′ (x)| mε.

1 .

 1 Dataset:shapes {S1, . . . , Sp} (binary images). 2. Preprocessing: compute eccentricity descriptors hS i . 3. Input: shape S. 4. Retrival: return i ⋆ = argmin i δ(hS, hS i ).

Fig. 15 .

 15 Fig. 15. Examples of shape recognitions. The shape on the left is the input S and the second shape in each raw is Si⋆ .

Fig. 16 .

 16 Fig.[START_REF] Kothe | Edge and junction detection with an improved structure tensor[END_REF]. Example of propagations with a priority P(x) = U (x) + λV (x) for λ = 0, 0.5, 0.9.

Theorem 8 (

 8 Reversed triangular inequality). For all (x, y) ∈ M, one has d(x, y) = sup z |d(x, z) -d(z, y)| . If one restricts the minimum to a small subset of landmark points {z 1 , . . . , z n } ⊂ M, one can define the following approximate distance dz1...zn (x, y) = sup k=1...n |d k (x) -d k (y)| ,

Fig. 17 .

 17 Fig. 17. Heuristically driven propagation in 2D with an increasing number of landmark points.

Fig. 18 .

 18 Fig. 18. Heuristically driven propagation on a 3D mesh with landmark points.

Fig. 19 .

 19 Fig. 19. Examples of farthest point sampling (the colormap indicates the distance function to the seeds).

1 .

 1 Initialization: x1 ←random, d(x) ← dM(x1, x), set i = 1. 2. Select point: xi+1 = argmax x d(x), ε = d(xi+1).

3 .

 3 Local update of the distance: d(x) ← min(d(x), dM(xi+1, x)). This update is restricted to the set of points {x \ dM(xi+1, x) < d(x)}. 4. Stop: If i < n or ε > ε0, set i ← i + 1 and go back to 2.

Fig. 20 .

 20 Fig. 20. Examples of sampling and triangulations with an isotropic metric H(x) = W (x) 2 Id2. The sampling is denser in the regions where the metric is small (dark).

Fig. 21 .

 21 Fig. 21. Example of Voronoi segmentations for an increasing number of seeding points.Theorem 11. For a large enough number of points, the Delaunay graph is a valid triangulation.

α 1 Fig. 22 .

 122 Fig. 22. Meshing of a square with a metric of decreasing anisotropy of a same synthetic tensor field. Top: Voronoi diagrams, tensor fields and points added by the algorithm (last image is the Euclidean case). Bottom: resulting meshes.

Fig. 23 .

 23 Fig. 23. Anisotropic meshing of a square with an increasing number of points.

Fig. 24 .Fig. 25 .

 2425 Fig. 24. Shape meshing with an increasing number of points. Left and center: uniform meshing, right: adaptive meshing.

Fig. 26 .

 26 Fig. 26. Adaptive remeshing with a linearly increasing density.

Table 2 :

 2 Shape retrieval process.

Table 3 :

 3 Farthest point sampling algorithm.Theorem 10 (Farthest seeding properties). The farthest point sampling {x 1

	The output sampling of the algorithm enjoys the property of being a well
	separated covering of the manifold.


Gabriel Peyré

Laurent Cohen

Geodesic Computations for Fast and Accurate Surface Remeshing and Parameterization

Keywords: Remeshing, geodesic computation, fast marching algorithm, mesh segmentation, surface parameterization, texture mapping, deformable models. 3D model Speed function Uniform Semi-adaptive Adaptive

In this paper, we propose fast and accurate algorithms to remesh and flatten a genus-0 triangulated manifold. These methods naturally fits into a framework for 3D geometry modeling and processing that uses only fast geodesic computations. These techniques are gathered and extended from classical areas such as image processing or statistical perceptual learning. Using the Fast Marching algorithm, we are able to recast these powerful tools in the language of mesh processing. Thanks to some classical geodesic-based building blocks, we are able to derive a flattening method that exhibit a conservation of local structures of the surface.

On large meshes (more than 500 000 vertices), our techniques speed up computation by over one order of magnitude in comparison to classical remeshing and parameterization methods. Our methods are easy to implement and do not need multilevel solvers to handle complex models that may contain poorly shaped triangles.

Introduction

The applications of 3D geometry processing abound nowadays. They range from finite element computation to computer graphics, including solving all kinds of surface reconstruction problems. The most common representation of 3D objects is the triangle mesh, and the need for fast algorithms to handle this kind of geometry is obvious. Classical 3D triangulated manifold processing methods have several well-identified shortcomings: mainly, their high complexity when dealing with large meshes, and their numerical instabilities.

To overcome these difficulties, we propose a geometry processing pipeline that relies on intrinsic information of the surface and not on its underlying triangulation. Borrowing from well-established ideas in different fields (including image processing and perceptual learning) we are able to process very large meshes efficiently.

Overview

In Section 2 we introduce some concepts we use in our geodesic computations. This includes basic facts about the Fast Marching algorithm, and a recently proposed greedy algorithm for manifold sampling.

To flatten each patch of a segmented surface, we will recall some recent advances in perceptual inference learning in Section 3. Combining these techniques with our geodesic computational framework will lead to an elegant solution to the flattening problem for large meshes.

In the conclusion, we will show the two algorithms in action, and see how we can texture large meshes faster than current techniques would otherwise allow. We will then give a complete study of the timings of each part of our algorithm, including a comparison with classical methods.

Related Work

Surface Remeshing and Finite Elements. Remeshing methods roughly fall into two categories:

• Isotropic remeshing: a surface density of points is defined, and the algorithm tries to position the new vertices to match this density. For example the algorithm of Terzopoulos and Vasilescu [START_REF] Terzopoulos | Adaptive Meshes and Shells: Irregular Triangulation, Discontinuities, and Hierarchical Subdivision[END_REF] uses dynamic models to perform the remeshing. Remeshing is also a basic task in the computer graphics community, and [START_REF] Surazhsky | Isotropic Remeshing of Surfaces: a Local Parameterization Approach[END_REF] have proposed a procedure based on local parameterization. • Anisotropic remeshing: the algorithm takes into account the principal directions of the surface to align locally the newly created triangles and/or rectangles. Finite element methods make heavy use of such remeshing algorithms [Kunert, 2002]. The algorithm proposed in [START_REF] Alliez | Anisotropic Polygonal Remeshing[END_REF] uses lines of curvature to build a quad-dominant mesh. The importance of using geodesic information to perform this remeshing task is emphasized in [START_REF] Sifri | Geodesic-based Surface Remeshing[END_REF].

Greedy solutions for sampling a manifold (see Section 2.2) have been used with success in other fields such as computer vision (component grouping, [Cohen, 2001]), halftoning (void-and-cluster, [Ulichney, 1993]) and remeshing (Delaunay refinement, [Ruppert, 1995]).

Flattening and Parameterization.

The flattening problem can be seen as a particular instance of parameterization. The work of [START_REF] Eck | Multiresolution Analysis of Arbitrary Meshes[END_REF] first introduces the harmonic formulation for the resolution of the mesh parameterization problem. Most of these classical methods come from graph-drawing theory, and [START_REF] Floater | Parameterization of Manifold Triangulations[END_REF] gives a survey of these techniques. Authors of [START_REF] Desbrun | Intrinsic Parameterizations of Surface Meshes[END_REF] give an in-depth study of the various energies that can be built to flatten a mesh. The flattening algorithm of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF] is based on methods for finding parameters that reduce a dataset's dimensionality. Such methods have been developed for the purpose of perceptual learning [Tenenbaum et al., 2000, Roweis and[START_REF] Roweis | [END_REF], and we will explain in Section 3 how to exploit these methods to handle a 3D mesh with a large amount of vertices.

Geodesic Remeshing

Fast Marching Algorithm

The classical Fast Marching algorithm is presented in [Sethian, 1999], and a similar algorithm was also proposed in [Tsitsiklis, 1995]. This algorithm is used intensively in computer vision, for instance it has been applied to solve global minimization problems for deformable models [START_REF] Cohen | Global Minimum for Active Contour Models: A Minimal Path Approach[END_REF].

This algorithm is formulated as follows. Suppose we are given a metric P (s)ds on some manifold S such that P > 0. If we have two points x 0 , x 1 ∈ S, the weighted geodesic distance between x 0 and x 1 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 geodesic 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 , where -→ n x is the exterior unit vector normal to the curve at 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:

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 .
This is a second-order equation that is solved as detailed for example in [Cohen, 2001]. 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.

In [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF]], a generalization to an arbitrary triangulation is proposed. This allows performing front propagations on a triangulated manifold, and computing geodesic distances with a fast and accurate algorithm. The only issue arises when the triangulation contains obtuse angles. The numerical scheme presented above is not monotone anymore, which can lead to numerical instabilities. To solve this problem, we follow [START_REF] Kimmel | Computing Geodesic Paths on Manifolds[END_REF]] who propose to "unfold" the triangles in a zone where we are sure that the update step will work. Figure 2 shows the calculation of a geodesic path computed using a gradient descent of the distance function.

A Greedy Algorithm for Uniformly Sampling a Manifold

A new method for sampling a 3D mesh was recently proposed in [START_REF] Peyré | Geodesic Remeshing Using Front Propagation[END_REF]] that follows a farthest point strategy based on the weighted distance obtained through Fast Marching on the initial triangulation. This is related to the method introduce in [Cohen, 2001]. A similar approach was proposed independently and simultaneously in [START_REF] Moenning | Fast Marching Farthest Point Sampling[END_REF]. It follows the farthest point strategy, introduced with success for image processing in [START_REF] Eldar | The Farthest Point Strategy for Progressive Image Sampling[END_REF] and related to the remeshing procedure of [Chew, 1993].

This approach iteratively adds new vertices based on the geodesic distance on the surface. Figure 3 shows the first steps of our algorithm on a square. The result of the algorithm gives a set of vertices uniformly distributed on the surface according to the geodesic distance.

Once we have found enough points, we can link them together to form a geodesic Delaunay triangulation. This is done incrementally during the algorithm, and leads to a powerful remeshing method. Figure 4 shows progressive remeshing of the bunny and the David. In order to have a valid triangulation, the sampling of the manifold must be dense enough (for example 100 points is not enough to capture the geometry of the ears of the bunny). A theoretical proof of the validity of geodesic Delaunay triangulation can be found in [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF], and more precise bound on the number of points is derived in [START_REF] Onishi | Estimation of the necessary number of points in Riemannian Voronoi diagram[END_REF]. Note that our algorithm works with manifolds with boundaries, of arbitrary genus, and with multiple connected components.

Adaptive Remeshing

In the algorithm presented in Sections 2.2, the fronts propagate at a constant speed which results in uniformly spaced mesh. To introduce some adaptivity in the sampling performed by this algorithm, we use a speed function F = 1/P (which is the right-hand side of the Eikonal equation ( 2)) that is not constant across the surface. Figure 5 shows the progressive sampling of a square using a speed function with two different values. The colors show the level sets of the distance function U to the set of points. The local density of vertices can also reflect some geometric properties of the surface. The most natural choice is to adapt the mesh in order to be finer in regions where the local curvature is larger. The evaluation of the curvature tensor is a vast topic. We used a robust construction proposed recently in [START_REF] Cohen-Steiner | Restricted Delaunay Triangulations and Normal Cycles[END_REF]. Let us denote by τ (x) def.

= |λ 1 |+ |λ 2 | the total curvature at a given point x of the surface, where λ i are the eigenvalues of the second fundamental form. We can introduce two speed functions F 1 (x)

def. = 1 + ετ (x) and F 2 (x) def. = 1 1+µτ (x)
, where ε and µ are two user-defined parameters. Figure 6 (a) shows that by using function F 1 , we avoid putting more vertices in regions of the surface with high curvature. The speed function F 1 can be interpreted as an "edge repulsive" function. On the other hand, function F 2 could be called "edge attractive" function, since it forces the sampling to put vertices in region with high curvature such as mesh corners and edges. Figure 6 (b) shows that this speed function leads to very good results for the remeshing of a surface with sharp features, which is obviously not the case for the "edge repulsive" speed function (Figure 6 (a)).

Fast Geodesic Parameterization

The flattening problem can be seen as a particular instance of the more generic problem of mesh parameterization. Given a genus-0 triangulated manifold S homeomorphic to a disc, it consists in finding a map f : S → U, where U ⊂ R 2 is a planar domain.

Geodesic Flattening and IsoMap

Recently, some nonlinear algorithms for dimensionality reduction have appeared in the community of perceptual manifold learning. The most notable are IsoMap [START_REF] Tenenbaum | A Global Geometric Framework for Nonlinear Dimensionality Reduction[END_REF] and Locally Linear Embedding (LLE) [START_REF] Roweis | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF]. Interestingly, the only echo of these techniques in the computer graphics community seems to be the multi-dimensional scaling approach to flattening of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF]. This method is closely related to IsoMap, and we will see that it shares its main drawbacks.

We start with a given set of points {x 1 , . . . , x n } on our manifold, and we seek f (x i ) = x i ∈ R 2 such that the mapping minimizes some measure of distortion. The most natural constraint is to try to keep the same distance between points, which is exactly what IsoMap is doing by requiring that d(x i , x j ) ≈ || x ix j ||, where d stands for (some approximation of) the geodesic distance on the manifold. The method of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF] is very close to this approach, since it uses the geodesic distance d computed via the Fast Marching algorithm presented in Section 2.1.

The major bottleneck of this method is that it needs to compute all pairwise distances d(x i , x j ). To overcome this difficulty, the authors of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF] proposed to restrict the computations to a small set of points, which gives rise to three questions:

• What should be done to speed up computation?

• How should we choose this small set of base points?

• How should we extend the map f from this small set of points to the rest of the mesh?

In the next subsection, we will show how the LLE algorithm can bring a important speed improvement that answers the first question. The answer to the two last questions will be given in Subsections 3.3 and 3.4 respectively, with an extension of LLE to triangulated manifolds.

Speeding Up Computation with LLE

The LLE algorithm is explained in detail in [START_REF] Roweis | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF]. The goal of the algorithm is to find a low-dimensional embedding in R d of a set of points {x 1 , . . . , x n } in R s , s > d. The only parameter of this algorithm is an integer K that measures the size of the neighborhood of each point. We will denote by N i the K-neighborhood of x i , that is to say N i def.

= {x m(1) , . . . , x m(K) }, where x m(j) is the j th closest point to x i for the Euclidean metric. We will briefly recall the two main steps of the procedure:

Step 1: First, for each point x i , we are looking for some weights w i,j that locally best reconstruct the manifold, from the set N i only, by minimizing

E 1 ({w i,j } j ) = ||x i - j w i,j x j || 2 .
(3)

We further enforce that w i,j = 0 if x j / ∈ N i , and that j w i,j = 1. This imposes that the reconstruction is both local and invariant under affine transformations. In a Euclidean setting, the minimization (3) requires the introduction of the Gram matrix C(x) defined by

(C(x)) i,j def. = x -x m(i) , x -x m(j) (4)
The solution of the minimization (3) is then

w i,j = 1 Ni (x j ) q (C(x i ) -1 ) k,q p,q (C(x i ) -1 ) p,q
, where x j = x m(k) .

The value of 1 Ni (x) is equal to 1 if x ∈ N i , and 0 otherwise.

Step 2: To reconstruct the manifold in low dimension (here in 2D), we want to solve a global minimization procedure, for

x i ∈ R 2 : minimize E 2 ({ x i }) = i || x i - j w i,j x j || 2 .
subject to x i = 0 and || x i || 2 = 1 to avoid a degenerate solution. To solve this problem, we need to form the matrix M def.

= (W -Id)

T (W -Id), where W is a sparse matrix containing all the weights. The eigenvector of M with lowest eigenvalue it the constant vector 1 which should be discarded. The d following eigenvectors give us the coordinates of our embedding in R d for each point.

The fact that we only need to perform computations on sparse matrices allows an improvement of one order of magnitude over dense procedures such as in [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF].

Geodesic LLE

In the previous section we saw the classical LLE algorithm in a Euclidean setting. To solve the flattening problem for a mesh, we need to extend these computations to the manifold setting. The following modifications allow such an extension.

Modification 1: The points {x 1 , . . . , x n } should be sampled as uniformly as possible on S. That is why we use the greedy sampling algorithm of Section 2.2 to select these points. To get an adaptive sampling, one could use a varying speed function, as shown in Figure 2.2 (see also [START_REF] Peyré | Geodesic Remeshing Using Front Propagation[END_REF]] for a curvature-based adaptation).

Modification 2: The K-neighborhood N i of each point should be computed using the geodesic distance and not the Euclidean one. This can be done very quickly using a local front propagation.

Modification 3: The matrix C(x) of equation ( 4) can not be computed anymore using dot products. Instead, following [START_REF] Roweis | Nonlinear Dimensionality Reduction by Locally Linear Embedding[END_REF] (pairwise LLE), we propose the following formula

-2C(x) i,j def. = d(x m(i) , x m(j) ) 2 - 1 K K k=1 d(x m(i) , x m(k) ) 2 - 1 K K k=1 d(x m(k) , x m(j) ) 2 + 1 K 2 K k,l=1 d(x m(k) , x m(l) ) 2 ,
that only uses geodesic distance information. This formula is equivalent to (4) in the Euclidean setting.

Extending the Map

The three modifications proposed in the previous section allow us to find the location of

x i = f (x i ) = (f 1 (x i ), f 2 (x i ))
T ∈ R 2 for each base point x i . To compute the whole map f , we need to interpolate the location of f (x) = (f 1 (x), f 2 (x))

T for each point x ∈ S, using the known locations f (x i ). This problem has been addressed very recently in [START_REF] Bengio | Out-of-Sample Extensions for LLE, Isomap, MDS, Eigenmaps, and Spectral Clustering[END_REF], by recasting it into a unified framework of eigenvector learning, common to many dimensionality reduction methods.

To extend f , we use the fact that vectors {f 1 (x i )} n i=1 and {f 2 (x i )} n i=1 are eigenvectors of the symmetric matrix M = (W -Id)

T (W -Id) (with eigenvalues λ 1 and λ 2 ). In the continuous setting, this matrix becomes a symmetric kernel M (x, y) for each point x, y in S. Matrix multiplication by M is then replaced by

ϕ → M ϕ(x) def. = S ϕ(y) M (x, y)dy, ( 5
)
where ϕ is any mapping from S to R. Using this remark, it is natural to suppose that the continuous maps f 1 and f 2 are eigenfunctions of the operator defined by equation ( 5) for the same eigenvalues λ 1 and λ 2 . This implies that we can compute them using a Nyström-like formula

f 1 (x) = 1 λ 1 S f 1 (y) M (x, y)dy ≈ 1 nλ 1 n i=1 f 1 (x i ) M (x i , y), (6)
and similarly for f 2 . Since the f 1 (x i ) are known, we just need to setup our kernel M . The only constraint is that for all y ∈ S, M (x i , y) should be easy to compute, e.g., it should only involve already computed distances such as d(x j , y) for x j ∈ N i . This can be done in a straightforward manner by first setting the weights for y:

w(x i , y) def. = 1 Ni (y) q (C(x i ) -1 ) k,q p,q (C(x i ) -1 ) p,q with y = x m(k) ,
and then defining the kernel:

M (x, y) def. = w(x, y) + w(y, x) - k w(x k , x)w(x k , y),
We can check that for base points, we retrieve the original matrix up to a substraction of the identity, i.e., M (x i , x j ) = δ i,j -M i,j . This shift is only here to avoid a singularity along the diagonal and does not modify the computation. This shows that we can extend the map f to a new point x using only some local distance information between x and its neighborhood in {x 1 , . . . , x n }. Furthermore, most of the time, this information is already available from previous front propagations performed to flatten {x 1 , . . . , x n }.

Figure 7 shows the flattening of one half of a human head. Even with a large patch that contains holes, our method gives very good results (no face flip) with only 100 base vertices.

Figure 8 shows the influence of the number of base points on the flattening. Even with only 20 points, the resulting embedding is nearly smooth except at the border of the mesh, and with 100 points, we get a perfectly smooth flattening.

Results and Discussion

Texturing of a Complex Model. To perform texture mapping on a complex 3D mesh, a segmentation step is required to first cut the model into disk-shaped charts. Although the study of this step is outside the scope of this paper, we note that the notion of Voronoi cells is often to perform mesh partition, as introduced in [START_REF] Eck | Multiresolution Analysis of Arbitrary Meshes[END_REF]. We choose to use a scheme based on a weighted geodesic distance [START_REF] Peyré | Surface Segmentation Using Geodesic Centroidal Tesselation[END_REF], since its continuous nature is clearly related to our flattening approach.

On Figure 9 one can see the whole pipeline in action. This includes first a centroidal tessellation of the mesh, then the extraction and flattening of each cell, and lastly the texturing of the model.

Computation Times. For our tests of the flattening procedure, we have chosen to use a fixed number of points (200 points), since the geometric complexity of the meshes was almost constant. The parameterization of [START_REF] Desbrun | Intrinsic Parameterizations of Surface Meshes[END_REF] is implemented using the boundary-free formulation (Neumann condition).

Table 1 shows the complexity of the algorithms mentioned in the paper, for a mesh of 10k vertices. The constant A is the number of steps in the gradient descent for the localization of the intrinsic center of mass, which is about A = 8 for 10k vertices. The constant B represent the number of base points, which is n/100 in our tests. This clearly shows the speed up that Geodesic LLE can bring over global methods such as [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF]. This is confirmed by the running times reported in table 2. For large meshes, the stability of our method is an advantage over the approaches based on large linear system such as [START_REF] Desbrun | Intrinsic Parameterizations of Surface Meshes[END_REF], for which it is difficult to ensure the convergence of the conjugate gradient. Discussion. The complete texturing of the David mesh (700,000 vertices) shown on Figure 9 clearly enlightens the strengths of our approach:

• The resulting flattening map is smooth, with no face flip (at least on this model). This can be seen on the close-up of the flattened domain. • The whole texturing procedure takes 740s, which shows an important speed up with respect to previous methods. • Our scheme is more local than the flattening procedure of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF], but it does not reach the per-vertex resolution of classical methods such as [START_REF] Desbrun | Intrinsic Parameterizations of Surface Meshes[END_REF]. This enables both fast computations and respect of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF]] E c = 0.8

M 1 M 2 M 3 [Desbrun et al., 2002] E c = 0.9 E c = 1.5 E c = 2.5 (conformal) E a = 1.2 E a = 3 E a = 10.4 [Desbrun et al., 2002] E c = 1.4 E c = 3.0 E c = 8.3 (authalic) E a = 0.6 E a = 1.1 E a = 3.5 [
The flattening (MDS) E a = 0.9 is not valid Our scheme

E c = 1.1 E c = 1.7 E c = 6.5 (GeodesicLLE)
E a = 0.9 E a = 1.6 E a = 5.5

Table 3. Area and angular distortion for various schemes the small scale variations (bumps or noise), which is not the case of [START_REF] Zigelman | Texture Mapping Using Surface Flattening via Multi-dimensional Scaling[END_REF].

Notice that there is no theoretical guarantee on the validity of the flattening. The only cases where face flips occurs is on patches with huge isoperimetric distortion. We believe that this is not a real issue since such degenerate cases can be easily detected and fixed (for example by subdividing the region). It is important to note that classical methods also face similar problems. In [START_REF] Levy | Least Squares Conformal Maps for Automatic Texture Atlas Generation[END_REF] a cut is performed to accelerate the convergence of the system resolution and ease the parameterization in regions with a sock-like shape. In practice however, our segmentation algorithm ensures that patches that need to be flattened do not contain high curvature variations and the whole process performs very well with no face flip. Distorsion Measures. To support our claim that our flattening scheme performs a trade-off between conservation of area and conservation of angle, we have performed some test (see table 3). We used 3 finger-like meshes M 1 , M 2 and M 3 with increasing isoperimetric distortion. On each face x we compute the eigenvalues (s 1 , s 2 ) of of the Jacobian of the parameterization map (linearly evaluated). Locally, conformality is characterized by s 1 = s 2 and conservation of area by s 1 s 2 = 1. As a conformal metric, we use

E c (M ) 2 = 1 A x s 1 (x) s 2 (x) + s 2 (x) s 1 (x) -2 2 A(x)
and as an equi-areal metric we use

E a (M ) 2 = 1 A x s 1 (x)s 2 (x) + 1 s 1 (x)s 2 (x) -2 2 A(x)
where A(x) is the area of a face x, A is the total area.

Conclusion

We have described new algorithms to perform the remeshing and the flattening of a genus-0 triangulated manifold. The main tool that allows having a fast algorithm is the fast marching on a triangulated mesh, together with some improvements we added. We have presented a fast algorithm for remeshing of a surface with a uniform or adaptive distribution. This is based on iteratively choosing the farthest point according to a weighted distance on the surface. We introduced a geodesic version of Locally Linear Embedding that is able to perform fast computations on a given set of points, and to extend the embedding to the rest of the mesh in a transparent manner. The resulting flattening is smooth and achieves a desirable trade-off between conservation of angle and area.

Figure 2 .

 2 Figure 2. Front Propagation (on the left), level sets of the distance function and geodesic path (on the right).

Figure 3 .

 3 Figure 3. An overview of the greedy sampling algorithm.

Figure 4 .

 4 Figure 4. Geodesic remeshing with an increasing number of points.

Figure 5 .

 5 Figure 5. Iterative insertion of points in a square.

Figure 6 .

 6 Figure 6. Uniform versus curvature-based sampling and remeshing.

Figure 7 .Figure 8 .

 78 Figure 7. The original model, texture on the flattened domain, and on the 3D mesh. 5 base points 10 base points 20 base points 100 base points

Figure 9 .

 9 Figure 9. Texturing of the David.

F

  . Marching Greedy sampl. 1 Lloyd iter. Zigelman02 Geodesic LLE Complexity n log(n) n log(n) 2 An log(n) Bn log(n) + B 3 n log(n) + B 2

Table 1 .

 1 Complexity of the algorithms

	Nbr.vertices [Zigelman et al., 2002] [Desbrun et al., 2002] Geodesic LLE
	1,000	7s	3s	5s
	10,000	55s	25s	28s
	100,000	440s	210s	150s
	700,000	2160s	1320s	740s

Table 2 .

 2 Comparison of flattening algorithms


Thomas Tillocher

Philippe Lefaucheux

Xavier Mellhaoui

Remi Dussart

Pierre Ranson
email: pierre.ranson@univ-orleans.fr

Mohamed Boufnichel

Lawrence J Overzet

Oxidation threshold in silicon etching at cryogenic temperatures

Keywords:

teaching and research institutions in France or abroad, or from public or private research centers.

1

The cooling of patterned silicon wafers is necessary to etch high aspect ratio structures in SF 6 /O 2 plasmas: this is the so called cryogenic process. At temperatures as low as -100℃, a SiO x F y passivation layer is formed on the sidewalls and acts as an etch inhibitor [1]. Since ion bombardment prevents its formation at the bottom of the patterns, anisotropic etching can then be achieved.

The formation of a similar layer has been previously reported by several authors during the etching of maskless silicon wafers [2][3][4][5]. They worked at room temperature in CF 4 /O 2 and SF 6 /O 2 gas mixtures. Their XPS measurements confirmed the presence of silicon, oxygen and fluorine in this layer. It also enabled the computation of its stoichiometry and its thickness as a function of the O 2 content in the plasma. As the O 2 proportion was increased in percentage, the surface was increasingly oxidized and the stoichiometry gradually changed from SiF 2 without oxygen, to SiOF 2 at middle proportions and finally to SiOF [2,3]. However, the film thickness did not have a monotonic evolution: it first decreased with increasing O 2 percentage and then began to increase again. The thickness of the SiO x F y layer, or in other words the level of oxidation, directly impacted the etch rate as well. Thus, above an O 2 proportion threshold, the layer became thick enough to strongly reduce the etch rate. This agreed well with the studies of Puech et al. who showed the existence of such an etch threshold at cryogenic temperatures for low oxygen concentrations [6,7]. They also noticed that increasing either the temperature or the self-bias makes the threshold occur at higher oxygen contents.

In this article, we present our studies of the threshold of oxidation using bare silicon wafers and mass spectrometry as a function of the O 2 proportion at room temperature and especially at low temperatures. We also highlight the role of both source power and chuck self-bias ; the increase of which can shift the threshold to higher oxygen proportions. These results are then discussed with the help of a short model modified from that in reference [8].

An industrial ICP (Inductively Coupled Plasma) reactor (Alcatel 601 E) only dedicated to cryogenic etching has been utilized for the experiments detailed in this letter. SF 6 /O 2 plasmas are created in the alumina tube of the source. It can be powered up to 3 kW.

The maximum gas flows are 1000 sccm for SF 6 and 70 sccm for O 2 . A magnetic field of 15 Gauss confines the plasma in both the source and the diffusion chamber. The reactor is equipped with an electrostatic chuck cooled with liquid nitrogen and warmed up with a heating element. It enables us to regulate the temperature between -150℃ and +40℃ with This chuck is independently rf-biased with a 500 W power supply, which is voltage-regulated.

A Quadrupole Mass Spectrometer (Balzers QMS 200) is mounted 5 cm above the bottom of the diffusion chamber [9]. Molecules are extracted from the plasma through a homemade limiter that reduces the detection of molecules recombined after their entrance in the spectrometer. Two signals of interest for our experiments are SiF + 3 and F + . Indeed, the SiF + 3 line dominantly comes from the fragmentation of SiF 4 which is the main etch product; because Si 2 F 6 is produced in negligible quantities [3]. Thus the SiF + 3 signal is reasonably proportional to the etch rate. The interpretation of the F + signal is not so obvious since it can result from the fragmentation of many species in the ionizer. Consequently, it does not only reflect the concentration of F radicals.

Experiments were first carried out at 3 Pa, 700 W source power and no bias. In this case, the O 2 /SF 6 (feed gas flow ratio) percentage was varied from 0 to 35%. The pressure in the spectrometer was approximately 7.10 -6 mbar. 6" bare or SiO 2 coated silicon wafers have been used. They are 400 µm thick and (100) oriented. The SiO 2 coating, if present, is nearly 0.8 µm thick. Fig. 1a shows the SiF + 3 peak intensity as a function of the oxygen content for processing bare silicon wafers at room temperature and -85℃ and a SiO 2 coated wafer at room temperature. In this last case, the signal is quite small because SiO 2 is etched much slower than silicon. It also weakly increases with oxygen proportion. For the etching on a bare silicon wafer at ambient temperature, the SiF + 3 signal first increases up to maximum at 5%. Then, it slowly decreases down to a very sharp drop at 25%. This is the threshold where oxidation begins to dominate etching. This is clearly shown on fig. 1a: the SiF + 3 signal becomes as low as that for SiO 2 hence the etching is strongly reduced. This is in agreement with the results in ref [2,3,5,10].

At -85℃, the trend is slightly different: the threshold is reached for an oxygen percentage around 2.5%. Then, the SiF + 3 signal increases more quickly than for the SiO 2 masked wafer. This suggests that the layer is not as strong as an oxide film: as this layer already contains fluorine, it can be etched more easily and rapidly. Moreover, the SiF + 3 signal maximum at -85℃ is not as high as it is at ambient temperature. This probably does not mean that the fluorine reactivity with silicon is reduced at cryogenic temperatures, but that less SiF 4 is desorbed. This is in agreement with a passivation mechanism detailed in ref [9] : O radicals can then react with SiF x sites and create a passivation layer. This is more efficient when the substrate temperature is lower, as can be inferred from fig. 2a. The oxidation threshold quasi linearly increases from -85℃ to 0℃ which suggests that atomic oxygen has a higher sticking coefficient on SiF x sites at lower temperature. This is consistent with studies which show a higher initial sticking coefficient of molecular oxygen on silicon at cryogenic temperatures [START_REF] Dabroski | Silicon Surfaces and Formation of Interfaces[END_REF][START_REF] Engel | [END_REF], but we have not been able to find similar references for atomic oxygen adsorption on fluorinated silicon surfaces.

These trends are consistent with the F + peak intensity under the same experimental conditions, displayed on fig. 1b. When the oxidation threshold is reached, etching is strongly reduced and the consumption of fluorine is reduced. As a result the F + signals sharply increase just after the thresholds determined from the SiF + 3 lines. The increase of the etch product signal (SiF + 3 ) at low oxygen concentrations can be explained in terms of SF 6 /O 2 chemistry. The addition of oxygen in SF 6 plasmas promotes fluorine creation [13,14]. Such an effect is clearly shown on fig. 1b: the F + peak intensity increases with oxygen when the substrate has a SiO 2 mask. Indeed, reaction of oxygen with SF x species reduces their recombination with F atoms (eq. 1) and also increases F production (eq. 2).

SF x + O wall --→ SOF x (x=1-4) (1)
SF x + O -→ SOF x-1 + F (x=1-5) (2)
The first increase of the etch rate can also be attributed to sulfur atoms which are reported to have a retarding effect on oxidation at room temperature [2,10]. For low oxygen contents, sulfur is present on the surface and prevents oxygen atoms from adsorbing. But for higher concentrations, sulfur is no longer efficient [2]. In our case, a very low temperature seems to reduce this effect since the oxygen adsorption probability is increased. This also explains why the SiO x F y film is created for lower percentages of oxygen.

Physical mechanisms have non negligible effects on the oxidation as well. Figure 2b shows the oxidation threshold as a function of the source power at a substrate temperature of -85℃. The other plasma parameters are not changed. The threshold seems to increase linearly with the power from 700 W to 2000 W.

The observation is the same on fig. 2c, which displays the oxidation threshold against the bias, at the same temperature: we also have a quasi linear evolution. After each of these experiments, columnar microstructures (CMS) had appeared on the substrates. This is an evidence of the formation of a strong passivation layer at high oxygen contents [15]. And these CMS appear for higher oxygen content when the self-bias is increased because oxidation is delayed.

Such results can be interpreted in terms of density of power (P i ) deposited on the surface by the positive ions expressed as [START_REF] Lieberman | Principles of plasma discharges and materials processing[END_REF]:

P i = Γ i E i = Γ i (ǫ i + q(V p -V b ) + 2kT e ) (3)
Γ i is the ion flux and is proportional to the ion density which can be increased with the source power. The second term, E i , is composed of the ionization energy ǫ i , the energy acquired in the sheath q(V p -V b ) and the electron energy. The q(V p -V b ) quantity can be controlled with the bias V b . Consequently, increasing either the source power or the selfbias contributes to increase either the ion flux or the ion energy (E i ) hence the density of power transferred to the surface by ions and especially to the SiO x F y layer. This is consistent with the disappearance of the passivation layer when the substrate is heated: the ion bombardment effectively warms up the surface by supplying enough energy to desorb the film [9]. Consequently, this removal mechanism is different from sputtering.

We now consider a modified version of the model proposed by Mogab et al. and also used by Campo et al [2,3,8]. It should take into account the mechanisms involved at cryogenic temperatures:

dθ O dt = k a (T )(1 -θ O )N O -k d θ O N F -k i (E i )θ O N i (4)
θ o is the oxygen coverage of the silicon , k a is the adsorption of oxygen on silicon rate constant and N o the oxygen concentration. k d is the rate constant associated to the desorption of oxygen by the etching and N f the fluorine density. This equation illustrates the fact that oxygen and fluorine are in competition for the adsorption sites on silicon through the first two terms when no bias is applied to the substrate. In such a case and at steady state, the oxygen coverage is given by:

θ O = 1 1 + k d N F kaN O (5)
The limiting case, θ o = 1, is obtained when oxidation dominates etching (k d N F ≪ k a N O ) hence for high oxygen concentrations or for higher values of k a . This is certainly achieved at very low temperatures since the oxidation is then more efficient (fig 2.a). That's why we consider that k a depends on the temperature. The last term in eq. 5 represents the ion contribution in the removal of the passivation layer. The rate constant, k i , which depends on the ion energy E i , is used to describe the desorption of oxygen by ion bombardment and N i is the ion density. θ o now becomes :

θ O = 1 1 + k d N F +k i N i kaN O (6)
An increase in the ion contribution can cause the oxygen coverage to decrease. As it has been previously shown, this can be achieved by increasing either the source power or the bias. This is in good agreement with the experimental results. As oxygen plays a key role in the formation of the passivation layer, we can assume that eq. 4 can also be used to describe its behavior.

These experiments have given new elements about the oxidation of silicon. SiO x F y films grow when the O 2 percentage is sufficiently large. An oxidation threshold occurs for lower oxygen proportions when the substrate is cooled. This threshold is reached for higher O 2 contents when either the source power or the self-bias is increased. We have presented a model which is in good agreement with both our results and a passivation mechanism for cryoetching proposed in ref [9].

FIG. 1 : 8 FIG. 2 :

 182 FIG. 1: (a) SiF + 3 signal peak intensity versus the oxygen percentage in three cases for 3 Pa, 700 W and without bias (b) F + signal peak intensity versus the oxygen percentage in the same conditions

The authors acknowledge STMicroelectronics (Tours, France) for helpfull discussions. This work is supported by STMicroelectronics and by the French Governement.