Jean-Daniel Boissonnat

Kan-Le Shi

Jane Tournois

Mariette Yvinec

Anisotropic Delaunay Meshes of Surfaces

Keywords: Mesh generation, anisotropic meshes, Delaunay triangulation, Delaunay refinement

Anisotropic simplicial meshes are triangulations with elements elongated along prescribed directions. Anisotropic meshes have been shown to be well suited for interpolation of functions or solving PDEs. They can also significantly enhance the accuracy of a surface representation. Given a surface S endowed with a metric tensor field, we propose a new approach to generate an anisotropic mesh that approximates S with elements shaped according to the metric field. The algorithm relies on the well-established concepts of restricted Delaunay triangulation and Delaunay refinement and comes with theoretical guarantees. The star of each vertex in the output mesh is Delaunay for the metric attached to this vertex. Each facet has a good aspect ratio with respect to the metric specified at any of its vertices. The algorithm is easy to implement. It can mesh various types of surfaces like implicit surfaces, polyhedra or isosurfaces in 3D images. It can handle complicated geometries and topologies, and very anisotropic metric fields.

Maillages de surfaces Delaunay anisotropes

Résumé : Les maillages anisotropes simpliciaux sont des triangulations dont les éléments sont étirés suivant certaines directions imposées. Les maillages anisotropes sont connus pour être bien adaptés à l'interpolation de fonctions ou à la résolution d'équations aux dérivées partiellles. Ces maillages peuvent aussi améliorer notablement la précision de l'approximation d'une surface. Etant donnée une surface S, munie d'un champs de tenseurs qui définit la métrique en tout point de la surface, nous proposons un nouvel algorithme pour générer un maillage anisotrope qui approxime S par des triangles dont les formes s'adaptent à la métrique locale. L'algorithme repose sur les concepts bien établis de triangulation de Delaunay restreinte et de raffinement de Delaunay et offre des garanties théoriques. L'étoile de chaque sommet dans le maillage est formée par des triangles de Delaunay pour la métrique du sommet central. Chaque triangle a un bon rapport d'aspect dans la métrique attachée à chacun de ces sommets. L'algorithme est facile à programmer. Il permet de mailler diffrents types de surfaces, comme des surfaces implicites, des polyèdres ou encores des isosurfaces dans des images 3D. L'algorithme peut traiter des surfaces de géométrie ou topologie complexe, il peut aussi prendre en compte des anisotropies très prononcées.

Mots-clés : Génération de maillages, maillages anisotropes, Triangulation de Delaunay, raffinement de Delaunay

Introduction

Anisotropic simplicial meshes are triangulations with elements elongated along prescribed directions. Anisotropic meshes have been shown to be particularly well suited for interpolation of functions and solving PDE's [START_REF] Simpson | On optimal interpolation triangle incidences[END_REF][START_REF] Richard | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF][START_REF] Mirebeau | Optimal meshes for finite elements of arbitrary order[END_REF]. They can also significantly enhance the accuracy of a surface representation if the anisotropy of the mesh conforms to the curvature tensor of the surface [START_REF] Heckbert | Optimal triangulation and quadric-based surface simplification[END_REF]. In this paper, we consider the problem of sampling and meshing a surface S endowed with a metric tensor field that varies over S. The metric tensor can be chosen arbitrarily provided that it satisfies a continuity condition over the surface S. The metric tensor may, for instance, be the curvature tensor of S or the Hessian matrix of a function defined over S.

Background

Due to its practical importance, several methods have been proposed in the literature to address the problem of approximating surfaces using anisotropic surface meshes. The case of a parametric surface S can be handled by working in the 2-dimensional space U of the parameters of S, using an appropriate metric field to account for the deformation induced by the mapping U → S. The problem is then 2-dimensional [BGH + 97, BH96, LTÜ99]. Several methods have also been reported to remesh a polyhedral surface according to some specified metric tensor. Alliez et al. [ACSD + 03] first estimate the principal direction fields and trace a network of lines of curvature from which they deduce a quad-dominant anisotropic mesh. Jiao et al. [START_REF] Jiao | Anisotropic mesh adaptation for evolving triangulated surfaces[END_REF] use the quadratic metric tensor of Heckbert and Garland [START_REF] Heckbert | Optimal triangulation and quadric-based surface simplification[END_REF] and propose a set of operations to allow anisotropic adaptation of a triangular mesh to static or evolving surfaces. The case of implicit surfaces can be handled by first polygonizing the surface in a standard way and then remeshing the obtained surface using one of the anisotropic methods mentioned above. Although this approach is clearly not optimal, very few methods have been proposed to directly generate anisotropic meshes on implicit surfaces. One exception is the work of Azernikov and Fischer [START_REF] Azernikov | Anisotropic meshing of implicit surfaces[END_REF] who propose a grid-based approach to mesh implicit surfaces according to their curvature tensor. The method uses a deformed grid and produces quad-meshes. The anisotropic grid is obtained by deforming a cartesian grid while maintaining the structured topology, which is a limitation when considering complex shapes.

In order to settle more rigorous foundations and to precisely characterize anisotropic triangulations, several authors have introduced anisotropic Voronoi diagrams (AVD). An anisotropic triangulation, called anisotropic Delaunay triangulation (ADT), is then defined as the dual of the AVD. Various AVD and ADT have been proposed. Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and voronoi diagrams for riemannian manifolds[END_REF] introduced Delaunay triangulations and Voronoi diagrams for Riemannian manifolds. This approach requires to compute geodesic paths and intrinsic balls, which may be quite complicated in practice. To overcome this issue, two approaches have been proposed that approximate the geodesic distance between a given site and any point by considering the metric constant and equal to the metric at the site [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] or at the point [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF][START_REF] Lévy | Lp centroidal voronoi tessellation and its applications[END_REF]. Although nice results have been reported, it is recognized in [START_REF] Du | Anisotropic centroidal voronoi tessellations and their applications[END_REF] that "more precise theoretical analysis are needed in identifying conditions on the metric, geometry, and the generator distributions that guarantee the well-defined duality between the AVDs and the ADTs". The Labelle and Shewchuk approach [LS03] is easier to analyze and, in the 2-dimensional case, the authors have proposed a refinement algorithm that can provably produce anisotropic meshes. This approach has been extended by Cheng et al. [START_REF] Cheng | Anisotropic surface meshing[END_REF] to produce anisotropic meshes of surfaces embedded in R 3 . The algorithm is however very complicated and no implementation has been reported. Although some progress has been recently reported in [START_REF] Canas | Orphan-free anisotropic Voronoi diagrams[END_REF], it remains unclear under which conditions an anisotropic Voronoi diagram admits a dual embedded triangulation for 3 or higher dimensional domains.

A different approach has been proposed by Boissonnat, Wormser and Yvinec [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. They introduced a new variant of anisotropic Delaunay meshes, called locally uniform anisotropic meshes. Such meshes are not defined by duality from some anisotropic Voronoi diagram. Instead, it is defined as a simplicial mesh in which the star of each vertex is Delaunay for the metric attached to the vertex. It is shown that, when the metric field satisfies some continuity property, such a mesh can be obtained for bounded polyhedral domains of R 2 or R 3 with no sharp edges. The extension of the approach to surfaces has not been considered so far and no implementation has been reported for 3-dimensional domains.

Our Contributions

Our work uses the notion of locally uniform anisotropic mesh defined in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], combined with the notion of Delaunay triangulation restricted to a surface. As in the previous approaches, we assume that the anisotropy is prescribed by a metric field that associates to each point p of the surface, a symmetric positive definite 3 × 3 matrix M p describing the metric at point p. Given a set of (punctual) sites V and v ∈ V , it is easy to compute the 3D Delaunay triangulation Del v (V ) for the metric M v of v. Indeed, Del v (V ) is just the image of an Euclidean Delaunay triangulation under a stretching transformation. For each site v ∈ V , we consider the star T v of v in Del v (V ), and the surface star S v . The star T v is the subcomplex of Del v (V ) formed by the tetrahedra that are incident to v, while the surface star S v is the subcomplex formed by the triangular facets incident to v and whose dual Voronoi edges intersect S. The surface star S v is thus a collection of triangles incident to p.

In general, there are inconsistencies among the surface stars of the sites : a facet τ , appearing in the surface stars of some of its vertices, may not appear in the surface stars of all of them. As a result, the set of surface stars of V does not form a triangulated surface.

In [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF][START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF], we showed that it is possible to refine the set V of sites until there is no more inconsistencies among the set T (V ) of 3D stars. This leads to a locally uniform anisotropic 3D mesh in which the star of each vertex is Delaunay with respect to the metric of this vertex. In the present paper, we are interested in generating surface meshes. The algorithm will incrementally compute a set V of sites on the surface and will maintain the associated surface star set S(V ). The set V is refined until there is no inconsistencies among the facets of S(V ). In contrast to [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], we do not consider inconsistencies among the 3D stars. When there is no more inconsistencies among the surface stars, the facets of these surface stars form a locally uniform anisotropic surface mesh that approximates S and conforms to the specified metric field. This algorithm is provably correct and much more efficient than the 3D anisotropic mesh generator of [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF].

In addition to conforming to the given metric field, this new method has several notable advantages over previous methods.

-The algorithm relies on the well-established concepts of restricted Delaunay triangulations and Delaunay refinement and comes with theoretical guarantees.

-It can handle complicated geometries and topologies, as well as very anisotropic metric fields.

-It can handle various types of surfaces, e.g. implicit surfaces, polyhedra, isosurfaces in 3D images.

-It is easy to implement and robust since it relies on the usual Delaunay predicates (applied to some stretched spaces).

Inria

-The facets of the output mesh have a good aspect ratio as measured with respect to the metric field.

The idea of maintaining independent stars for each vertex of a mesh has been first proposed by Shewchuk [START_REF] Shewchuk | Star splaying: an algorithm for repairing Delaunay triangulations and convex hulls[END_REF] for maintaining triangulations of moving points. This idea has also been used to build the dual of an anisotropic Voronoi diagram as defined by Labelle and Shewchuk [START_REF] Schoen | Robust, guaranteed-quality anisotropic mesh generation[END_REF]. The idea of stitching stars restricted to a surface can be found in the work on the tangential complex [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF]. The context is however different from ours since no anisotropy was considered and the main motivation was to reconstruct a k-manifold embedded in a d-dimensional space without constructing any global d-dimensional triangulation, which is not a critical issue in 3 dimensions (see however the discussion in Section 5.3). A major technical difference with our paper is that the stars considered in [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF] are all contained in a global Euclidean d-dimensional Delaunay triangulation, which is not the case here. Combining ideas from the present paper and [START_REF] Boissonnat | Manifold reconstruction using Tangential Delaunay Complexes[END_REF] would allow to define and generate anisotropic meshes of smooth submanifolds embedded in high dimensional spaces.

Preliminaries

Anisotropic Metric

An anisotropic metric in R 3 is defined by a symmetric positive definite quadratic form represented by a 3 × 3 matrix M . The distance between two points a and b, as measured using a metric M , is defined as

d M (a, b) = (a -b) t M (a -b).
This definition can be extended to M -lengths and M -areas.

In the following, we often use the same notation, M , for a metric and the associated matrix in a given basis. Given the symmetric positive definite matrix M , we denote by F M any matrix such that det(F M ) > 0 and

F t M F M = M . Note that d M (a, b) = (a -b) t F t M F M (a -b) = F M (a -b) (1)
where the notation . stands for the Euclidean norm.

Given some metric M , an M -sphere C M (c, r), with center c and radius r, is defined as the set of points p such that d M (c, p) = r. Likewise, an M -ball B M (c, r), is defined as the set of points p such that d M (c, p) ≤ r. Note that an M -sphere is an Euclidean ellipsoid, with its axes aligned along the eigenvectors of M .

Given a simplex τ in R 3 and a metric M , we define the M -circumsphere C M (τ ) as the smallest M -sphere that circumscribes τ . The M -circumball B M (τ ) is the M -ball bounded by C M (τ ) and the M -circumradius r M (τ ) of a simplex τ is the radius of C M (τ ).

Distortion between Metrics

The following definition is due to Labelle and Shewchuk [START_REF] Labelle | Anisotropic voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]. Given two metrics M and N , and their square-roots F M and F N , the distortion between M and N is defined as

γ(M, N ) = max{ F M F -1 N , F N F -1 M },
where . is the matrix norm operator associated with the Euclidean norm, i.e. for a 3 × 3 square matrix A, A = sup x∈R d

Ax

x . Observe that γ ≥ 1 and

γ = 1 iff M = N . A fundamental RR n°8400
property of the distortion γ(M, N ) is to relate distances d M and d N . Specifically, for any x, y, we have

1/γ(M, N ) d M (x, y) ≤ d N (x, y) ≤ γ(M, N ) d M (x, y).
In the rest of the paper, we consider a compact closed surface S of R 3 endowed with a continuous metric field {M p , p ∈ S} defined over S. For convenience, we will simply write X p instead of X Mp when a quantity X is measured using metric M p . Hence, we will write for instance F p for F Mp and d p (a, b) for d Mp (a, b).

Given two points p and q of S, the relative distortion between p and q is defined as γ(p, q) = γ(M p , M q ). Note that Γ = sup x,y∈S γ(x, y) is finite since S is compact.

Let γ 0 be some positive constant called the distortion bound. The distortion radius at p, dr(p, γ 0 ) is the upper bound on distances ℓ such that for all q and r in S,

max(d p (p, q), d p (p, r)) ≤ ℓ ⇒ γ(q, r) ≤ γ 0 .
It can be shown that the distortion radius satisfies the following continuity property [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF] :

1 γ(p, q) (dr(p, γ 0 ) -d p (p, q)) ≤ dr(q, γ 0 ) dr(q, γ 0 ) ≤ γ(p, q) (dr(p, γ 0 ) + d p (p, q)) (2)

Sizing Field

To control the distorsion of the metric on mesh elements, we use a sizing field that is pointwise smaller than the distorsion radius. In addition to controlling the distortion of the metric, it is important to control the topology of the output mesh. To do so, we adapt the notion of local feature size to the anisotropic setting. The (anisotropic) local feature size at a point p ∈ S, denoted by lfs(p), is defined as the M p -distance from p to the M p -medial axis of S. It can be proved that the anisotropic local feature size enjoys a continuity property analoguous to (2).

To control both distortion and topology, we use a sizing field sf(p, γ 0 ) defined as sf(p, γ 0 ) = min(lfs(p), dr(p, γ 0 )). Plainly, the sizing field sf(p, γ 0 ) satisfies a continuity property analog to (2). In fact, the user may include his own sizing requirement by using as sizing field any function provided that it is pointwise smaller than sf(p, γ 0 ) and satisfies a continuity property analog to (2).

Stars and Refinement

Anisotropic Star Set

Given a metric M and a set of punctual sites V , the Delaunay triangulation Del M (V ) for metric M is the triangulation of V such that the M -circumball of each tetrahedron is empty. By empty, we mean that the interior of the M -circumball contains no points of V . By (1), the Delaunay triangulation Del M (V ) of a finite set of points V for metric M is simply obtained by computing the Euclidean Delaunay triangulation of the stretched image F M (V ), and then stretching back the result using F -1 M . For each site v in V , we consider the Delaunay triangulation Del v (V ) of V for metric M v . We define the star T v of site v as the subcomplex of Del v (V ) formed by the tetrahedra that are incident to v. The collection of all the stars T v , v ∈ V , is called the star set of V . We denote it by T (V ).

Inria

We now adapt to the anisotropic setting the definition of a restricted Delaunay triangulation. The surface star of v, denoted by S v , is the subcomplex of T (v) formed by the facets of T (v) that are incident to v and whose dual Voronoi edge intersects S. The collection of the surface stars S v , for v ∈ V , is called the surface star set of V . We denote it by S(V ).

Let N be a metric and τ be a triangle with its three vertices p 0 , p 1 and p 2 on S. We define the N -axis of τ as the set of centers of all the N -circumballs of τ . If p 0 , p 1 and p 2 are sufficiently close to each other, and if τ is well-shaped with respect to metric N , the N -axis of τ intersects S. We call surface N -circumball of τ , and note B N (τ ), the smallest N -circumball of τ centered on S. A facet τ ∈ T v incident to v belongs to S v iff its surface M v -circumball is empty.

The facet τ with vertices p 0 , p 1 and p 2 on S has three surface circumballs B i (τ ) associated respectively to the metrics M i of each of its vertices. We define the distortion γ(τ ) of τ as the maximal distortion between any pair of points of S which are both inside the same surface circumball B i (τ ).

Inconsistencies

Two surface stars S v and S w are said to be inconsistent if some facet incident to v and w appears in only one of the two stars S v and S w . Any facet that appears in the stars of some of its vertices, but not in the stars of all of them, is also said to be inconsistent (see Figure 1). Our algorithm incrementally inserts new sites in V to resolve inconsistencies among the surface stars.

Quasi-Cosphericity

Let γ 0 > 1 be a bound on the distortion and M be a metric. Following [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF], we now introduce the notion of (γ 0 , M )-cosphericity and show its link to inconsistent simplices.

Let τ be a facet of some surface star and p ∈ V \ τ . We call U = (τ, p) a (γ 0 , M )-cospherical configuration iff there exists two metrics N and N ′ such that

1. γ(M, N ) ≤ γ 0 , γ(M, N ′ ) ≤ γ 0 and γ(N, N ′ ) ≤ γ 0 ; RR n°8400
2. p is included in the interior of the surface N ′ -circumball B N ′ (τ ) and not in the interior of the surface N -circumball B N (τ ).

If M is clear from the context, we simply say that U is a γ 0 -cospherical configuration and if both M and γ 0 are understood, we say that U is a quasi-cospherical configuration. See Figure 2 for an illustration where S is the plane.

a c b d N' (τ) N (τ) τ Figure 2: The configuration (τ, d) is quasi-cospherical because point d is outside the N -circumball B N (τ ) but inside the N ′ -circumball B N ′ (τ ).
The following lemma ensures that, when all simplices in the stars have small distorsion, inconsistencies only result from the occurence of quasi-cospherical configurations.

Lemma 1 (Inconsistencies and quasi-cosphericities) Let τ be an inconsistent facet that appears in star S v but not in star S w where v and w are two vertices of τ . If γ(τ ) < γ 0 , then there exists a site q ∈ S w such that the configuration (τ, q) is (γ 0 , M v )-cospherical.

Proof. Take N = M v and N ′ = M w . Because the distortion of τ is less than γ 0 , we have γ(v, w) ≤ γ 0 , and therefore the distorsions γ(M v , N ), γ(M v , N ′ ) and γ(N, N ′ ) are all less than γ 0 . Assume first that all the vertices of τ appear in S w . Then if τ ∈ S w , the surface M w -circumball B w (τ ) of τ contains some sites of S w in its interior. Consider now the case where some vertex v ′ of τ does not appear in S w . This means that v ′ w is not an edge of Del w (V ). Hence, any surface M w -ball circumscribing edge v ′ w, and in particular B w (τ ), contains some sites of S w . Thus, in both cases, there is some site q in S w such that q is included in the interior of B w (τ ) and not in the interior of B v (τ ). It follows from Definition 3.3 that (τ, q) is a (γ 0 , M v )-cospherical configuration. ✷ 4 Meshing Algorithm

Algorithm Outline

Let S be a given surface we want to mesh. The algorithm constructs the set of sites V by inserting new sites on S in a greedy way. While there remain inconsistent or bad facets in the surface star set, the algorithm selects one such facet τ and kills τ by inserting a new site on S called the refinement point of τ . Inserting a new point is performed by the following Insert procedure that maintains the star sets T (V ) and S(V ). We say that a tetrahedron σ of a star T w conflicts with a point p if p belongs to the M w -circumball of p.

Insert(p)

Inria 1. insert p in each star T w that contains a simplex in conflict with p and update the surface star S w ;

2. create the new stars T p and S p .

The refinement algorithm consists in applying the following rules with a priority order : Rule (i) is applied only if no Rule (j), j < i, can be applied. The algorithm ends when no rule applies anymore. The algorithm relies on procedure Insert to insert a new point in the data structures, and on procedure Pick_valid to select the location of the new site (see the next section).

If τ is a facet in some star S v , we write B v (τ ) for the surface M v -circumball of τ , and c v (τ ) and r v (τ ) for the center and the radius of this ball. The quality of facet τ with respect to metric M v is characterized by the M v -radius-edge ratio, ρ v (τ ), which measures the ratio between the M v -circumradius of τ and the M v -length of the edge of τ of shortest M -length.

The refinement rules depend on the constants α, ρ 0 and γ o .

1. Sizing field -Distortion :

If ∃v ∈ V and τ ∈ S v such that r v (τ ) ≥ α sf(c v (τ ), γ 0 ), then Insert(c v (τ ));
2. Shape -Radius-edge ratio :

If ∃v ∈ V and τ ∈ S v such that ρ v (τ ) > ρ 0 , then Insert (Pick_valid(τ, M v )); 3. Inconsistency : If ∃v ∈ V, τ ∈ S v such that τ is inconsistent, then Insert (Pick_valid(τ, M v ));
Assume that the algorithm terminates. Then, no rule applies anymore. In particular rule (3) does not apply which means that the star sets S(V ) includes no inconsistent simplex. All surface stars are consistent and can be stitched together to form a triangulated surface T with the property that the star of any vertex v in T is Delaunay for metric M v . Hence, T is a locally uniform anisotropic mesh of surface S. The next step is therefore to establish that the algorithm terminates (Section 4.3).

Picking Region and Hitting Sets

The simplest idea to kill a bad triangle τ in some surface star S v is to insert a new site at the center of the surface M v -circumball of τ . However, this simple strategy may lead to endless occurrences of quasi-cosphericities and inconsistencies. We present in this section an alternative strategy, inspired by the strategy introduced by Li and Teng [START_REF] Li | Generating well-shaped delaunay meshed in 3d[END_REF][START_REF] Li | Generating well-shaped d-dimensional delaunay meshes[END_REF] to avoid slivers in isotropic meshes and also used in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. The basic idea is to relax the choice of the refinement point of a bad facet. Instead of using the center of the surface circumball, the refinement point is picked from a small region around this center, and carefully chosen so as to avoid the formation of new quasi-cosphericities.

Let δ < 1 be a constant called the picking ratio. The M v -picking region of a facet τ in star S v , denoted byP v (τ ), is the intersection of the M v -ball B v (c v (τ ), δr v (τ )) with the surface S.

In fact, it is not possible, when choosing a refinement point in the picking region P v (τ ) of a triangle τ of S v to avoid the formation of any new quasi-cospherical configurations. The Pick_valid procedure will only avoid the creation of small quasi-cospherical configurations where the meaning of small depends on an additional parameter β ≥ 1 and is made precise below.

Assume we are choosing a new site p in the M v -picking region of some triangle τ in surface star S v . A subset σ of three sites of the current set V is said to hit p if the configuration U = (σ, p) is both :

1. (γ 0 , M )-cospherical for some metric M such that γ(M p , M ) ≤ γ 0 , 2. small, meaning that the smallest M -circumradius of all the triangles with vertices in U is less than βr v (τ ). A subset σ of three sites of V that hits some point p in the picking region P v (τ ) is called a hitting set of P v (τ ). We associate to each hitting set σ a so-called forbidden region that consists of the points of P v (τ ) hit by σ. A point of the picking region P v (τ ) is called a valid refinement point when it does not belong to any forbidden region. See Figure 3.

Note that the definition of a valid refinement point depends on the constants δ and β that respectively defines the size of the picking regions and bounds the minimum size of new quasicospherical configurations. It also depends on the constant γ 0 that defines quasi-cosphericity. The next lemma is fundamental for our algorithm.

Lemma 2 (Picking lemma) For any positive values of the parameters α ≤ 1, β ≥ 1, δ < 1 and ρ 0 , it is possible to choose a distortion parameter γ 0 close enough to 1, so that the valid Inria refinement points form a subregion of positive area inside any picking region considered by the algorithm.

The proof, only sketched here, is given in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF]. It uses a volume argument based on the two following facts. 1. The number of hitting sets of a picking region is bounded. 2. For given values of parameters α, β, δ and ρ 0 , the volume of the forbidden region associated to a hitting set tends to zero when γ 0 tends to 1.

To find a valid refinement point in the M v -picking region P v (τ ) of some bad k-simplex τ , the insertion algorithm calls the following Pick_valid procedure that randomly selects a point in the picking region P v (τ ) until a valid refinement point is found. This procedure depends on the constants γ o , δ and β, to be fixed later in Section 4.3.

Pick_valid(τ, M v ) :
1. Pick randomly a point p in the picking region P v (τ ) 2. If there exists a subset of three sites in V that hits p, then discard p and go back to 1.

Return p.

Lemma 2 ensures that the picking process will succeed with a positive probability at each trial. Observe also that Pick_valid takes care of small configurations only. Hence, checking if a point of some picking region is valid requires only to consider sites that are close to the picking region. Furthermore, the number of triplets we need to consider as possible hitting sets is bounded by a constant (see the proof of Lemma 2 in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF]).

Termination of the Algorithm and Quality of the Mesh

We now prove that, under some conditions on the parameters ρ 0 , γ 0 , δ and β, the number of inserted vertices is bounded, which implies that the algorithm terminates. Let us define the intersite distance between two sites p and q as d(p, q) = min(d p (p, q), d q (p, q)), and the insertion radius r(p) of a site p as the minimum intersite distance between p and the current set of sites, i.e. the set of sites that have been inserted before p. Recall that Γ = sup x,y∈S γ(x, y) is finite since S is compact. The proofs of the following lemmas, available in [START_REF] Boissonnat | Anisotropic Delaunay Mesh Generation[END_REF], are only sketched here.

Lemma 3 (Insertion radius lemma) If the two conditions

(1 -δ)ρ 0 Γγ 2 0 ≥ 2 and (1 -δ)β Γγ 2 0 (1 + δ) ≥ 2
hold, then there exists a constant C such that the insertion radius r(p) of any vertex p of the mesh is at least C sf(p).

The proof uses the fact that the new point is inserted close to the center of the surface circumball of a facet in a surface star, and therefore relatively far from the current set of sites. The precise proof is a case analysis that studies the three rules of the algorithm. From this lemma, we can bound from below the intersite distance between any two points p and q in the final mesh. We consider whether p or q has been inserted last. Lemma 3 and the continuity property (2) of sf then allow to prove the following lemma.

Lemma 4 (Smallest intersite distance) If the conditions of Lemma 3 are satisfied, the smallest intersite distance between a site p and all the other vertices of the final mesh is at least :

d(p, q) ≥ C (1 + C)Γ sf(p).
n ≤ γ 2 0 π 1 + 2Γ(1 + C) C 2 S d M p sf(p) 2 .
When the algorithm terminates, the surface stars are consistent and can be stitched so as to make a triangulated surface mesh Ŝ. Each triangular facet of Ŝ is Delaunay for the metric of any of its three vertices, conforms to the sizing field sf (Rule 1) and has a radius-edge ratio smaller than ρ 0 (Rule 2) when measured using the metric of any of its vertices. The radius-edge ratio being a fair quality measure for triangles [START_REF] Richard | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF], the last property ensures that the facets of the output triangulated surface are well-shaped with respect to the local anisotropy.

In addition, we can ensure that the topology of Ŝ is the same as the topology of S. This can be done by choosing a sizing field that is sufficiently small with respect to the local feature size of the surface. This condition, together with the bound on the radius-edge ratio of the facets, ensures that Ŝ is homeomorphic (and even ambient isotopic) to S. See for example [START_REF] Amenta | A simple algorithm for homeomorphic surface reconstruction[END_REF][START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF].

Implementation

The algorithm can handle different kinds of surfaces provided that we are given an oracle able to :

(1) compute a few initial points on the surface, (2) detect and compute the intersections between a given line segment and the surface.

A small set of initial points is enough. We create the initial stars and put the bad facets of the surface stars in a global queue containing all the facets to be refined. The refinement process can then start. Our implementation uses the cgal library [CGA].

The algorithm parameters

Our anisotropic surface mesh generation algorithm depends on five parameters : α, β, δ, ρ 0 and γ 0 . A priori, the parameters values should be chosen to satisfy the conditions of Lemma 3 and γ 0 should be set close enough to 1 to ensure the existence of valid refinement points in any picking region (Lemma 2). Note that, if the conditions of Lemma 3 are satisfied for a given choice of α, β, δ, ρ 0 and γ 0 , they are also satisfied for the same values of α, β, δ, ρ 0 and a lower γ 0 .

Conditions in Lemmas 2 and 3 are sufficient to ensure the termination of the meshing algorithm. However they are not necessary. These conditions arise from our proof of termination and our will to keep the proof simple makes them over-restrictive. In particular, the conditions in Lemma 2 depend on a global distortion bound Γ. This makes the proof easier but a more thorough proof would involve local distortion bounds. In practice, in our experiments we do not try to evaluate a global distortion bound, nor to respect the conditions of Lemmas 2 and 3. In fact, the meshing algorithm runs smoothly and yields good meshes for a quite large range of parameter values.

Inria

Parameter α should be chosen ≤ 1. This parameter is needed to prove the picking lemma. In practice, this parameter is not critical and has been set to 1 in all our experiments. Parameters β and δ, which define the picking region and the validity of a refinement point, have little influence on the quality of the final mesh. Those parameters have to be set in such a way that procedure Pick_valid. finds a valid refinement point in a reasonable amount of time. By definition, we must have δ < 1, and the second condition in Lemma 3 suggests to choose β > 2. In all experiments reported later in this paper, we use β = 2.5 and δ = 0.3. Further experiments showed that varying the value of β in the range [2, 5] and the value of δ in the range [0.1, 0.5] has basically no effect on the results (number of vertices, computing time, and quality of the output mesh). More information on these experiments are given in the appendix.

Parameter ρ 0 , on one hand, should be as small as possible to ensure that the simplices of the final mesh are well-shaped according to the metric field. In the isotropic case, the Delaunay refinement algorithm for smooth surface [START_REF] Chew | Guaranteed-quality mesh generation for curved surfaces[END_REF] is known to terminate for any value ρ 0 > 1. The first condition in Lemma 3 suggests to take ρ 0 > 2 . In all our experiments, we use ρ 0 = 3.

The size of the simplices in the final mesh mostly depends on the parameter γ 0 and on the sizing field . As defined above, the sizing field takes into account the variation of the metric (distortion radius) and the geometry of the surface (local feature size). Estimating the local feature size is a notoriously difficult problem, for which existing solutions can be used [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF][START_REF] Antani | Mesh sizing with additively weighted Voronoi diagrams[END_REF]. In this paper, the emphasis is put on the metric variation and, in our experiments, we assume the local feature size to be uniformly bounded from below by some parameter r 0 that we use to specify the desired size of the mesh simplices. To take into account the metric variation, we replace the distortion radius by a direct and simpler control on the distortion of the simplices. Specifically, Rule 1 is replaced in practice by the following rule :

1. Sizing field -Distortion :

If ∃v ∈ V and τ ∈ S v such that either r v (τ ) ≥ r 0 or γ(τ ) ≥ γ 0 , then Insert(c v (τ ));

Parameter γ 0 has to be greater than 1 and close enough to 1. The influence of the value of parameter γ 0 is complex. On one hand, choosing for γ 0 a value very close to one triggers many applications of the Sizing field -Distorsion rule and leaves only a few inconsistencies to be removed by the Inconsistency rule. Such a choice also keeps the area of the forbidden regions small, which eases picking valid refinement points. On the other hand, such a choice may lead to an overly dense mesh. Raising the value of γ 0 results in a larger number of inconsistencies to be removed by the Inconsistency rule, an increasing number of calls to procedure Pick_valid, which therefore slows down the meshing process. In all our experiments, γ 0 is chosen between 1.3 and 1.8. Also, for efficiency reasons, we approximate γ(τ ) by the maximum distortion between two vertices of τ . The influence of parameter γ 0 on the resulting mesh density and algorithm efficiency is shown in Figure 7. In summary, our implementation depends on parameters β, δ, ρ 0 , γ 0 and r 0 . In all our experiments, β, δ, ρ 0 have been set to fixed values.

Computing the Metric Field

The algorithm requires to compute the metric field at any refinement point on S. It can either be prescribed by the user or approximate the curvature tensor of the surface at the refinement point. We detail below several cases to illustrate the versatility of our approach.

Three-Dimensional Scalar Fields

Three-dimensional scalar fields are common in physical or mathematical applications. If E denotes such a scalar field defined in R 3 and S a surface, we may want to mesh S so as to provide the best approximation of E on S. The triangles should then be elongated orthogonally to the gradient ∇E. We define the metric field as follows.

F M = U • diag{1/(1 + ϕ ∇E ), 1, 1} • U t where U = [ ∇E ∇E , U 1 , U 2 ],
U 1 and U 2 being two arbitrary unit vectors that form an orthogonal frame with ∇E. If the gradient is zero, the metric field is isotropic. Parameter ϕ controls how the facets are stretched with respect to the norm of the gradient.

Figure 4 presents an anisotropic Delaunay mesh of a sphere endowed with a sine-shaped shock scalar field. The triangles are elongated orthogonally to the gradient direction and the density depends both on the norm of the gradient and on the changing rate of the gradient.

Curvature Tensor on an Implicit Surface

Let S be an implicit surface defined by f (x, y, z) = 0. At each point p of S, we denote by N the normal vector -∇f / ∇f and by H the Hessian of f . The principal curvatures c max , c min are the non-zero eigenvalues of the matrix [START_REF] Hughes | Differential geometry of implicit surfaces in 3-space-a primer[END_REF]. Let U max , U min and N be the normalized eigenvectors of C, and let U = [U max , U min , N ]. The metric at point p is defined by M p = F t p F p , where:

C = P • H ′ • P where H ′ = H/ ∇f , P = I 3 -N • N t , see e.g.
F p = U • ∆ • U t ,
with ∆ = diag{e max , e min , e n }, e max = max{ǫ, c max }, e min = max{ǫ, c min }. Here, ǫ is a small positive constant that prevents the eigenvalues of the metric to vanish when c max and/or c min vanish, e.g., at planar or parabolic points. e n is a global constant to be discussed in the next subsection, and ∆ = diag{e max , e min , e n } stands for a diagonal matrix whose entries are {e max , e min , e n }.

Inria

(x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 . In
Figure 5 shows an anisotropic Delaunay mesh of the "chair" implicit surface.

Curvature Tensor on Polyhedral Surfaces

To remesh a polyhedral surface S that may be an approximation of a smooth surface S, we can estimate the curvature tensor of S from S. This is done in two steps.

1. Curvature estimation at the vertices of S We first estimate the curvature tensor at each vertex v ∈ S as follows. We compute the Euclidean k-nearest vertices of v, and then apply the jet-fitting algorithm of [START_REF] Cazals | Estimating differential quantities using polynomial fitting of osculating jets[END_REF] to Figure 6: Result on the polyhedral surface "fertility", a model from AIM@SHAPE repository. We took r 0 = 0.1, γ 0 = 1.5. The mesh has 12,480 vertices and 24,972 facets.

estimate the two principal curvatures c max and c min , the corresponding unit vectors U max and U min , and the normal vector N at v. From the estimated curvatures, we compute e max and e min as above. We choose e n as explained in subsection 5.4 and store at each vertex v a resulting estimated metric tensor M (v) = U ∆ 2 U t , with U = [U max , U min , N ] and ∆ = diag{e max , e min , e n }.

Blending metric tensors

To estimate the metric tensor at any prescribed query point on S, we smoothly blend the metric tensors computed in the previous step using the following formula :

M p = v∈S ω p (v)M (v), with ω p (v) = µ p (v) v∈S µ p (v), µ p (v) = e -p-v 2 /r 2 B
where r B is a parameter that controls the locality of the blending. By default, r B can be computed by averaging the distance to the k-nearest points of each vertex in S.

M p is a positive combination of symmetric tensors and therefore has no negative eigenvalues. This property can be proved by contradiction. Assume that M p has some negative eigenvalues. Then there exists a vector x that satisfies

x t M p x = v∈S ω p (v)x t M (v)x < 0.
Since all M (v) are symmetric tensors, x t M (v)x ≥ 0 and, since according to the definition, ω p (v) ≥ 0, we have x t M p x ≥ 0, which provides a contradiction.

Figure 6 shows an anisotropic Delaunay mesh of a polyhedral surface conforming to the estimated curvature metric field.

Computing Surface Stars

Our algorithm needs to maintain the surface star S v of each vertex v in the mesh. As described above, S v can be extracted from the star T v of v in Del v (V ). Yet, maintaining stars T v is not a fully satisfactory solution. Indeed, since the vertices of the mesh all lie on the surface, most of the tetrahedra in stars T v spread through space and their diameters do not decrease significantly when the sample density increases. As a consequence, the algorithm is not local and the computing time is not optimized. To remedy this problem, we add a few vertices close to the medial axis of the surface and maintain for each vertex v a triangulation T ′ v that does not necessarily include the star T v but nevertheless has the same restricted star S v as T v . The vertices close to the medial axis are computed from the initial points following Amenta and Bern's poles method [START_REF] Amenta | Surface reconstruction by voronoi filtering[END_REF]. The triangulation T ′ v is defined by the following insertion rule: a site p is inserted in the triangulation T ′ v of a site v only if p belongs to the so-called conflict zone of S v . The conflict zone of a surface star S v , denoted by Z v , is the union of the surface M v -circumballs B v (τ ) for all the facets τ of S v . The correctness of this variant of the algorithm follows from the fact that the conflict zone of star S v can only shrink when inserting a new site p, and from the fact that the surface star S v has to be updated iff p belongs to Z v .

When the sample density increases, the mesh facets become small and have small surface circumballs. At each insertion of a new site p, we use the conflict zones as filters to quickly report the vertices v whose small triangulations T ′ v have to be updated and the vertices w that have to be inserted in the triangulation T ′ p created for the new site. To each conflict zone Z v , we associate its (Euclidean) axis-parallel bounding box that we simply call the bounding box of v and denote by B v . We maintain in a data structure Q the bounding boxes of the current set of vertices V and we maintain V in a kd-tree Kd(V ). When inserting a new site p, we first query the data structure Q and report all the boxes that contain p. If p belongs to some box B v , we further check whether p belongs to the conflict zone Z v and, Inria in the affirmative, we insert p in the triangulations T ′ v . We then build the triangulation T ′ p of the new site. T ′ p is initialized as the star of p in the 3D Delaunay triangulation of p and of all the vertices v whose triangulations T ′ v have been updated by inserting p. Then, we compute a bounding box for the current conflict zone Z p , and query the kd-tree Kd(V ) to compute all the vertices in this box. Each vertex in this box is inserted in T ′ p if it belongs to the current conflict zone Z p .

Sliver Prevention

When the sampling is dense enough, quasi-cospherical configurations are tetrahedra with wellshaped facets (because of Rule 2) that have small volumes and are almost tangent to the surface. We call such tetrahedra slivers. We can reduce the number of slivers by enlarging the eigenvalue of the metric in the direction normal to the surface (eigenvalue e n defined in the preceding subsection). Indeed, the volume of those tetrahedra in the new metric will be increased, hence reducing the number of slivers and inconsistencies. Observe that this modification of the metric does not significantly affect the shape of the facets of the mesh since those facets are almost tangent to the surface. Experimental evidence indicates that e n should be of the same order of magnitude as the global maximum of c max .

Results

Performance and Complexity

Our experiments reveal that our algorithm is roughly linear behavior with respect to the number of vertices of the final mesh. As is well-known, for most point sets and insertion orders, inserting a new point in a Delaunay triangulation takes a constant time. Moreover, our use of the conflict zones (Section 5.3) ensures that a constant number of stars are visited when inserting a new vertex and that a constant number of vertices are considered when creating the star of the new vertex.

Figure 7 shows computation times for meshing a torus defined by an implicit function. All meshes are generated with r 0 = 0.1 and for different values of γ 0 . Each curve in Figure 7 consists roughly of two linear portions with different slopes. This corresponds to the fact that, in a first stage, only Rule 1 is applied. During the second stage, Rules 2 and 3 and the Pick_valid function are activated, which slows down the selection of new points.

Mesh Quality

Our meshing algorithm provides surface meshes of guaranteed quality since the algorithm controls the radius-edge ratio of each facet, as measured by the metric of any facet vertex. Figure 8 shows the histogram of the squared radius-edge ratios of the facets. For each facet, three radius-edge ratios were computed, one for each of the metrics of the vertices of the facet.

The quality of the meshes produced by our algorithm can still be improved by locally optimizing the positions of the vertices. This is done by adapting the optimization method of Alliez et al. [START_REF] Alliez | Variational tetrahedral meshing[END_REF]. Each vertex v is moved to the barycenter of the circumcenters of the facets in star S v , weighted by their respective areas. In our anisotropic setting, the circumcenters are M v -circumcenters, the weights are the M v -areas of the facets of S v , and the barycenter b v is an M v -barycenter. In order to preserve the locally uniform property of the anisotropic mesh, the motion should not introduce inconsistent, badly-shaped nor over-distortioned facet. Hence, instead of directly moving v to b v , we test new positions on the line segment (v, b v ) and check

v ′ = b v . If v ′ is a valid refinement point, we move v to v ′ , otherwise, we set v ′ = (v + v ′ )/2 while d v (v, v ′ ) > h.
Each star is optimized in turn. In our implementation, the stars are sorted by decreasing values of d v (v, b v ) and maintained in a priority queue. The first element is optimized first. The algorithm terminates when no more star can be optimized. Figure 8 shows that Figure 8: The distributions of the squared radius-edge ratios of the mesh facets, respectively before (blue curve) and after (red curve) the optimization step. The radius-edge ratio of each mesh facet is computed three times, using each time the metric attached to one of the three vertices. Model : "fertility", r 0 = 1.0, γ 0 = 1.5.

the distribution of radius-edge ratios is improved after the optimization step. The additional computing time required for this optimization is negligible.

Accuracy-Size Trade-off

Using an anisotropic metric field derived from the curvature tensor leads to surface meshes that enhance the trade-off between the accuracy of the approximation and the size of the mesh.

To measure the accuracy of a mesh, we use the mean error which we compute by uniformly sampling the mesh and averaging the Euclidean distance between each sample point on the mesh and its projection on the surface.

Table 1 provides a comparison between isotropic meshes and curvature adapted anisotropic meshes of tori. Each torus has a small radius r = 1, while the big radius R takes different values. In both the isotropic and the anisotropic cases, the meshes have been obtained by a refinement process where the size criterion (r 0 ) has been replaced by a bound on the quality of the approximation over a facet. Hence, both meshes are guaranteed to approximate the surface up to the same error bound. The bound on the approximation error is 0.001. Each column of the table gives the number of vertices in the isotropic mesh of a piece of the torus, the number of vertices in the anisotropic mesh of the same piece of the torus, and the ratio κ between these two numbers. Anisotropic meshes have been computed using the following parameter values ρ 0 = 3.0, γ 0 = 1.8, β = 2.5, δ = 0.3 and r 0 = 1.0. Figure 9 shows a piece of the mesh of the torus with R = 50 that was generated by our algorithm and highlights the high anisotropy ratio that can be reached.

As can be seen, the number of vertices N a of the anisotropic mesh can be much smaller than the number of vertices N i of the corresponding isotropic mesh with the same accuracy. The table also highlights the fact that the ratio between N a and N i decreases with the ratio r/R.

Conclusion

We have introduced a novel approach to generate anisotropic meshes of surfaces. Given a surface S endowed with a metric tensor that varies over S, our algorithm produces a mesh that approximates S and whose elements conform to the metric field. The resulting mesh provides a better approximation of the surface (if one takes as metric tensor the curvature tensor of S) or a better approximation of a function f defined over S (if we take as metric tensor the hessian of f ). Our algorithm is simple and comes with theoretical guarantees without relying on heuristics. As demonstrated by our experiments, it can handle complex shapes and metric fields. The algorithm is generic in the sense that it can handle various types of surfaces such as implicit surfaces or triangulated surfaces. Although we have not presented such results, it could be easily extended to mesh isosurfaces in 3D images.

As future work, we intend to consider surfaces with sharp features following the approach of [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF]. We also intend to extend our method to evolving surfaces [START_REF] Li | Generating well-shaped delaunay meshed in 3d[END_REF], and to 3D domains bounded by surfaces and equipped with a 3D metric field, combining the results of this paper with [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF].

More on parameters β and δ.

We provide here a few additional experimental results, to support our claim that parameters β and δ have little influence on the quality of the final mesh.

Table 2 shows computing times and numbers of vertices needed to mesh a torus (r = 1, R = 10) with r 0 = 0.1, γ 0 = 1.5 and different values of β (2,3,4,5) and δ (0.1, 0.3, 0.5). As can be Inria We took ρ 0 = 3.0, γ 0 = 1.8, r 0 = 1.0, β = 2.5, δ = 0.3. The approximation error is smaller than 0.001. seen, changing the value of β has basically no effect. The situation is different with parameter δ. As expected, choosing a small δ induces small picking regions and makes the finding of a valid refinement point harder. Hence, the picking process fails more often, leading to more vertex insertion trials and to an increased computing time.

Figure 10 shows the quality of the meshes obtained in the above experiment. Specifically, the figure shows the distribution of the squared radius-edge ratios of the mesh facets as in Fig. 8. We have selected three histograms to highlight the fact that the overall quality of the mesh is not impacted by varying β nor δ. The first histogram corresponds to our default values (δ = 0.3, β = 2). The two others correspond to one default value and one "large" value, specifically (δ = 0.3, β = 5), and (β = 2, δ = 0.5).

Figure 1 :

 1 Figure 1: Example of inconsistent stars in 2D : stars S v and S w are inconsistent because edge [vw] belongs to S v but not to S w .

Figure 3 :

 3 Figure 3: The picking region of triangle τ . Set σ = {q, r, s} is a hitting set for the picking region P v (τ ). It defines a forbidden region, shown with red dashes, to be avoided by the refinement point p of τ .

  It follows from the lemma that the M p -balls B p (p, l(p)), where l(p) = C 2(1+C)Γ sf(p), have disjoint interiors. A simple volume calculation then provides an upper bound on the number of vertices of the final mesh. Lemma 5 (Number of inserted sites) If the conditions of Lemma 3 are satisfied, the number n of sites inserted by the algorithm is at most

Figure 4 :

 4 Figure 4: A unit sphere endowed with a metric field derived from the sine-shaped shock function E(x, y, z) = tanh( 1 λ (2x -sin(5y))) + x 3 + xy 2 . λ = 0.6, ρ 0 = 3.0, r 0 = 0.1, γ 0 = 1.4, β = 2.5 and δ = 0.3. The mesh has 18323 vertices and 36642 facets.

Figure 5 :

 5 Figure 5: The "chair" implicit surface with a curvature-related metric field. The equation of the chair is f(x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 .In this

  Figure 5: The "chair" implicit surface with a curvature-related metric field. The equation of the chair isf (x, y, z) = x 2 + y 2 + z 2 -ak 2 2 -b (z -k) 2 -2x 2 (z + k) 2 -2y 2 . In this example, a = 0.8, b = 0.4, k = 1.0, ρ 0 = 3.0, r 0 = 0.1, γ 0 = 1.3, β = 2.5 and δ = 0.3. The mesh has 24164 vertices and 48336 facets.

Figure 7 :

 7 Figure 7: Computing time to mesh an implicit torus with r 0 = 0.1 and different values of γ 0 between 1.3 and 1.8.

Figure 9 :

 9 Figure 9: Anisotropic Delaunay mesh of a piece of torus with big radius 50 and small radius 1.We took ρ 0 = 3.0, γ 0 = 1.8, r 0 = 1.0, β = 2.5, δ = 0.3. The approximation error is smaller than 0.001.

Table 2 :

 2 Number of vertices and timings to mesh a torus with different values of β and δ.

	number of vertices β = 2	β = 3	β = 4	β = 5
	time
	δ = 0.1	5687	5652	5660	5643
		511.5 s 197.2 s 254.3 s 209.1 s
	δ = 0.3	5227	5104	5130	5065
		93.5 s 101.5 s 121.2 s 114.7 s
	δ = 0.5	4926	4979	4903	4973
		72.9 s	77.4 s	81.0 s	72.4 s

RR n°8400

Inria

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

This work has been partially supported by the Agence Nationale de la Recherche (ANR) under project GIGA (Geometric Inference and Geometric Approximation). The second author was partially supported by Chinese 973 Program (2010CB328001), Chinese 863 Program (2012AA040902) and the NSFC (61035002, 61272235) * Geometrica