Jean-Daniel Boissonnat

Arijit Ghosh

Triangulating Smooth Submanifolds with Light Scaffolding

Keywords: Algorithms, Certification, and Programming, Software and Architecture Manifold triangulation, meshing, manifold learning, manifold sampling, computational geometry, computational topology

We propose an algorithm to sample and mesh a k-submanifold M of positive reach embedded in R d . The algorithm first constructs a crude sample of M. It then refines the sample according to a prescribed parameter ε, and builds a mesh that approximates M. Differently from most algorithms that have been developped for meshing surfaces of R 3 , the refinement phase does not rely on a subdivision of R d (such as a grid or a triangulation of the sample points) since the size of such scaffoldings depends exponentially on the ambient dimension d. Instead, we only compute local stars consisting of k-dimensional simplices around each sample point. By refining the sample, we can ensure that all stars become coherent leading to a k-dimensional triangulated manifold M. The algorithm uses only simple numerical operations. We show that the size of the sample is O(ε -k ) and that M is a good triangulation of M. More specifically, we show that M and M are isotopic, that their Hausdorff distance is O(ε 2 ) and that the maximum angle between their tangent bundles is O(ε). The asymptotic complexity of the algorithm is T (ε) = O(ε -k 2 -k ) (for fixed M, d and k).

Triangulation efficace de variétés lisses

Résumé : On propose un algorithme pour échantillonner et mailler une sousvariété M de dimension k plongée dans R d . Après avoir construit un échantillon grossier, l'algorithme raffine l'échantillon et le maillage selon un paramètre ε. L'algorithme ne construit pas de subdivision de R d mais seulement des triangulations locales (stars) de dimension k autour de chaque point de l'échantillon. On montre qu'en raffinant l'échantillon, on peut rendre toutes les stars cohérentes et ainsi obtenir une variété triangulée M qui approche M. L'algorithme n'utilise que des opérations numériques simples, la taille de l'échantillon produit est O(ε -k ). On montre que M et M ont le même type topologique, que leur distance de Hausdorff est O(ε 2 ) et que l'angle entre leurs espaces tangents est O(ε). La complexité asymptotique de l'algorithm est T (ε) = O(ε -k 2 -k ) (pour M, d k fixés).

Mots-clés :

Triangulation de variétés, génération de maillage, échantillonnage, géométrie algorithmique, topologie algorithmique

Introduction

We intend to sample and mesh a k-manifold M of positive reach embedded in R d . Manifolds of positive reach have been introduced by Federer [START_REF] Federer | Curvature Measures[END_REF][START_REF] Federer | Geometric measure theory[END_REF] and include in particular C 2 manifolds. By mesh, we mean an embedded polyhedral approximation of M made up of simplices. We are especially interested in the case where the dimension k of M is much smaller than d, and intend to design an algorithm whose complexity depends on k rather than on d. Applications can be found in scientific computing for solving partial differential equations where the domain of interest has the structure of a manifold, in dynamical systems for computing the topology of space attractors, and in statistics and machine learning to approximate statistical manifolds.

Related work

The problem of triangulating manifolds has a long history in the mathematical literature. In differential topology, seminal contributions are due to Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF], Cairns [START_REF] Cairns | A simple triangulation method for smooth manifolds[END_REF], Munkres [START_REF] Munkres | Elementary Differential Topology[END_REF], Whitehead [START_REF] Whitehead | On C 1 -complexes[END_REF] to name a few. Although these papers are not of an algorithmic nature, they introduce and study several interesting concepts that have been extensively used in Computational Geometry recently such as Voronoi diagrams restricted to a manifold, ε-sample of a manifold, fat (or thick) triangulations. However, these papers do not discuss the geometric quality of the approximation nor the size of the sample. The optimal sampling and approximation of convex bodies is also a long standing problem in convex optimization with major contributions by Gruber [START_REF] Gruber | Asymptotic estimates for best and stepwise approximation of convex bodies I[END_REF][START_REF] Gruber | Optimum quantization and its applications[END_REF] and Dudley [START_REF] Dudley | Metric entropy of some classes of sets with differentiable boundaries[END_REF]. Recently, Clarkson [START_REF] Clarkson | Building triangulations using ε-nets[END_REF] extended this line of work to non-convex smooth manifolds of arbitrary dimensions. However, his algorithm follows an intrinsic point of view which makes it difficult to use in practice since it requires to compute geodesic distances on the manifold which may be quite complicated in practice [START_REF] Peyré | Geodesic computations for fast and accurate surface remeshing and parameterization[END_REF]. Other, more practical algorithms for approximating convex bodies, including the well-known sandwich algorithm, have been analyzed by Kamenev [START_REF] Kamenev | The initial convergence rate of adaptive methods for polyhedral approximation of convex bodies[END_REF]. We are not aware of similar studies for non convex manifolds except for the case of surfaces embedded in R 3 which has been extensively studied in the Computational Geometry literature. See [START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF] for a recent survey. These methods start by computing some subdivision of the embedding space (such as a grid or a triangulation of the sample points) and their direct extension to higher dimensions would face an exponential dependence on d. A step in this direction is the extension of the celebrated Marching Cube algorithm to manifolds of higher dimensions [START_REF] Min | Simplicial Isosurfacing in Arbitrary Dimension and Codimension[END_REF][START_REF] Bhaniramka | Isosurface construction in any dimension using convex hulls[END_REF]. Continuation methods do not use any subdivision of the ambient space and are close in spirit to our approach. They construct a triangulated approximation of a k-dimensional submanifold in a greedy way and extend the current k-dimensional triangulated domain by adding a neighborhood of a boundary point. Some experimental results can be found in [START_REF] Henderson | Multiple parameter continuation: computing implicitly defined k-manifolds[END_REF] but no theoretical analysis of continuation methods is available.

Our approach

In this paper, we will follow the extrinsic approach but show that we can avoid using any d-dimensional data structure (except in the initialization step). This paper extends a technique developped for anisotropic mesh generation [START_REF] Boissonnat | Locally Uniform Anisotropic Meshing[END_REF] and builds on results from our companion paper investigating the related problem of manifold reconstruction [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF].

We assume that the manifold M to be meshed has a positive reach and that we know a lower bound on the reach. In addition we assume that we can compute, for each point p ∈ M, the k-dimensional tangent space T p of M at p.

The algorithm starts with a sufficiently dense sample of M and then refines the sample and builds a mesh that approximates M so as to satisfy a prescribed sampling rate ε. The size of the initial sample does not depend on ε but only on M. For each sample point p ∈ P, we compute its k-dimensional star in the restriction of the d-dimensional Delaunay triangulation of the sample P to the tangent space T p at p. Such a star can be computed in the k-dimensional flat T p once we have projected P onto T p .

In general, the stars do not glue coherently and it may well happen that q is a vertex in the star of p while p is not a vertex in the star of q. The crucial observation is that by refining the sample, we can ensure that all the stars become coherent leading to a k-dimensional mesh M. For ε small enough, we show that the size of the sample is O(ε -k ) and that M is a good approximation of M. Specifically, we show that M and M are isotopic, that their Hausdorff distance is O(ε 2 ) and that the maximum angle between their tangent bundles is O(ε). Our bound on the Hausdorff distance matches the lower bound of Clarkson [START_REF] Clarkson | Building triangulations using ε-nets[END_REF] (up to a multiplicative constant that depends on M). The bound on the distance between the tangent bundles seems to be new.

To refine the mesh according to a sampling parameter ε, we need an oracle to query the manifold and to compute new points on M. This is a critical issue with respect to practical efficiency. In our algorithm, we only need to compute a point in the (0-dimensional) intersection of M with a (dk)-flat. Except from the oracle and the projection of points onto k-dimensional flats (the tangent spaces at the points of P), all computations are performed in those k-flats. As a consequence, the asymptotic complexity of the algorithm is O(ε -k 2 -k ) for fixed k, d, and M. Hence, while our approach is extrinsic, the ambient dimension appears only in the constant hidden in the big-O.

The present work combines four main ideas that have been introduced separately before : the general mechanism of Delaunay refinement [START_REF] Chew | Guaranteed-Quality Delaunay Meshing in 3D[END_REF][START_REF] Ruppert | A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation[END_REF][START_REF] Boissonnat | Provably Good Sampling and Meshing of Surfaces[END_REF], the concept and properties of Delaunay triangulations restricted to a manifold [START_REF] Amenta | A Simple Algorithm for Homeomorphic Surface Reconstruction[END_REF][START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF][START_REF] Edelsbrunner | Geometry and Topology for Mesh Generation[END_REF], the notion of tangential Delaunay complex [START_REF] Boissonnat | A coordinate system associated with points scattered on a surface[END_REF][START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF], and a perturbation technique due to Li to remove flat simplices [START_REF] Li | Generating Well-Shaped d-dimensional Delaunay Meshes[END_REF]. Several of the structural results we need are borrowed from [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF]. However, the algorithm in [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF] takes as input a given set of points while here the algorithm has to construct the sample as well as the mesh, which makes the algorithm different and its analysis more delicate. This paper also aims at clarifying the basic operations that are required to triangulate a manifold.

Organization of the paper

We recall in Section 2 the definition of tangential Delaunay complex. This complex is embedded in R d but is not in general a k-dimensional triangulation due to the presence of so-called inconsistent configurations to be studied in Section 2. To remove inconsistent configurations, we propose an algorithm that refines the complex. The algorithm is described in Section 3 and analyzed in Section 4. Lastly, in Section 5, we show that the output of the algorithm is a good approximation of M.

Notations

In the paper, M denotes a compact closed k-manifold of positive reach embedded in R d and P a finite set of points on M. The tangent space at x ∈ M is denoted by T x and the normal space by N x . For a point p in R d and r ≥ 0, B(p, r) ( B(p, r)) denotes the d-dimensional open (closed) ball centered at p of radius r, and

B M (p, r) ( BM (p, r)) denotes B(p, r) ∩ M ( B(p, r) ∩ M).
If U and V are two affine spaces with dim U ≤ dim V . The angle between U and V is defined as

∠(U, V ) = max u∈U min v∈V ∠(u, v),
where u and v are vectors in U and V respectively.

Definitions and preliminaries

This section recalls some definitions and results borrowed from [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF]. For completeness, proofs are given in the appendix.

Sampling conditions

Let π : R d → M map each point of R d to its closest point on M. The reach of M, denoted by rch(M), is defined as the supremum of all δ such that any point x of R d lying at distance less than δ from M has a unique image π(x). In this paper, we assume that M has strictly positive reach.

As shown by Federer [START_REF] Federer | Curvature Measures[END_REF], rch(M) is (strictly) positive when M is of class C 2 or even C 1,1 , i.e. the normal bundle is defined everywhere on M and is Lipschitz continuous.

The reach is used for defining sampling conditions on the manifold.

A point sample P is said to be a t-sample if, for any point x ∈ M, there exists a point p ∈ P such that px ≤ t rch(M).

We call minimal interpoint distance of P the smallest distance between any two points of P. An ε-sample P is called κ-sparse if the minimal interpoint distance of P is at least κ rch(M). A κ-sparse t-sample of M is simply called a (t, κ)-sample.

We will use the following results [START_REF] Federer | Curvature Measures[END_REF][START_REF] Giesen | Shape dimension and intrinsic metric from samples of manifolds[END_REF].

Lemma 1 1. For any point q ∈ M such that pq = t rch(M) for some 0 < t < 1, sin ∠(pq, T p ) ≤ t/2.

2.

Let q be a point in T p such that pq = t rch(M) for some 0 < t ≤ 1/4.

Let q ′ be the point on M closest to q. Then qq ′ ≤ 2t pq .

Properties of Simplices

A j-dimensional simplex (or j-simplex for short) τ is the convex hull of j + 1 affinely independent points p 0 , . . . , p j . We write τ = [p 0 , . . . , p j ] and, for convenience, we may confound a simplex and the set of its vertices. We write c τ for the circumcenter of τ (i.e. the center of its minimum enclosing d-ball), aff(τ ) for the j-dimensional affine hull of τ , N τ for the (dj)-dimensional affine space normal to aff(τ ) and passing through c τ (which lies in aff(τ )).

For any j-simplex τ , we denote by r τ the circumradius of τ (i.e. the radius of its minimum enclosing d-ball), by L τ (∆ τ ) the length of its shortest (longest) edge, by ρ τ = r τ /L τ the radius-edge ratio of τ , and by vol(τ ) the j-dimensional volume of τ . For a vertex p of τ , we write τ p = τ \{p} for the (j -1)-dimensional face of τ opposite to p, and D p (τ ) for the distance from p to the affine hull aff(τ p ) of τ p . In addition, we define the fatness of τ as

Θ τ = 1 if j = 0 vol(τ )/∆ j τ if j > 0 (1)
The following lemma is proved in Appendix A.

Lemma 2 (Properties of simplices) Let τ = [p 0 , . . . , p j ] be a j-simplex and p be a vertex of τ .

1. Θ τ ≤ 1 j! . 2. j! Θ τ ≤ Dp(τ ) ∆τ ≤ j 2 j-1 ρ j-1 τ Θτ Θτ p . 3. The distance of p from the (j -1)-sphere ∂B(c τp , r τp ) ∩ aff(τ p ) is less than b(ρ τ ) D p (τ ), where b(ρ τ ) = 1 + 1 1- √ 1-1/4ρ 2 τ .
The following lemma is due to Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF].

Lemma 3 Let τ = [p 0 , . . . , p j ] be a j-dimensional simplex and let H be a jdimensional affine flat such that τ is contained in the offset of H by η (i.e. any point of τ is at distance at most η from H). If u is a unit vector in aff(τ ), then there exists a unit vector u H in H such that sin ∠(u, u H ) ≤ 2η (j-1)! Θτ Lτ .

We deduce from the above lemma the following important corollary that bounds the angle between a simplex and the tangent space at a vertex of the simplex. See also Lemma 1 in [START_REF] Fu | Convergence of curvatures in secant approximations[END_REF] and Lemma 16 in [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]. Proof. It suffices to apply Lemma 3 with H = T p and to use

η = ∆ 2 τ /2 rch(M) (from Lemma 1 (1)) and r τ /ρ τ = L τ ≤ ∆ τ ≤ 2 r τ . Hence sin ∠(T p , aff(τ )) ≤ 2η (j -1)!Θ τ L τ ≤ 2ρ τ ∆ τ Θ τ rch(M) ≤ 2ρ τ ∆ τ Θ τ rch(M) .

Tangential Delaunay complex and inconsistent configuration

Let P be a finite set of points on M and Del(P) be the d-dimensional Delaunay triangulation of P, i.e. the collection of all the simplices with vertices in P that admit an empty circumscribing d-dimensional ball. A ball (or more generally any domain of R d ) is called empty if its interior contains no point of P. Let in addition Del pi (P) be the Delaunay triangulation of P restricted to the tangent space T pi , i.e. the collection of all the simplices with vertices in P that admit an empty circumscribing d-dimensional ball centered on T pi . Equivalently, the simplices of Del pi (P) are the simplices of Del(P) whose Voronoi dual face intersect T pi . Observe that Del pi (P) is in general a k-dimensional complex and can always be made k-dimensional by applying some infinitesimal perturbation on P. We will assume that the points of P are in general position in the rest of the paper, meaning that all Del pi (P) are k-dimensional triangulations. Finally, write star(p i ) for the star of p i in Del pi (P), i.e. the set of simplices that are incident to p i in Del pi (P). See Figure 1.

We recall the definition of tangential Delaunay complex and some known results from [START_REF] Boissonnat | A coordinate system associated with points scattered on a surface[END_REF][START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF].

Definition 1 (Tangential Delaunay complex) We call tangential Delaunay complex the simplicial complex Del T M (P) = {τ, τ ∈ star(p), p ∈ P}.

RR n°7660

Plainly, Del T M (P) is a subcomplex of Del(P). The following easy lemma is crucial since it shows that computing the tangential Delaunay complex reduces to computing n weighted Delaunay triangulations in k-dimensional flats if n denotes the cardinality of P. See Appendix B for a proof.

We denote by π i : P → T pi the orthogonal projection of P onto T pi and by Π i : P → T pi × R the 1-1 mapping that associates to a point p the weighted point defined by

Π i (p) = (π i (p), -π i (p) -p 2 ).
Lemma 4 Del pi (P) is the pullback by Π i of the k-dimensional weighted Delaunay triangulation of Π i (P).

Let τ = [p 0 , . . . , p k ] be a k-simplex with vertices in M and let denote by B pi (τ ) the d-dimensional ball circumscribing τ that is centered on T pi . The corresponding center and radii are denoted by c pi (τ ) and r pi (τ ). The following lemma, which is a variant of Lemma 10 in [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF], bounds the size of the simplices of Del T M (P) as a function of the sampling density. See Appendix C for a proof.

Lemma 5 Let P be an ε-sample of a manifold M with ε ≤ 1/8. Then we have:

1. Vor(p) ∩ T p ⊆ B(p, 4ε rch(M)).
2. for any k-simplex τ ∈ star(p), r p (τ ) ≤ 4ε rch(M).

3. for all edges pq ∈ Del T M (P), pq ≤ 8ε rch(M).

In general, the tangential Delaunay complex is not a triangulated k-manifold. This is due to the presence of so-called inconsistent simplices. Refer to Figure 2.

Let τ be a k-simplex in the star of p i which is not in the star of p j . Equivalently, the Voronoi (dk)-dimensional face Vor(τ ) dual to τ intersects T pi (at a point c pi (τ )) but does not intersect T pj . Observe that c pi (τ ) is the center of an empty d-dimensional ball B pi (τ ) circumscribing τ . Let c pj (τ ) denote the intersection of aff(Vor(τ )) with T pj . Differently from B pi (τ ), the d-dimensional ball B pj (τ ) centered at c pj (τ ) that circumscribes τ contains a subset P j (τ ) of points of P in its interior. Accordingly, the line segment [c pi (τ ) c pj (τ )] intersects the interior of some Voronoi cells (in particular, the cells of the points of P j (τ )). We denote by p l the point of P \ τ whose Voronoi cell is hit first by the segment [c pi (τ ) c pj (τ )], when oriented from c pi (τ ) to c pj (τ ). We now formally define an inconsistent configuration.

c p i (τ ) T p j c p j (τ ) M i φ
Figure 2: The figure shows an example of an inconsistent simplex τ = [p i , p j ] that belongs to star(p i ) and does not belong to star(p j ). Vor(τ ) intersects T pi but not T pj .

Write i φ for the first point of Vor(p l ) hit by the oriented segment [c pi (τ ) c pj (τ )]. i φ is the center of an empty d-dimensional ball that circumscribes φ. Hence φ is a Delaunay (k + 1)-simplex and i φ is the point on [c pi (τ ) c pj (τ )] that belongs to Vor(φ), the Voronoi face dual to φ. Since we assumed that the points are in general position, an inconsistent configuration cannot belong to the tangential Delaunay complex (which does not contain faces of dimension greater than k). Observe also that some of the subfaces of an inconsistent configuration may not belong to the tangential Delaunay complex.

We will use the same notations for inconsistent configurations as for simplices, e.g. r φ and c φ for the circumradius and the circumcenter of φ, ρ φ and Θ φ for its radius-edge ratio and fatness respectively. We also write R φ = i φp i , where

p i is a vertex of φ. Note that R φ = i φ -p i ≥ c φ -p i = r φ .
The following important lemma bounds the radius and fatness of an inconsistent configuration.

Lemma 6 Let φ be an inconsistent configuration witnessed by p, q and r, and let τ = φ \ {r}. Assume that r φ < rch(M)/4, and write θ = max θ x where θ x = ∠(aff(τ ), T x ) and x is a vertex of τ (sin θ ≤ 4 ρτ rτ Θτ rch(M) by Corollary 1). We have

1. R φ ≤ rτ cos θ . 2. Θ φ ≤ rτ rch(M) cos θ 1 + 4ρτ Θτ .
Proof.

RR n°7660

(a) (b) 1. We have c p (τ ) = Vor(τ )∩T p , c q (τ ) = aff(Vor(τ ))∩T q , and r p (τ ) = c p (τ )-p and r q (τ ) = c q (τ )q . Since θ = max x θ x where θ x = ∠(aff(τ ), T x ) and x is a vertex of τ , we have r p ′ (τ ) ≤ r τ / cos θ and c p ′ (τ )c τ ≤ r τ tan θ, for p ′ ∈ {p, q}. As i φ ∈ [c p (τ ), c q (τ )], we have i φ -c τ ≤ r τ tan θ. Then, by Pythagoras theorem, we have

p q r p q r s v s v u u T p T q T r T s T v T u
R φ = r 2 τ + i φ -c τ 2 ≤ r τ √ 1 + tan 2 θ = r τ / cos θ.
2. We will now bound Θ φ = vol(φ) ∆ k+1 φ . We use vol(φ) = Dr (φ) vol(φr) k+1 and bound D r (φ) and vol(τ ).

Using the fact that ∆ φ ≤ 2r φ ≤ 2rτ cos θ from 1, we have

D r (φ) = dist(r, aff(τ )) = sin ∠(pr, aff(τ )) × p -r ≤ (sin ∠(pr, T p ) + sin ∠(aff(τ ), T p )) × ∆ φ ≤ p -r 2 rch(M) + 2ρ τ ∆ τ Θ τ rch(M) ∆ φ ≤ ∆ 2 φ 2 rch(M) 1 + 4ρ τ Θ τ (2)
From the definition of fatness of a simplex and Lemma 2 (1), we get

vol(τ ) = Θ τ ∆ k τ ≤ ∆ k τ k! . (3)

RR n°7660

Using inequalities ( 2) and ( 3), and

∆ τ ≤ ∆ φ ≤ 2r φ ≤ 2 rτ cos θ , we get Θ φ = vol(φ) ∆ k+1 φ = D r (φ) vol(τ ) k + 1 × 1 ∆ k+1 φ ≤ ∆ 2 φ 2 rch(M) 1 + 4ρ τ Θ τ × ∆ k τ (k + 1)! ∆ k+1 φ ≤ r τ rch(M) cos θ 1 + 4ρ τ Θ τ 3 Algorithm
We now describe our meshing algorithm. The algorithm assumes that we know the dimension k of M and that we can get the tangent space T p at any point p ∈ M. In addition, we assume to know a positive lower bound on the reach of the manifold M. We write it also rch(M) for simplicity.

The algorithm takes as input parameters ε, ρ 0 ≥ 1/2, Θ 0 < 1/2. The sampling parameter ε will be used in Section 5 to bound the size of the sample and the approximation error. The two constants ρ 0 and Θ 0 are used below to define good simplices and slivers (a kind of flat simplices). Additional parameters and conditions will be specified in Section 4.

The algorithm first constructs an initial sample P 0 of M of constant size. Then, it upsamples P 0 by inserting new points on M in a greedy way so as to satisfy a sampling condition expressed in terms of parameter ε, and making sure that all the stars are consistent.

We now detail the main features of the algorithm.

Primitive operations

We assume that the manifold is generic in the sense that the intersection of any (dk)-flat with the manifold is a bounded set of points. The only primitive operation of our algorithm that involves M, namely ints(M, F ), computes the intersection of M with a (dk)-flat F . This primitive operation can be implemented for various representations of M : e.g. when M is given implicitly as a system of dk algebraic equations, computing ints(M, F ) reduces to solving a 0-dimensional system of d-variate algebraic equations.

We will also need to pick random points in Euclidean balls of R k .

RR n°7660

3.2 Computing the initial sample P 0

The construction of the initial sample P 0 can be done in various ways. We can use the continuation method of [START_REF] Henderson | Multiple parameter continuation: computing implicitly defined k-manifolds[END_REF] or use a simpler grid. We sketch the grid method which is easy to implement although the construction requires 2 O(d log d) time. Take a uniform d-dimensional grid with cells of diameter rch(M)/16 and pick the intersection points between the manifold and the (dk)-faces of the grid to build a set S ⊂ M which is an 1/32-sample of M. To make the sample sparse, we do the following:

1. Set P 0 = ∅ and S = S;

2. Take a point p from S, insert p in P 0 , and remove from S the points that belong to B(p, rch(M)/32).

Repeat

Step 2 until S = ∅.

The subsample P 0 ⊆ S ⊂ M, is a 1/32-sparse 1/32-sample of S, which in turn is a 1/16-sample of M. Therefore P 0 will be a 1/32-sparse 1/16-sample of M.

Good simplices and slivers

We adapt the following definitions from [START_REF] Li | Generating Well-Shaped d-dimensional Delaunay Meshes[END_REF].

Definition 3 (Good simplex) A simplex τ is a good simplex if ρ τ ≤ ρ 0 and min σ⊆τ, dim(σ)>0 Θ 1 dim(σ) σ ≥ Θ 0 ,
where dim(σ) denotes the dimension of the simplex σ.

Definition 4 (Sliver

) A j-simplex τ is called a sliver if j > 1, ρ τ ≤ ρ 0 , Θ τ < Θ j
0 , and all of its proper subfaces are good simplices.

The next lemma follows from Lemma 6 (2). It relates inconsistent configurations and slivers.

Lemma 7 Let φ be an inconsistent configuration witnessed by p, q and r, and

let τ be the k-dimensional simplex φ \ {r}. Assume that ρ φ ≤ ρ 0 , r τ ≤ ε rch(M)
and that the subfaces of φ are good simplices. Then, if

ε ≤ Θ k+1 0 1 + 4ρ0 Θ k 0 2 + 16ρ 2 0 Θ 2 0
then φ is a sliver.

RR n°7660

Proof.

From Lemma 6 (2), we have

Θ φ ≤ rτ rch(M) cos θ 1 + 4ρτ
Θτ where sin θ ≤ 4ρτ rτ Θτ rch(M) . Using the fact that r τ ≤ r p (τ ) ≤ εrch(M), ρ τ ≤ ρ 0 , and

Θ τ ≥ Θ k 0 , we deduce that, if ε ≤ Θ k+1 0 1 + 4ρ0 Θ k 0 2 + 16ρ 2 0 Θ 2 0 then Θ φ < Θ k+1 0 . The lemma follows.
Hence, if ε is small enough, removing all slivers of dimensions at most k + 1 will result in removing inconsistencies from Del T M (P).

This remark motivates the following definition. The completed complex will be maintained by the algorithm and slivers will be removed from this complex.

Definition 5 (Completed complex) The completed complex C(P) consists of the following simplices and their subfaces:

(i) The k-simplices of Del T M (P).
(ii) The inconsistent configurations φ witnessed by p, q and r, such that 1.

r p (τ ) ≤ ε rch(M) and 2. τ = φ \ {p l } is a good simplex.

Picking region and good points

A new point to be inserted is chosen so as to remove a bad simplex σ of C(P).

It will be taken from the so-called picking region of σ which we define now. We introduce two new parameters, β > 1 and δ ∈ [0, 1).

Definition 6 (Picking region Π(σ, δ)) We consider the following two cases:

1. If σ = τ is a k-dimensional simplex in star(p), then the picking region of τ is defined as Π(τ, δ) = B(c p (τ ), δ r p (τ )) ∩ M. 2. If σ = φ is an inconsistent configuration, then the picking region of φ is defined as Π(φ, δ) = B(i φ , δ R φ ) ∩ M.
Definition 7 (Tiny sliver) A simplex τ is called a tiny sliver with respect to a simplex σ if τ is a sliver and r τ ≤ β r σ .

Definition 8 (Good point) A point x in a picking region Π(σ, δ) is called a good point if inserting x does not create any j-dimensional sliver that is both incident to x and tiny with respect to σ, j ≤ k + 1.

The algorithm makes use of the following two functions: Observe that S3 prevents to create simplices incident to x that are tiny with respect to τ .

RR n°7660
If σ is an inconsistent configuration φ, we proceed as follows. Let φ be witnessed by p, q and r. According to the definition of an inconsistent configuration, the k-dimensional simplex τ = φ \ {r} belongs to star(p) and not to star(q). We implement good-pick(φ, δ) as in Case 1 except that we pick random points from the k-dimensional ball B(c p (τ ), r) ∩ T p where r = δ R φ + i φc p (τ ) .

In Section 4, we will prove the existence of good points in Π(σ, δ).

Refinement Algorithm

We can now give the details of the algorithm.

input a finite 1/32-sparse 1/16-sample P 0 of M, and parameters ε, ρ 0 , Θ 0 , β and δ (the parameters should be chosen so that they satisfy the conditions given in Theorem 1 in Section 4)

output a sample P of M and M = Del T M (P)

The refinement algorithm consists of applying the following rules. Rule (i) is only applied if Rule (j) with j < i cannot be applied. Each rule kills a simplex σ (i.e. removes σ from a star) by inserting a new point in its picking region. To insert (or remove) a point means here to update P, as well as the completed complex C(P). We call new simplex a simplex of C(P) that is created when inserting a new point.

Notice that all the new k-simplices in Del T M (P) and all the new (k+1)-simplices in C(P) will be incident to the newly inserted point p. Observe however that a simplex (possibly a sliver) that existed in the Delaunay triangulation Del(P) but not in C(P) before the insertion of p may become a (subface of a new) simplex of C(P) after the insertion of p.

Rule 1 Big simplices : if there exists a k-simplex τ in star(p)

s.t. r p (τ ) > ε rch(M), insert x = pick(c p (τ ), p).
Rule 2 Bad radius-edge ratio :

a If there exists a k-simplex τ in star(p) such that r τ > ρ 0 L τ , insert x = pick(c p (τ ), p).
b Similarly, if φ is an inconsistent configuration witnessed by p, q and l, such that r φ > ρ 0 L φ , insert x = pick(i φ , p).

Rule 3 Type-1 sliver : If there exists a k-simplex τ of Del T M (P) that is a sliver or has a subsimplex that is a sliver, insert x = good-pick(τ, δ).

Rule 4 Type-2 sliver : If an inconsistent configuration φ ∈ C(P) is a sliver or has a subsimplex that is a sliver, insert x =good-pick(φ, δ).

Once the algorithm terminates, all slivers and inconsistent configurations have been killed. Hence, all stars are consistent and a simple sweep allows to merge all the stars into the final mesh M = Del T M (P).

Analysis of the algorithm

To prove that the algorithm terminates, we first bound the volume of the socalled forbidden regions. This will be helpful in proving that there exist good points in the picking regions. Termination of the algorithm is then proved by showing that the interpoint distance remains bounded from below. Lastly, we analyze the time complexity of the algorithm.

The following lemma is proved in Appendix D.

Lemma 8 Let p be a point on M. There exist constants ξ and A that depend only on k such that, for all t ≤ ξ and r = t rch(M), we have

0 < 1 -A t ≤ vol(B M (p, r)) φ k r k ≤ 1 + A t
where φ k is the volume of the k-dimensional unit Euclidean ball.

RR n°7660

Forbidden regions

For a given j-simplex (1 ≤ j ≤ k) with vertices on M, the forbidden region F µ of µ is defined as

F µ = {x ∈ M | µ ∪ {x} forms a (j + 1)-dimensional sliver}.
Remember that µ must be a good simplex by definition of a sliver. We will now bound the volume of F µ .

Lemma 9 Let µ be a good j-dimensional simplex with 2 ≤ j ≤ k with vertices on M, r µ ≤ t rch(M).

If (i) t ≤ ξ/2 (ξ is defined in Lemma 8), (ii) 4ρ0 Θ k 0 + 2 t < Θ 0 and (iii) (B + 1)Θ 0 ≤ 1 for some B that depends on k and ρ 0 , then vol(F µ ) ≤ D Θ 0 r k µ
, where D depends also on k and ρ 0 .

Proof. Let x ∈ F µ and x * be the point closest to x on ∂B(c µ , r µ ) ∩ aff(µ). We denote by τ the (j + 1)-dimensional simplex τ = µ ∪ {x}, and τ is a (j + 1)dimensional sliver since x ∈ F µ . From Lemma 2 (2) and ( 3), the facts that Let p be a vertex of µ. Let c be the point closest to c µ on T p and c * be the point closest to c on M (see Figure 4). From Corollary 1, we have

∆ τ ≤ 2r τ ≤ 2ρ τ L τ ≤ 2ρ τ L µ ≤ 4ρ τ r µ , ρ τ ≤ ρ 0 and Θτ Θµ < Θ 0 (since τ is a (j + 1)-dimensional sliver), we get x -x * ≤ b(ρ τ ) dist(x, aff(µ)) ≤ (j + 1)2 j ρ j τ b(ρ τ ) Θ τ Θ µ ∆(τ ) ≤ (j + 1)2 j+2 ρ j+1 0 b(ρ 0 ) Θ 0 r µ ≤ (k + 1)2 k+2 max(1, ρ 0 ) k+1 b(ρ 0 ) Θ 0 r µ def = BΘ 0 r µ (4) p T
c -c µ ≤ sin ∠(T p , aff(µ)) × r µ ≤ 4ρ 0 t Θ j 0 × r µ ≤ 4ρ 0 t Θ k 0 × r µ def = C t r µ , (5)

RR n°7660

the last inequality follows from the fact that Θ 0 < 1. From Lemma 1 (2) we have

c -c * ≤ 2 c -p 2 rch(M) ≤ 2 r 2 µ rch(M) ≤ 2t r µ . (6)
Using the fact that c µx * = r µ and inequalities (4), ( 5) and ( 6), we get

c * -x ≤ c * -c + c -c µ + c µ -x * + x * -x ≤ r µ (1 + (BΘ 0 + (C + 2) t)) < r µ (1 + (B + 1)Θ 0 ) ,
the last inequality follows from hypothesis (ii), which implies (C + 2)t ≤ Θ 0 . We can similarly prove that cx ≥ r(1 -(B + 1)Θ 0 ).

Writing δ = (B + 1)Θ 0 ≤ 1 (from hypothesis (iii)), we deduce from the inequalities above that c

* -x ∈ [r µ (1 -δ), r µ (1 + δ)]. Therefore, the forbidden region F µ is included in B M (c * , r µ (1 + δ)) \ B M (c * , r µ (1 -δ))
. We now use Lemma 8 to bound the volume of F µ . Observe that Lemma 8 can be applied since

r µ (1 + δ) ≤ 2r µ ≤ 2t rch(M) ≤ ξ rch(M) (as t ≤ ξ/2 and δ ≤ 1). We have vol(F µ ) φ k ≤ vol(B M (c * , r µ (1 + δ)) \ B M (c * , r µ (1 -δ))) φ k ≤ (1 + A(1 + δ)t) r k µ (1 + δ) k -(1 -A(1 -δ)t) r k µ (1 -δ) k = r k µ (1 + δ) k -(1 -δ) k + A t r k µ (1 + δ) k+1 + (1 -δ) k+1 ≤ 2 k δ r k µ + A(2 k+1 + 1) t r k µ ( 7
)
the last inequality follows from the fact that

(1 + x) k -(1 -x) k ≤ 2 k x for x ∈ [0, 1].
From hypothesis (ii) and the fact that Θ 0 < 1, we have

t < Θ 0 4ρ0 Θ k 0 + 2 ≤ Θ 0 4ρ 0 + 2 .
Using this inequality and inequality [START_REF] Boissonnat | Provably Good Sampling and Meshing of Surfaces[END_REF], yields the result.

Proof of termination

To prove that the refinement algorithm terminates, we prove that the distance between any two points inserted by the algorithm is bounded away from 0, which is sufficient since we assumed that M is compact.

Remember that there are two types of simplices that are refined by the algorithm. Let σ denote a k-simplex in Del T M (P) or a (k + 1)-dimensional simplex in C(P). A new point that is inserted in the picking region of σ is said to refine σ. We denote by e(σ) the minimal distance between such a new point and the current sample.

We assume without loss of generality that rch(M) = 1 for the rest of this section. We here give the hypotheses that will be used in this section.

H1. β ≥ 2 1-δ H2. ρ 0 ≥ 4 1-δ H3. Θ 0 < Eδ k N k+1 β k D , 1 B+1 H4. ε < min ξ 4β , 8ξ 1+31δ+32β , δ 8(C+1) H5. ε ≤ Θ0 2β(C+2)
In the hypotheses ξ is the constant defined in Lemma 8, B is defined in Lemma 9, C = 4ρ0 Θ k 0 , D is defined in Lemma 9, E will be defined in Lemma 13, N will be defined in the proof of Lemma 15, and ε, δ, β, ρ 0 and Θ 0 are parameters of the algorithm.

Observe that once δ is fixed in [0, 1), β and ρ 0 can be fixed so as to satisfy H1 and H2. Then, we can fix Θ 0 so that H3 is satisfied, and lastly we can fix ε. H5 provides a trade-off between improving the quality of the simplices (by fixing a high Θ 0 ) and minimizing the size of the sample.

Lemma 10 Let p be a point on M and q be a point on T p such that p-q ≤ 1/4. Then qpick(q, p) ≤ 2 pq 2 .

Proof.

Let A = B(p, r) ∩ M where r = 2 pq . Let f : A → T p be the orthogonal projection map of A to T p . It is proved in [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF] (Lemma 5.3) that B(p, r cos θ) ∩ T p ⊆ f (A) where sin θ = pq ≤ 1/4. pick(q, p) returns the point x closest to q in M ∩ F where F is a (dk)dimensional flat passing through q and parallel to N p . Since q ∈ B(p, r cos θ)∩T p and B(p, r cos θ) ∩ T p ⊆ f (A), x ∈ A. Therefore, from Lemma 1 (1) and the fact that px ≤ 2 pq , we have

q -x = p -x sin(T p , px) ≤ 2 p -q 2 .
The following Lemmas 11, 12 and 15 will bound the minimum interpoint distance.

Lemma 11 (Rule 1) If τ is a k-simplex of star(p) for which Rule 1 is applied, i.e. r p (τ ) > ε, then e(τ ) ≥ r p (τ )/2 > ε/2.

Proof.

Let x = pick(c p (τ ), p) be the point inserted by Rule 1 to refine τ . Since P is a 1/16-sample of M, it follows from Lemma 5 that r p (τ ) ≤ 1/4.

Using c p (τ )p = r p (τ ), r p (τ ) ≤ 1/4, and Lemma 10, we get

c p (τ ) -x ≤ 2 r 2 p ≤ r p (τ )/2. Therefore x ∈ B(c p (τ ), r p (τ )/2). For any vertex v inserted before x, we have v -x ≥ r p (τ ) -c p (τ ) -x ≥ r p (τ )/2 > ε/2 .

RR n°7660

Lemma 12 (Rule 2) Under Hypotheses H2 and H4, for a simplex σ being refined by Rule 2, i.e. σ ∈ C(P) with ρ σ > ρ 0 , we have e(σ) ≥ r σ /2 > ρ 0 L σ /2 > 2L σ .

Proof.

1. Consider first the case where σ = τ is a k-simplex of star(p) for some p. Let x = pick(c p (τ ), p) be the point inserted for refining τ . Using the fact that P is a 1/16-sample of M, and arguments similar to the ones used in the proof of Lemma 11, for any vertex v inserted before x, we have

v -x ≥ r p (τ )/2 ≥ r τ /2 > ρ 0 L τ /2 ≥ 2L τ .
The last inequality follows from the fact ρ 0 ≥ 4 (Hypothesis H2).

2. Consider now the case where σ = φ is an inconsistent configuration in C(P) witnessed by p, q and r, and let τ = φ \ {r} be a k-dimensional simplex. By definition of an inconsistent configuration, τ belongs to star(p). Since φ belongs to C(P), we have r p (τ ) ≤ ε (by the definition of C(P)) and from Corollary 1, sin ∠(aff(τ ), T p ) ≤ 4ρ 0 ε/Θ k 0 = Cε, as τ is a good simplex. Let x = pick(i φ , p) be the point inserted by Rule 2 to refine φ. Let i ′ denote the projection of i φ onto T p and i ′′ = pick(i ′ , p).

From Hypothesis H4, we have ε ≤ δ 8(1+C) which implies Cε < 1/2 (a crude bound for simplicity). Using the same arguments as in the proof of Lemma 6 and sin ∠(aff(τ ), T p ) ≤ Cε, we have

c τ -c p (τ ) , c τ -i φ ≤ r τ tan ∠(aff(τ ), T p ) ≤ Cε r τ √ 1 -C 2 ε 2 ≤ 2Cε r τ (8)
and

r φ ≤ R φ ≤ r τ + c τ -i φ ≤ (1 + 2Cε) r τ ≤ 2r τ . (9)
Using the facts that pi ′ ≤ R φ ≤ 2 r τ , r τ ≤ r p (τ ) ≤ ε, and Lemma 10, we have

i ′ -i ′′ ≤ 2 p -i ′ 2 ≤ 4εR φ . (10)
Since i ′ is the projection of i φ onto T p , hence

i φ -i ′ ≤ i φ -c p (τ ) ≤ i φ -c τ + c τ -c p (τ ) ≤ 4Cεr τ , (11)
the last inequality follows from inequality (8).

Since the line segments i φ i ′ , i ′ i ′′ are parallel to N p , hence the line segment i φ i ′′ is parallel to N p . From the definition of x = pick(i φ , p), inequalities ( 10) and [START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF] and ε ≤ δ 8(1+C) , we have

i φ -x ≤ i φ -i ′′ ≤ i φ -i ′ + i ′ -i ′′ ≤ 4εR φ + 4Cεr τ ≤ 4ε(1 + C)R φ ≤ R φ /2.
Let v be a vertex that has been inserted before x. Since B(i φ , R φ ) is empty, vi φ ≥ R φ and therefore we have

v -x ≥ R φ /2 ≥ r φ /2 > ρ 0 L φ /2 ≥ 2L φ .

RR n°7660

The last inequality again follows from the fact that ρ 0 ≥ 4.

It follows that the shortest interpoint distance is not decreased when Rule 2 is applied.

To prove a similar result for Rule 3 and 4 we use a volume argument. The next lemma provides a lower bound on the volume of the picking regions.

Lemma 13 (Volume of Π(σ, δ)) Under Hypotheses H1 and H4, if σ is a simplex to be refined by either Rules 3 or 4, we have

vol(Π(σ, δ)) ≥ E δ k r k σ
where E is a constant > 0 and depends only on k.

Proof. 1. Consider first the case where σ = τ is a k-simplex of star(p); then Π(τ, δ) = B M (c p (τ ), δ r p (τ )). Let c be the point of M closest to c p (τ ). Since τ is being refined by Rule 3, hence r p (τ ) ≤ ε. Therefore, from Lemma 1 (2), we get

c p (τ ) -c ≤ 2ε r p (τ ) . (12)
From Hypothesis H4, we have ε ≤ From the above inequality and Lemma 8, we then have vol(Π(τ, δ)) ≥ vol(B M (c, 3δε/4))

≥ 3 k 4 k 1 - 3Aδε 4 φ k δ k r k p (τ ) ≥ 3 k 4 k 1 - 3Aξ 32 φ k δ k r k τ .
The last inequality follows from the facts that δ ≤ 1, r p (τ ) ≥ r τ and ε ≤ ξ

8

(since from Hypothesis H4 we have ε ≤ ξ 4β and from Hypothesis H1 we have β > 2).

2. Consider now the case where σ = φ is an inconsistent configuration witnessed by p, q and r. Then Π(φ, δ) = B M (i φ , δR φ ). From the definition of inconsistent configurations, the k-dimensional simplex τ = φ \ {r} is in star(p). Since φ belongs to C(P), we have r p (τ ) ≤ ε and, since τ is a good simplex, we have from Corollary 1 that sin ∠(aff(τ ), T p ) ≤ 4ρ 0 ε/Θ k 0 = Cε. As in the proof of Lemma 12 we use the crude bound Cε < 1/2 which follows from Hypothesis H4; therefore inequalities (8) and ( 9) follow. Also, using Cε < 1/2 and inequality (8), we have

r p (τ ) ≤ r τ + c p (τ ) -c τ ≤ (1 + 2Cε)r τ ≤ 2r τ . ( 13
)
Let c denote the point of M closest to c p (τ ). As in the first part of the proof, we have cc p (τ ) ≤ 2ε r p (τ ). Using inequalities (8), ( 9) and ( 13), and the fact that ε ≤ δ 8 (1+C) , we get B M (c, r) ⊆ Π(φ, δ) where

r = δR φ -i φ -c ≥ δr φ -c -c p (τ ) -c p (τ ) -c τ -c τ -i φ ≥ δr φ -2εr p (τ ) -4Cε r τ ≥ δr φ -4(1 + C)εr τ ≥ δr φ 2 .
Moreover, we deduce from inequality ( 9)

δ r φ 2 ≤ δR φ 2 ≤ δr τ ≤ δε
We then deduce, using Lemma 8,

vol(Π(φ, δ) ≥ vol(B M (c, δ r φ /2)) ≥ 1 2 k (1 -Aδε) φ k δ k r k φ ≥ 1 2 k 1 - Aξ 8 φ k δ k r k φ .
The last inequality again follows from the facts that δ < 1 and ε ≤ ξ/8.

Lemma 14 Let B p be a ball of radius R centered at a point p ∈ M and let V be a maximal set of points of B p ∩ M such that the smallest interdistance between the points is not less than 2r.

If R + r ≤ ξ, |V | = 1+Aξ 1-Aξ R r + 1 k .
Proof. Denote by B x the ball centered at x ∈ V of radius r. Plainly, for any x ∈ V , B x ⊂ B + p = B(p, R+r) and for any x, y ∈ V , B x ∩B y = ∅. By Lemma 8,

vol(B + p ∩ M) ≤ (1 + A (R + r)) φ k (R + r) k and vol(B x ) ≥ (1 -A r) φ k R k .
It follows that the number of points of V is at most

1 + A(R + r) 1 -Ar R + r r k ≤ 1 + Aξ 1 -Aξ R + r r k .
Lemma 15 (Rules 3 & 4) Under Hypotheses H1 to H5, application of Rule 3 or 4 on a simplex σ does not decrease the interpoint distance to less than (1-δ)ε 4 and does not create any tiny slivers.

Proof. The proof is by induction. Specifically, we prove that the algorithm maintains the following invariants Invariant 1 When refining a simplex σ using Rules 3 or 4, no tiny slivers of dimension ≤ k + 1 with respect to σ are created in Del(P).

Invariant 2

The interpoint distance remains greater than (1-δ)ε 4 .

We first consider the case when σ = τ is a k-simplex in star(p) to be refined by application of Rule 3. The case of an inconsistent configuration to be refined by Rule 4 is similar. Note that, since Rule 1 has not been applied, r p (τ ) ≤ ε.

Invariant 1 First observe that Invariant 1 is maintained if τ is refined by inserting a good point in Π(τ, δ).

We now prove the existence of good points in Π(τ, δ). We first show that the set of points of P that can form tiny slivers with respect to τ in Del(P) upon insertion of a point x ∈ Π(τ, δ) are at distance at most δr p (τ )+ 2βr τ < (δ + 2β)ε from c p (τ ). Indeed, recall that a tiny sliver with respect to τ has a circumradius less than β r τ and that a point x ∈ Π(τ, δ) is a good point if x does not form a tiny sliver (of dimension ≤ k + 1 and with respect to τ ) with the sample points. Hence, it is enough to consider the points of P that belong to IB = B(c p (τ ), δ r p (τ ) + 2β r τ ) since all the new simplices in Del(P) upon insertion of x are incident to x. This proves the claim.

From Lemma 13, we have Π(τ, δ) = ∅. Let c ∈ Π(τ, δ), then c p (τ )c ≤ δr p (τ ) .

We deduce

IB = B(c p (τ ), δ r p (τ ) + 2β r τ ) ⊆ B (c, (2δ + 2β) ε) def = IB + . (14)
We now apply Lemma 14 to bound the number n of sample points in

IB + . Set R = (2δ + 2β) ε, r = (1-δ)ε 8
and observe that R + r = (1+15δ+16β)ε 8 ≤ ξ by Hypothesis H4. We then get

n ≤ N def = 1 + Aξ 1 -Aξ 32(δ + β) 1 -δ + 1 k
Let µ be a good simplex with vertices in IB with r µ ≤ βr τ ≤ β ε. From Hypothesis H3, H4 and H5, Θ 0 < 1 B+1 , βε ≤ ξ/2 and (C + 2)βε < Θ 0 . We can apply Lemma 9, from which we deduce

F µ ≤ D Θ 0 r k µ ≤ D Θ 0 β k r k τ .
Consider the j-simplices, j ≤ k + 1, that are included in IB ⊆ IB + . The total number of such simplices is at most N k+1 . Hence, the total volume of the forbidden regions associated to all those simplices is at most

W = N k+1 × D Θ 0 β k r k τ . ( 15
)
On the other hand, from Lemma 13 and the fact that ε ≤ δ 8(1+C) (Hypothesis H4), we know that vol(Π(τ, δ)) ≥ Eδ k r k τ . By Hypothesis H3, the volume W of all the forbidden regions is less than vol(Π(τ, δ)), the volume of the picking region of τ . This proves the existence of good points in the picking region Π(τ, δ) of τ .

Invariant 2 We will now show that Invariant 2 is also maintained.

Let τ ′ ⊆ τ denote a simplex of τ that is a sliver. We denote by p(τ ′ ) the simplex whose killing gave birth to τ ′ . Let us now prove that the interpoint distance τ p(τ ) remains at least (1-δ)ε 4 after the insertion of x from the picking region Π(τ, δ) of τ . Let x * denote the point whose insertion killed p(τ ′ ). Observe that x * is a vertex of τ ′ , and also of τ as τ ′ ⊆ τ . See Figure 5. We distinguish the following cases.

x x * p Π(τ, δ)
Case 1 p(τ ′ ) is a big simplex killed by application of Rule [START_REF] Abresch | Injectivity Radius Estimates and Sphere Theorems[END_REF]. According to Lemma 11, the lengths of the edges incident to x * in τ ′ are greater than ε/2. The distance between x and the other points is thus greater than

(1 -δ)r p (τ ) ≥ (1 -δ)∆ τ ′ 2 ≥ (1 -δ)ε 4 .
The last inequality follows from the fact that ∆ τ ′ will be greater than the lengths of the edges of τ ′ incident to x * , which in turn are > ε/2 since the radius of p(τ ′ ) is greater than ε and we insert the new point in Π(p(τ ′ ), 1/2).

Case 2 p(τ ′ ) is a simplex with a bad radius-edge ratio killed by application of Rule [START_REF] Amenta | A Simple Algorithm for Homeomorphic Surface Reconstruction[END_REF]. From Lemma 12, we have ∆ τ ′ ≥ r p(τ ′ ) /2 > ρ 0 L p(τ ′ ) /2 and the distance between x and the other points is greater than

(1 -δ)r p (τ ) ≥ (1 -δ)∆ τ ′ 2 > (1 -δ)ρ 0 L p(τ ′ ) 4 ≥ L p(τ ′ ) ≥ (1 -δ)ε 4 .
The last two inequalities follow from Hypothesis H2 and the induction hypothesis respectively.

Case 3 p(τ ′ ) has been killed by application of Rule ( 3) or ( 4). The radius r τ ′ is bigger than β r p(τ ′ ) since, by the induction hypothesis, no tiny slivers have been created until this point. If r τ ′ > βr p(τ ′ ) then the distance between x and the other points is thus greater than

(1 -δ)r p (τ ) ≥ (1 -δ)r τ ≥ (1 -δ)r τ ′ ≥ (1 -δ)β r p(τ ′ ) ≥ (1 -δ)β L p(τ ′ ) 2 ≥ L p(τ ′ ) ≥ (1 -δ)ε 4 .
The last two inequalities follow from H1 and induction hypothesis respectively.

In all cases, the invariants are maintained after refinement of τ . This completes the proof of the lemma. The case of an inconsistent simplex φ to be refined by Rule ( 4) is similar.

We sum up the results of the section in the following theorem.

Theorem 1 Under Hypotheses H1 to H5, the algorithm terminates. If, in addition,

H6. ε ≤ Θ k+1 0 1 + 4ρ0 Θ k 0 2 + 16ρ 2 0 Θ 2 0
the algorithm removes all inconsistent configurations from Del T M (P).

Proof. Termination of the algorithm is a consequence of Lemmas 11, 12 and 15. The additional Hypothesis H6 and Lemma 7 show that all inconsistent configurations have been removed since we removed all slivers from the augmented complex C(P).

Combinatorial complexity analysis

We assume that Hypotheses H1 to H6 are satisfied. Hence the algorihm terminates and M = Del T M (P) has no inconsistencies. Before we prove the results, we define the normalized volume of M as follows:

V(M) = vol(M) rch(M) k (16)
We also assume in this section that δ ≤ 1/2.

Theorem 2

The number of points inserted by the algorithm is at most

|P| = 2 O(k) V(M) ε k , ( 17
)
where the constant of proportionality in the big-O is an absolute constant.

Proof.

Let L P denote the smallest interpoint distance of the point set P. From Lemmas 11, 12 and 15 and δ ≤ 1/2, the minimum interpoint distance in P satisfies

L P ≥ (1 -δ)ε rch(M) 4 ≥ ε rch(M) 8 .
Hence, for any p, q (p = q) in P, we have B(p, r)∩B(q, r) = ∅ where r = ε rch(M)

16

. Using the fact that ε ≤ ξ 8 (from Hypotheses H1 and H4) and Lemma 8, we have

vol(B M (p, r)) ≥ (1 -Aε 16 )φ k r k ≥ (1 -Aξ 128 )φ k r k
. By a packing argument, we get

|P| ≤ 16 k vol(M) 1 -Aξ 128 φ k ε k rch(M) k = 2 O(k) V(M) ε k . RR n°7660
The following lemma, which is a direct application of Propositions 6.2 and 6.3 from [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF], will be used in the proof of Lemma 17.

Lemma 16 Let p, q ∈ M with pq ≤ rch(M)

2

, then sin ∠(T p , T q ) ≤ 2 p-q rch(M) .

We will use Lemmas 17 and 18 to calculate the time and space complexity of the algorithm in Theorem 3.

Lemma 17 Let p ∈ M. Then, |B(p, rch(M) 2 ) ∩ P| ≤ 2 O(k) ε k
where the constant in the big-O is an absolute constant.

Proof.

We will first show that B M (p, rch(M)

2

) can be covered by 2 O(k) balls (where the constant in the big-O is an absolute constant) of radius rch(M) 6 centered on M. Then we will show that |B(x, rch(M) 6

) ∩ P| is less than 2 O(k) ε k
(the constant again in the big-O is an absolute constant) for any point x on M. Combining the two results, we will get our lemma.

x ∈ S 1 to T p is dist(x, T p ) = p -x × sin ∠(px, T p ) ≤ rch(M) 8 . ( 18
)
Using the fact that the radius of the balls B x (x ∈ S 1 ) is rch(M)

6

, and the above inequality [START_REF] Edelsbrunner | Geometry and Topology for Mesh Generation[END_REF], we get that the k-dimensional balls B x has squared radius

rch(M) 2 6 2 -dist(x, T p ) 2 ≥ 1 6 2 - 1 8 2 rch(M) 2 def = r 2 2 .
We will now bound |S 1 | using a packing argument. As the balls B x , x ∈ S 1 , are disjoint and contained in B, therefore

|S 1 | ≤ r k 1 r k 2 def = N 1 = 2 O(k) ,
where the constant in the big-O is an absolute constant.

Since S 1 is a maximal set of points such that xy ≥ rch(M) 3 for all x, y ( = x) ∈ S 1 , we claim that

B M (p, rch(M) 2 ) ⊆ x∈S1 B(x, rch(M) 3 ) . (19)

RR n°7660

Otherwise if there exist x ∈ B M (p, rch(M)

2

) \ x∈S1 B(x, rch(M) 3 ) then x -x ≥ rch(M) 3
for all x ∈ S 1 . We have reached a contradiction since we have assumed that S 1 is a maximal set of point in B M (p, rch(M) 2

) such that xy ≥ rch(M) 3 for all x, y (x = y) ∈ S 1 .

We have shown that B M (p, rch(M) 2 ) (equation ( 19)) can be covered by

2 O(k) balls centered in B M (p, rch(M) 2
) with radius rch(M)

3

. Following the same method, we can show that B(x, rch(M)

3

) can be covered by 2 O(k) (the constant in the big-O is an absolute constant) balls centered in B(x, rch(M)

3

) of radius rch(M)

6

. Therefore from inequality [START_REF] Federer | Curvature Measures[END_REF] and the above bound, we get that B M (p, rch(M)

2

) can be covered by 2 O(k) balls of radius rch(M)

6 centered in B M (p, rch(M) 2 ).
2. We will now show that for all q ∈ M, |B M (q, rch(M) 6

) ∩ P| ≤ 2 O(k)
ε k . As in the proof of Lemma 2, we have from Lemmas 11, 12 and 15 and δ ≤ 1/2, B(x, ε rch(M)

16

) ∩ B(y, ε rch(M)

16

) = ∅ for all x, y ( = x) ∈ P.

Let r = rch(M) 6 + εrch(M)

16

, and f : B M (q, r) → T q be the projection map of B M (q, r) to T q .

We will bound the volume of f (B(x, εrch(M)

16

)), where x ∈ B M (q, rch(M)

6

) ∩ P, by using the same arguments used in the proof of Lemma 5.3 in [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF].

Claim 1

The projection map f satisfies the following: (i) f is injective, and (ii) the derivative df is nonsingular for all x ∈ B M (q, r).

Proof.

1. Let α be the angle made by the segment [x 1 , x 2 ] with T q , where x 1 , x 2 ∈ B M (q, r). Using Lemmas 1 and 16 and the fact that ε < 1, we have

sin α ≤ sin ∠(x 1 x 2 , T x1 ) + sin ∠(T x1 , T q ) ≤ x 1 -x 2 2 rch(M) + 2 x 1 -q rch(M) ≤ 1 6 + ε 16 + 1 3 + ε 8 < 1 (20)
This implies f is injective. Otherwise there will exist two points x 1 , x 2 ∈ B M (q, r) such that the line segment [x 1 , x 2 ] is orthogonal to T q , but this is not possible from inequality (20).

2. If df is singular at some point x ∈ B M (q, r), then the line segment

[x, f (x)] lies in T x . As f is the projection map onto T q , therefore [x, f (x)] is parallel to N q . Since the segment [x, f (x)
] is orthogonal to T q and lies on T x , we have ∠(T q , T x ) = π/2. But from Lemma 16 and ε < 1, we have

sin ∠(T x , T q ) ≤ 2 q -x rch(M) ≤ 1 3 + ε 8 < 1.
We have reached a contradiction.

We will now bound the vol(f (B x )) where B x = B M (x, εrch(M)

16

), for all x ∈ B M (q, rch(M) 3 ) ∩ P. Let θ x is the maximal angle made by any secant s = [x, y]

RR n°7660

with T q where y ∈ Bx = BM (x, εrch(M)

16

). From Lemmas 1, 16, and ε < 1, we get

sin(θ x ) ≤ max y∈Bx sin ∠(xy, T x ) + sin ∠(T x , T q ) ≤ max y∈Bx x -y 2rch(M) + 2 q -x rch(M) ≤ ε 32 + 1 3 < 0.80 (21)
Since f is nonsingular at x and therefore locally invertible, hence there exists a ball of radius r centered on x such that f -1 (B(f (x), r)∩T q ) ⊆ B x . Let r x denote the maximal radius such that for all r < r x , we have

f -1 (B(f (x), r x )∩T q ) ⊆ B x . By definition r x is such that f -1 (B(f (x), r x ) ∩ T q ) ⊆ B x .
This can happen only when there exist a point y ∈ Bx = BM (x, ε 16 ) such that either f is singular at y or else y ∈ B x . As we have shown in Claim 1 (ii) that f is nonsingular at all points in B M (q, r) ⊃ B x , hence x ∈ Bx \ B x . Which implies that xy = εrch(M) 16 and the angle made by the segment [x, y] with T q is ≤ θ x (by definition of θ x ). Hence r x ≥ εrch(M)

16 cos θ x . Therefore vol(f (B x )) ≥ vol(B(f (x), r x ) ∩ T q ) = φ k ε k rch(M) k 16 k cos k θ x . (22)
Since the balls

B x = B M (x, εrch(M)

16

) for all x ∈ B M (q, rch(M)

3

) are disjoint and f is injective, we get from inequalities ( 22) and ( 21)

vol(f (∪ x∈S B x )) = x∈S vol(f (B x )) = |S| ε k rch(M) k 2 O(k)
where S = B M (x, rch(M)

3

) ∩ P. Using the fact that f ( x∈S B x ) ⊆ B(q, r) ∩ T q , we have

|S| ≤ vol(B(q, r) ∩ T q ) ε k rch(M) k 2 O(k) = 2 O(k) ε k .

Lemma 18

The expected number of times pick() is called within good-pick(σ, δ) is 1 T , where

T = 1 -O Θ 0 δ k .
The constant in the big-O depends on k and ρ 0 .

Proof.

In the algorithm, good-pick() is called either by Rule 3 to refine a k-simplex in Del T M (P) or by Rule 4 to refine an inconsistent configuration in C(P). We will consider the two cases separately.

1. We will first consider the case when σ = τ is a k-dimensional simplex in star(p) ⊂ Del T M (P). Since τ is a k-dimensional simplex in star(p), hence it is being refined by Rule 3 and r p (τ ) ≤ ε rch(M). Let B = B M (c p (τ ), δr p (τ )), and let f : M → T p denote the orthogonal projection map of M onto T p . As in the proof of Lemma 17, we can prove that the map f restricted to the set B is injective and the derivative d f is nonsingular for all points in B, which implies that f is an open map when restricted to the set B.

Function goodpick(τ, δ) picks a random point x ∈ B p = B(c p (τ ), δr p (τ )) ∩ T p and checks whether the two following conditions are satisfied: (C1) x ′ = pick(x, p) is in B, and (C2) x ′ does not form a j-dimension sliver (2 ≤ j ≤ k+1) with other sample points contained in the ball B(c p (τ ), δr p (τ ) + 2β r τ ). If both conditions (C1) and (C2) are satisfied, then return x.

We will now bound the volume of the set

S 1 = {x ∈ B p | pick(x, p) ∈ B},
i.e. the set of points in B p that do not satisfy (C1).

Claim 2 S 1 ⊆ B p \ f (B) Proof. Let x ∈ S 1 . Since x ∈ S 1
, implies that pick(x, p) is either empty, i.e. (dk)-flat, H x , passing through x and parallel to N p does not intersect M or x ′ = pick(x, p) ∈ B. We claim that there does not exist a point in B whose image under the map f is x. Otherwise if there exist a point y ∈ B such that f (y) = x, then this would imply that the line segment [x, y] lies in H x . This would imply that H x ∩ M is not empty and

x -y 2 = y -c p (τ ) 2 -x -c p (τ ) 2 (Pythagoras theorem) < δ 2 r p (τ ) 2 -x -c p (τ ) 2 (since y ∈ B) ≤ x ′ -c p (τ ) 2 -x -c p (τ ) 2 (since x ′ ∈ B) = x -x ′ 2 (Pythagoras theorem) (23)
We have reached a contradiction since H x ∩ M =, and yx < x ′x but by definition x ′ = pick(x, p) is the point closest to x in H x ∩ M. This implies that x ∈ B p \ f (B) and the claim follows.

From the Claim 2, we have

vol(S 1 ) ≤ vol(B p \ f (B)) = vol(B p ) -vol( f (B)) = φ k δ k r k p (τ ) -vol( f (B)). (24)
We will upper bound vol(S 1 ) by lower bounding vol( f (B)).

Let p ′ be the point closest to c p (τ ) on M. From Lemma 1 (2) we have pp ′ ≤ 2ε pc p (τ ) = 2ε r p (τ ). Therefore, B ′ = B(p ′ , r) ⊆ B where r = (δ -2ε)r p (τ ).

As in the proof of Lemma 17, using the fact that f is an open map when restricted to B, we can show that

vol( f (B ′ )) ≥ φ k r k cos k θ, ( 25
)
where θ is the maximal angle made by any secant s = [p ′ , x] with T p where x ∈ B′ = B(p ′ , r). Using Lemmas 1 and 16, and ε < 1, we get

sin θ ≤ (δ -2ε)ε 2 + 2ε + 4ε 2 < 3 √ ε. (26)

RR n°7660

Therefore using inequalities ( 25) and ( 26), we get

vol( f (B ′ )) ≥ φ k r k cos k θ ≥ φ k 1 - 2ε δ k δ k r k p (τ )(1 -9ε) k 2 ≥ φ k δ k r p (τ ) k 1 - 2kε δ 1 - 9kε 2 = 1 -O ε δ φ k δ k r p (τ ) k , (27)
where the constant in the big-O depends on k. Then, from equation ( 24) and inequality ( 27), we have

vol(S 1 ) = O ε δ × φ k δ k r k p (τ ) . (28)
We will now bound the volume of the set

S 2 = x ∈ B p x ′ = pick(x, p) forms a j-dimensional sliver (2 ≤ j ≤ k + 1
) with sample points in B(c p (τ ), δr p (τ ) + βr τ ) ,

i.e. the set of points in B p whose answer to (C2) is "yes".

Let S be the set of sample points in B(c p (τ ), δr p (τ ) + βr τ ). We have shown in the proof of Lemma 15, that |S| ≤ N . Let µ be a good simplex with vertices in S with r µ ≤ βr τ . Then, from Lemma 9, we have

vol(F µ ) ≤ DΘ 0 r k µ ≤ DΘ 0 β k r k p (τ ) . ( 29
)
The number of simplices of dimension ≤ k that can be formed with vertices from S is less than N k+1 . Let the union of the forbidden regions of all the good simplices of dimension ≤ k with vertices in S be denoted by W . Then, using the fact that |S| ≤ N and inequality (29), we have

vol(W ) ≤ N k+1 × DΘ 0 β k r k p (τ ). (30)
We need the following claim to upper bound the volume of S 2 .

Claim 3 S 2 ⊆ f (W ) ∩ B p Proof. Let x ∈ S 2 . Since x ′ = pick(x, p) ∈ H x ∩M
, where H x is a (d-k)-flat passing through x and parallel to N p , then f (x ′ ) = x. As x ∈ S 2 , we also have

x ′ ∈ W . Combining the facts that S 2 ⊂ B p , f (x ′ ) ∈ W and f (x ′ ) = x, we get x ∈ f (W ) ∩ B p . Therefore S 2 ⊆ f (W ) ∩ B p .
From Claim 3 and the fact that f is a projection map, we have

vol(S 2 ) ≤ vol( f (W ) ∩ B p ) ≤ vol( f (W )) ≤ vol(W ). (31)

RR n°7660

Combining inequality ( 28), ( 30) and ( 31), we get that the probability of x ′ = pick(x, p) satisfying conditions (C1) and (C2) for any random point in x ∈ B p is greater than

T def = vol(B p \ S 1 ∪ S 2 ) vol(B p ) ≥ vol(B p ) -vol(S 1 ∪ S 2 ) vol(B p ) ≥ vol(B p ) -vol(S 1 ) -vol(S 2 ) vol(B p ) ≥ 1 -O ε δ - N k+1 β k D Θ 0 φ k δ k = 1 -O Θ 0 δ k , (32)
the constant in the big-O depends only k, ρ 0 and β, since δ ≤ 1/2, ε ≤ Θ 0 (from Hypothesis H5), D depends on k and ρ 0 (Lemma 9), and N = 2 O(k) (N depends only on β and δ, see Lemma 15). Therefore the expected number of times we have to pick random points x ∈ B p s.t x ′ = pick(x, p) satisfies both the conditions (C1) and (C2) is less than

∞ i=1 i(1 -T ) i-1 T = 1 T .
2. We can similarly show that the result holds for the case when σ = φ ∈ C(P) is an inconsistent configuration. As in the proof of Theorem 2, we can show that the smallest interpoint distance between the points of P is L P ≥ ε rch(M)/8. Using Lemma 17, we have for all x ∈ P, |B(x, rch(M)/2)

∩ (P ∪ {p})| ≤ 2 O(k) ε k . ( 33
)
The star of point x ∈ B(p, rch(M)/2) ∩ (P ∪ {p}) can be calculated by projecting all the points in B(x, rch(M)/2) ∩ (P ∪ {p}) on T x and calculating the weighted Delaunay triangulations of these projected points (Lemma 4). The time complexity for modifying the stars of all points in B(p, rch(M)/2) is this

d2 O(k) ε k + 2 O(k 2 ) ε k 2
Using the same arguments, we get that the time complexity for modifying the inconsistencies is

d 2 O(k) V(M) ε k + d 2 O(k) ε k + 2 O(k 2 ) ε k 2 ,
where the first term is for calculating the points in B(x, rch(M)/2) ∩ (P ∪ {p}), see equation [START_REF] Ruppert | A Delaunay Refinement Algorithm for Quality 2-Dimensional Mesh Generation[END_REF].

By Theorem 2, the algorithm inserts 2 O(k) vol(M)/ε k many points. From Lemma 18, we get the expected number of times pick() is called within good-pick() is 1 T , where

T = 1 -O Θ 0 δ k ,
where the constant in the big-O depends on k and ρ 0 . Hence, the total time complexity of the refinement algorithm is

(V(M) + 1) d 2 O(k) V(M) T ε 2k + 2 O(k 2 ) V(M) T ε k 2 +k .

Topological and geometric guarantees

We assume that the conditions of Theorem 1 are satisfied. Therefore M has no slivers and no inconsistencies. Let π : R d → M map each point of R d to its closest point of M. The following lemma has been proved in [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF] (except for item 5 which is a direct consequence of item 2).

Theorem 4 (Properties of M) For ε sufficiently small, we have the following :

1. M is a piecewise-linear manifold without boundary.

2. Map π restricted to M provides an isotopy between M and M.

3. ∀x ∈ M, x -π -1 (x) = O(ε 2 rch(M))
, where the constant in the big-O depends on k, ρ 0 and Θ 0 .

RR n°7660 4. ∀x ∈ M, ∠N x N τ = O(ε)
, where τ is a k-simplex of M containing the point π -1 (x).

5. The output sample P is an (ε + O(ε 2 ), Ω(ε))-sample of M, where the constant in the big-O depends on k, ρ 0 and Θ 0 , and the constant in the big-Ω depend on δ.

Conclusion

We have shown how to sample and triangulate a k-dimensional submanifold of R k up to a prescribed sampling rate ε using a variant of Delaunay refinement. The submanifold is assumed to be compact, closed and of positive reach, but not necessarily oriented. The requirement rch(M) > 0 can be somehow relaxed and Lipschitz manifolds can be triangulated in very much the same way as manifolds of positive reach, as already shown for surfaces in [START_REF] Boissonnat | Provably Good Sampling and Meshing of Lipschitz Surfaces[END_REF].

We assumed to know the reach of M (or, at least, a positive lower bound) and to be able to compute the tangent space at any point p ∈ M. If M is described by a set of equations, computing the reach of M reduces to solving a 0-dimensional system of equations [START_REF] Boissonnat | Provably Good Sampling and Meshing of Surfaces[END_REF]. Remarkably, our algorithm can be proved to tolerate some uncertainty in the estimation of the tangent spaces.

The algorithm is simple and relies only on simple computations performed in affine subspaces. In order to walk around the curse of dimensionality, we do not triangulate the ambient space and only maintain a k-dimensional data structure, the so-called tangential Delaunay complex. This leads to an algorithm that uses a restricted set of simple numerical operations and whose asymptotic complexity is O(ε -k 2 -k ) for fixed M, d and k.

We have shown that the size of the sample is O(ε -k ) and that the output mesh M is a good approximation of M from a geometric and a topological points of view. Specifically, we showed that the Hausdorff distance between M and M is O(ε 2 rch(M)) and that the maximal angle between their normal bundles is O(ε). The constant hidden in the big-O depends on the normalized volume of M (defined in Section 4.3). Up to the multiplicative constant that depends on M, those bounds match Clarkson's results [START_REF] Clarkson | Building triangulations using ε-nets[END_REF] (note that Clarkson's bound is for the Hausdorff distance only).

If one knows at each point x of M the local feature size lfs(x), it is easy to modify the algorithm so that the constants depend on M dx lfs k (x) ≤ V (M). This constant could even be improved if one combines the algorithm of this paper with the related technique developped for anisotropic mesh generation [START_REF] Boissonnat | Locally Uniform Anisotropic Meshing[END_REF]. Provided that one can estimate the second fundamental form at any point of M, such an extension would allow to construct anisotropic meshes that locally conform to the local metric of M and approximate M with a better convergence rate.

A Proof of Lemma 2

Proof of Lemma 2. 1. Without loss of generality we assume that τ = [p 0 , . . . , p j ] is embedded in R j . From the definition of fatness we have

∆ j τ Θ τ = vol(τ ) = | det(p 1 -p 0 . . . p j -p 0 )| j! ≤ ∆ j τ j! .
2. Using the bound from Lemma 2 (1) and the definition of fatness, we get

D p (τ ) = j vol(τ ) vol(τ p ) ≥ j Θ τ ∆ j τ ∆ j-1 τ (j-1)! ≥ j! Θ τ ∆ τ . We deduce, using r τ /ρ τ = L τ ≤ ∆ τ ≤ 2r τ , D p (τ ) ∆ τ = j vol(τ ) ∆ τ vol(τ p ) = j Θ τ ∆ j-1 τ Θ τp ∆ j-1 τp ≤ j Θ τ ∆ j-1 τ Θ τp L j-1 τp ≤ j 2 j-1 ρ j-1 τ × Θ τ Θ τp .
3. Let p * be the point closest to p on ∂ B(c τp , r τp ) ∩ aff(τ p ) and p ′ be the point closest to p on aff(τ p ). We denote by H the distance of c τ from aff(τ p ), Q = p ′p * , and by t the angle between qc τp and qc τ , where q is a vertex of τ p (see Figure 6). Then D p (τ ) = pp ′ and r τp = r τ cos t, which implies that

cos t = rτ p rτ ≥ Lτ p 2rτ ≥ Lτ 2rτ = 1 2ρτ
. We also have H = r τ sin t = r τp tan t.Note that the points c τ , c τp , p, p * and p ′ lie on a 2-dimensional affine space. We have to consider the following cases: (a) p ′ ∈ B(c τ , r τ ), and c τ and p lie on opposite sides of aff(c µ p * ). We have

Dp(τ ) Q ≥ rτ -H rτ p ≥ rτ -H rτ ≥ 1
sin t, see Figure 6. The distance from p to p * is less than

D p (τ ) + Q ≤ 1 + 1 1 -sin t D p (τ ) .
Q ≤ Q(2 r τp + Q) 2 r τp = p ′ -p p ′ -p ′′ 2 r τp ≤ D p (τ ) tan t RR n°7660
The distance from p to p * is less than

D p (τ ) + Q ≤ (1 + tan t) D p (τ )
The lemma follows by observing that Proof of Lemma 4. By Pythagoras theorem, ∀x ∈ T pi ∩ Vor(p j ) we have

x -p j 2 ≤ x -p k 2 ⇔ x -p ′ j 2 + p j -p ′ j 2 ≤ x -p ′ k 2 + p k -p ′ k 2 ,
where p ′ = π i (p) for p ∈ {p j , p k }. Hence Vor Tp i (P) is the power diagram (or weighted Voronoi diagram) of the weighted points (π i (p j ), w j ) ∈ T pi where w j =p jπ i (p j ) 2 . Therefore Del pi (P), which is dual to this power diagram, is a k-dimensional weighted Delaunay triangulation embedded in T pi . Plainly, Del Tp i (P) is obtained from this triangulation by lifting back the π i (p j ) onto the p j .

C Proof of Lemma 5

Proof of Lemma 5. Proof is similar to Lemma 8 in [START_REF] Boissonnat | Manifold Reconstruction using Tangential Delaunay Complexes[END_REF]. Assume for a contradiction that there exists a point x ∈ Vor(p) ∩ T p s.t. px > 4ε rch(M). Let q be a point on the line segment [px] s.t. pq = 2ε rch(M). Let q ′ be the nearest to q on M. From Lemma 1 (2), we have qq ′ ≤ 8ε 2 rch(M). Since P is an ε-sample, there exists a point t ∈ P, s.t. q ′t ≤ ε rch(M). We thus have qt ≤ qq ′ + q ′t ≤ 8ε 2 rch(M) + ε rch(M) < 2ε rch(M) ,

the last inequality follows from the fact that ε ≤ 1/8.

From Eq. 34 we get is p = t, as pq = 2ε rch(M) and tq < 2ε rch(M).

From Figure 10, we can see that ∠ptx > Π/2. This implies that

xp 2xt 2 > pt 2 > 0 , the last inequality follows from the fact that p ∈ t. This implies x ∈ Vor(p), which contradicts our initial assumption. We conclude that Vor(p) ∩ T p ⊆ B(p, 4ε rch(M)) (i). (ii) and (iii) are easy consequences of (i).

D Proof of Lemma 8

D.1 Geodesic curves and balls

The geodesic distance dist g (p, q) between points p, q ∈ M is inf |γ pq | where the infimum is taken over all the geodesic curves γ pq connecting p and q. The geodesic ball of radius r at a point p ∈ M is defined as B g (p, q) = {q : dist g (p, q) ≤ r}.

We get from Proposition 6.3 in [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF] the following lemma.

Lemma 19 Let p, q ∈ M, pq ≤ t rch(M), and t ≤ 1 2 . Then, dist g (p, q) ≤ p-q 1-t and B M (p, r) ⊆ B g (p, r/(1t)). Using the three above lemmas and Theorem 5, we get a lower bound of inj(M) in terms of rch(M) as stated in the following lemma.

Lemma 23 Let M be a k-dimensional submanifold of R d with reach rch(M).

Then inj(M) ≥ ξ 1 rch(M), where ξ 1 only depends on k.

Proof. From Lemma 20, we have, for all p ∈ M and independent vectors u and v in T p , -λ 0 ≤ K(p, u, v) ≤ λ 0 , where λ 0 = √ 2 rch(M) . Let r = t rch(M) with t ≤ 1/2 and observe that 2r √ λ 0 ≤ √ 2.

1. We can apply Lemma 21 to get The same lower bound plainly holds for inj(M) = inf p∈M inj(p).

V k -λ0 (2r) ≤ φ k (1 + 4a λ 0 r 2 ) k-1 (2r) k ≤ 2 k (1 + 2a) k-1 φ k r k def = ζ ′ r k . ( 42

D.3 Proof of Lemma 8

Once the injectivity radius of M is bounded, we can apply the following theorem from Differential Geometry that bounds the volume of geodesic balls. Refer to [START_REF] Gray | Tubes[END_REF].

Theorem 6 (The Bishop-Günther inequalities) Let M be a complete kdimensional Riemannian manifold and assume that r ≤ inj(M). Assume that there exists two constants λ and Λ such that λ ≤ K(p, u, v) ≤ Λ for all p ∈ M and independent vectors u and v in T p . Then

V k Λ (r) ≤ vol(B g (p, r)) ≤ V k λ (r) .
We can now prove Lemma 8 using Lemma 19, Lemma 20 and Theorem 6.

Proof of Lemma 8. Let r = t rch(M) and t ≤ min(

≤ φ k 1 + a λ 0 r 2 (1 -t) 2 k-1 r k (1 -t) k
(by Lemma 21)

= φ k 1 + 2a t 2 (1 -t) 2 k-1 r k (1 -t) k . ( 44
)
Observe that from inequalities (43) and (44), we deduce that there exists ξ and A that depends only on k such that for t ≤ ξ, we have

0 < 1 -A t ≤ vol(B M (p, r)) φ k r k ≤ 1 + A t .

Figure 1 :

 1 Figure 1: M is the black curve. The sample P is the set of small circles. Tangent space at p is denoted by T p . The Voronoi diagram of the sample is in grey. The edges of the Delaunay triangulation Del(P) are the line segments between small circles. In bold, star(p) = {pr, pq}. Corollary 1 (Tangent approximation) If the vertices of τ belong to M, p is a vertex of τ , and ∆ τ < rch(M), then sin ∠(T p , aff(τ )) ≤ 2ρτ ∆τ Θτ rch(M) ≤ 4ρτ rτ Θτ rch(M) .

Figure 3 :

 3 Figure 3: (a) In the figure manifold M is the black curve, the sample P is the set of small circles, tangent space at a point x ∈ P is denoted by T x , Voronoi diagram of the sample is in grey and Del T M (P) is the line segments between the sample points, in dashed lines, are the inconsistent simplices in Del T M (P). (b) In the figure the line segments denote Del T M (P) and the grey triangles denote the inconsistent configurations.

Figure 4 :

 4 Figure 4: For the proof of Lemma 9.

δ 8 ( 8 (

 88 1+C) , and ε ≤ δ C+1) implies ε < δ/8. From inequality[START_REF] Cheeger | Finite propagation speed, kernel estimates for functions of the Laplace operator, and the geometry of complete riemannian manifolds[END_REF] and the fact that ε ≤ δ/8, we get Π(τ, δ) ⊇ B M (c, (δ -2ε) r p (τ )) ⊇ B M (c, 3δε/4).

Figure 5 :

 5 Figure 5: For the proof of Lemma 15.

1 . 2 ) 6 ) 2 + 3 . 6 )

 126236 Let S 1 be the maximal set of points in B M (p, rch(M) such that xy ≥ rch(M) 3 for all x, y ( = x) ∈ S 1 . By definition for all x ∈ S 1 , the balls B x = B(x, rch(M) are disjoint. Also, these balls are contained in B = B(p, r 1 ), where r 1 = rch(M) Let us consider the k-dimensional balls B x = B x ∩ T p = B(x, rch(M) ∩ T p for all x ∈ S 1 , and B = B(p, r) ∩ T p . The balls B x are disjoint since the balls B x are disjoint. From Lemma 1 (1), the distance of

Theorem 3 1 .

 31 Under Hypotheses H1 to H5, the time complexity for updating C(P ∪ {p}) from C(P) when a new point p is inserted to the current sample P by the algorithm is O(ε -k 2 ). The expected time complexity of the algorithm is O(ε -k 2 -k ) for fixed M, d and k. Proof. Initialization. Assume without loss of generality that M is enclosed in a d-dimensional box of unit length. We partition the unit box into a grid with unit length rch(M) 32 √ d and intersect it with the manifold M to obtain the initial 1/16-sample of M. Since the complexity of the grid is 2 O(d log d) rch d (M) hence the number of points in the initial point sample, denoted by P 0 , is 2 O(d lg d) rch d (M) where the constant in the big-O is an absolute constant. The time complexity to get a subsample of the initial sample which is 1/32-sparse 1/16-sample of M is O(d |P 0 | 2 ).

2 .

 2 Refinement. When a new point p is inserted by the algorithm, Del T M (P ∪ {p}) is updated by creating the star of p and modifying the stars of all the points in B(p, rch(M)/2) ∩ P. Inconsistent configurations are only considered when all big simplices in Del T M (P) have been removed by application of Rule 1. Hence, by Lemma 6, we only have to consider inconsistent configurations with diameter at most 2R φ ≤ 4r τ ≤ 4ε rch(M) (inequality[START_REF] Boissonnat | Locally Uniform Anisotropic Meshing[END_REF]) and therefore, to update C(P ∪ {p}), it suffices to look at the stars of the points in B(p, 4ε rch(M)) ∩ (P ∪ {p}).

  Figure 6:

≥ 1 ,

 1 Figure 7: (b) p ′ ∈ B(c τ , r τ ), and c τ and p lie on the same side of aff(c τp p * ). We have Dp(τ ) Q ≥ rτ +H rτ p ≥ rτ rτ p ≥ 1, see Figure 7. The distance from p to p * is less than D p (τ ) + Q ≤ 2D p (τ ).

  Figure 8: (c) p ′ ∈ B(c τ , r τ ). From Figure 8 we can see that p should lie on the right hand side of aff(p * p). Using the facts that Q × (2 r τp + Q) = p ′p × p ′p ′′ from the Intersecting Secants Theorem, p ′ -p ′′ ≤ 2H, and H = r τp tan t, we have

Figure 9 :

 9 Figure 9: Refer to the proof of Lemma 4.

Figure 10 :

 10 Figure 10: Refer to Lemma 5. x ′ is a point on the line segment such that px ′ = 4ε rch(M), L = 4ε rch(M), ∠pax ′ = π/2 and ∠ptx ≥ ∠ptx ′ > π/2.

Lemma 20 ( 2 . 3 .dx≤ k φ k r 0 ( 1 + 2 ,

 2023012 [START_REF] Dey | Topology from Data via Geodesic complexes[END_REF]) If M is a submanifold of R d with reach rch(M), thensup p,u,v |K(p, u, v)| ≤ 2 rch 2 (M) def = λ 0 . Lemma 21 For λ ≥ 0 and r ≤ 1 √ λ , we have V k -λ (r) ≤ φ k (1 + a λ r 2 ) k-1 r k, where a is an absolute constant.Proof. 1. It is known (see[START_REF] Berger | Geometry 2. Universitext[END_REF]) thatV k λ (r) = k φ k r 0 s(λ, x) k-1 dx (0 .and φ k is the volume of the k-dimensional unit Euclidean ball. For 0 ≤ x ≤ 1, Observing that r √ λ ≤ 1 by assumption, we deduce from equation (40) and inequality (41)V k -λ (r) = k φ a λ x 2 ) k-1 x k-1 dx ≤ φ k (1 + a λ r 2 ) k-1 r k Lemma 22 ([31]) Let M be a k-dimensional submanifold of R d with reach rch(M) and let p be a point on M. Then, for r < rch(M) we have vol(B M (p, r)) ≥ φ k r k cos k θ where θ = arcsin r 2 rch(M) .

) 2 . 3 .= ξ 1

 231 By Lemma 19, we have for any pointp ∈ M, B M (p, (1t) r) ⊆ B g (p, r). It follows that vol(B g (p, r)) ≥ vol(B M (p, (1t) r)) ≥ vol(B M (p, 3r4)) ≥ φ k 3 k r k cos k θ ′ 4 k def = ζ r kwhere θ ′ = arcsin 3r/42 rch(M) < arcsin3 16 . Since r ≤ rch(rch(M) .

  It returns a good point x in Π(σ, δ). (Here σ can be a k-dimensional simplex of Del T M (P) or a (k + 1)-dimensional inconsistent configuration of C(P).)To implement pick(x, p), we use the primitive ints(M, F ) to get the set of intersection points (generically finite) and then return the intersection point closest to x.

1. pick(x, p) : The function pick(x, p) takes as input two points x ∈ R d and p ∈ M. The function returns a point closest to x from the set F ∩ M, where F is the (dk)-dimensional flat passing through x and parallel to N p . 2. good-pick(σ, δ) : This function takes as input a simplex σ and δ ∈ [0, 1).

We implement good-pick(σ, δ) as follows. If σ is a k-simplex τ in star(p), we apply the following procedure;

S1. Pick a random point y ∈ B(c p (τ ), δ r p (τ ))∩T p and calculate x = pick(y, p).

S2.

If x ∈ B(c p (τ ), δ r p (τ )) then go to S3 else go back to S1 and start over.

S3. We check if x forms a j-dimensional sliver τ 1 (2 ≤ j ≤ k + 1) with other sample points contained in the ball B(c p (τ ), δr p (τ ) + 2βr τ ). If not, x is a good point and we return x. Otherwise, we go back to S1 and start over.

  M (p, r)) ≥ φ k r k 1 -r 2 4rch(M) 2 < 2r = 2t rch(M) ≤ ξrch(M) ≤ inj(M) (fromLemma 23), we can apply Theorem 6. We can also apply Lemma 21 since r

	2. Since r ≤ r 1-t √ λ0 1-t = t √ 2 1-t ≤ 1.
	vol(B M (p, r)) ≤ vol(B g (p, r/(1 -t)) ≤ V k -λ0 (r/(1 -t))				(by Lemma 19) (by Theorem 6)
	the constant defined in Lemma 23.			1 1+ √	2 , ξ1 2 ) ≤ 1 2 , where ξ 1 is
	1. From Lemma 22, we have
	vol (B k 2
	= φ k r k 1 -	t 2 4	k 2	≥ φ k r k 1 -	k 8	t 2	(43)
	RR n°7660

RR n°7660

Acknowledgment

The authors thank Mariette Yvinec for her careful reading of earlier drafts and insightful comments.

This research has been partially supported by the Agence Nationale de la Recherche (project GAIA 07-BLAN-0328-04).

Proof. From Proposition 6.3 of [START_REF] Niyogi | Finding the Homology of Submanifolds with High Confidence from Random Samples[END_REF], we have for pq ≤ rch(M)/2 dist g (p, q) ≤ rch(M) × 1 -1 -2 pq rch(M) .

Using the fact that pq ≤ t rch(M) and inequality [START_REF] Whitney | Geometric Integration Theory[END_REF], we get

The second statement of the lemma is a direct consequence of the first one.

D.2 Injectivity radius and reach

Let γ be a geodesic curve starting at a point p ∈ M. A cut point on γ is the first point of γ where γ stops minimizing the distance to p. The cut locus CL(p) of a point p is the set of cut points of all geodesic curves of M starting at p. The injectivity radius inj(p) at point p is defined as

The injectivity radius inj(M) of M is defined as

In this section, we will bound the injectivity radius inj(M) in terms of the reach rch(M) of the manifold. We need first to recall the definition of the sectional curvature of a manifold. Given a point p ∈ M and two linearly independent vectors u, v ∈ T p , the sectional curvature is defined as

where , is the metric tensor, R() is the Riemann curvature tensor and

The following theorem is due to Cheeger et al. [START_REF] Cheeger | Finite propagation speed, kernel estimates for functions of the Laplace operator, and the geometry of complete riemannian manifolds[END_REF]Theorem 4.7]. See also [START_REF] Abresch | Injectivity Radius Estimates and Sphere Theorems[END_REF].

Theorem 5 Assume that M is a connected, complete Riemannian k-manifold such that λ ≤ K(p, u, v) ≤ Λ for all p ∈ M and independent vectors u and v in

where V k λ (̺) denotes the volume of a ball of radius ̺ in the k-dimensional space M k λ with constant sectional curvature λ.

In order to apply this theorem, we need to bound K(p, u, v), V k λ (2r) and vol(B g (p, r)). This will be done in Lemmas 20, 21 and 22 respectively.