Jean-Daniel Boissonnat
email: jean-daniel.boissonnat@inria.fr

Arijit Ghosh

Manifold reconstruction using tangential Delaunay complexes

Keywords: Tangential Delaunay complex, weighted Delaunay triangulation, manifold reconstruction, manifold learning, sampling conditions, sliver exudation

We give a provably correct algorithm to reconstruct a k-dimensional smooth manifold embedded in d-dimensional Euclidean space. The input to our algorithm is a point sample coming from an unknown manifold. Our approach is based on two main ideas : the notion of tangential Delaunay complex defined in [BF04, Flö03, Fre02], and the technique of sliver removal by weighting the sample points [CDE + 00]. Differently from previous methods, we do not construct any subdivision of the d-dimensional ambient space. As a result, the running time of our algorithm depends only linearly on the extrinsic dimension d while it depends quadratically on the size of the input sample, and exponentially on the intrinsic dimension k. To the best of our knowledge, this is the first certified algorithm for manifold reconstruction whose complexity depends linearly on the ambient dimension. We also prove that for a dense enough sample the output of our algorithm is isotopic to the manifold and a close geometric approximation of the manifold.

Introduction

Manifold reconstruction consists of computing a piecewise linear approximation of an unknown manifold M ⊂ R d from a finite sample of unorganized points P lying on M or close to M. When the manifold is a two-dimensional surface embedded in R 3 , the problem is known as the surface reconstruction problem. Surface reconstruction is a problem of major practical interest which has been extensively studied in the fields of Computational Geometry, Computer Graphics and Computer Vision. In the last decade, solid foundations have been established and the problem is now pretty well understood. Refer to Dey's book [Dey06], and the survey by Cazals and Giesen in [START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF] for recent results. The output of those methods is a triangulated surface that approximates M. This triangulated surface is usually extracted from a 3-dimensional subdivision of the ambient space (typically a grid or a triangulation). Although rather inoffensive in 3-dimensional space, such data structures depend exponentially on the dimension of the ambient space, and all attempts to extend those geometric approaches to more general manifolds have led to algorithms whose complexities depend exponentially on d [BGO09, CL08, CDR05, NSW08].

The problem in higher dimensions is also of great practical interest in data analysis and machine learning. In those fields, the general assumption is that, even if the data are represented as points in a very high dimensional space R d , they in fact live on a manifold of much smaller intrinsic dimension [SL00]. If the manifold is linear, well-known global techniques like principal component analysis (PCA) or multi-dimensional scaling (MDS) can be efficiently applied. When the manifold is highly nonlinear, several more local techniques have attracted much attention in visual perception and many other areas of science. Among the prominent algorithms are Isomap [TdSL00], LLE [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], Laplacian eigenmaps [START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF], Hessian eigenmaps [DG03], diffusion maps [LL06, NLCK05], principal manifolds [START_REF] Zhang | Principal manifolds and nonlinear dimension reduction via local tangent space alignment[END_REF]. Most of those methods reduce to computing an eigendecomposition of some connection matrix. In all cases, the output is a mapping of the original data points into R k where k is the estimated intrinsic dimension of M. Those methods come with no or very limited guarantees. For example, Isomap provides a correct embedding only if M is isometric to a convex open set of R k and LLE can only reconstruct topological balls. To be able to better approximate the sampled manifold, another route is to extend the work on surface reconstruction and to construct a piecewise linear approximation of M from the sample in such a way that, under appropriate sampling conditions, the quality of the approximation can be guaranteed. First breakthrough, along this line, was the work of Cheng, Dey and Ramos [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF] and this was followed by the work of Boissonnat, Guibas and Oudot [START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF]. In both cases, however, the complexity of the algorithms is exponential in the ambient dimension d, which highly reduces their practical relevance.

In this paper, we extend the geometric techniques developed in small dimensions and propose an algorithm that can reconstruct smooth manifolds of arbitrary topology while avoiding the computation of data structures in the ambient space. We assume that M is a smooth manifold of known dimension k and that we can compute the tangent space to M at any sample point. Under those conditions, we propose a provably correct algorithm that construct a simplicial complex of dimension k that approximates M. The complexity of the algorithm is linear in d, quadratic in the size n of the sample, and exponential in k. Our work builds on [START_REF] Boissonnat | Manifold reconstruction in arbitary dimensions using witness complexes[END_REF] and [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF] but dramatically reduces the dependence on d. To the best of our knowledge, this is the first certified algorithm for manifold reconstruction whose complexity depends only linearly on the ambient dimension. In the same spirit, Chazal and Oudot [CO08] have devised an algorithm of intrinsic complexity to solve the easier problem of computing the homology of a manifold from a sample.

Our approach is based on two main ideas : the notion of tangential Delaunay complex introduced in [START_REF] Boissonnat | A coordinate system associated with points scattered on a surface[END_REF][START_REF] Flötotto | A coordinate system associated to a point cloud issued from a manifold: definition, properties and applications[END_REF]Fre02], and the technique of sliver removal by weighting the sample points [CDE + 00]. The tangential complex is obtained by gluing local (Delaunay) triangulations around each sample point. The tangential complex is a subcomplex of the d-dimensional Delaunay triangulation of the sample points but it can be computed using mostly operations in the k-dimensional tangent spaces at the sample points. Hence the dependence on k rather than d in the complexity. However, due to the presence of so-called inconsistencies, the local triangulations may not form a triangulated manifold. Although this problem has already been reported [Fre02], no solution was known except for the case of curves (k = 1) [START_REF] Flötotto | A coordinate system associated to a point cloud issued from a manifold: definition, properties and applications[END_REF]. The idea of removing inconsistencies among local triangulations that have been computed independently has already been used for maintaining dynamic meshes [She05] and generating anisotropic meshes [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. Our approach is close in spirit to the one in [START_REF] Boissonnat | Locally uniform anisotropic meshing[END_REF]. We show that, under appropriate sample conditions, we can remove inconsistencies by weighting the sample points. We can then prove that the approximation returned by our algorithm is ambient isotopic to M, and a close geometric approximation of M.

Our algorithm can be seen as a local version of the cocone algorithm of Cheng et al. [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]. By local, we mean that we do not compute any d-dimensional data structure like a grid or a triangulation of the ambient space. Still, the tangential complex is a subcomplex of the weighted d-dimensional Delaunay triangulation of the (weighted) data points and therefore implicitly relies on a global partition of the ambient space. This is a key to our analysis and distinguishes our method from the other local algorithms that have been proposed in the surface reconstruction literature [START_REF] Cohen-Steiner | A greedy Delaunay Based Surface Reconstruction Algorithm[END_REF][START_REF] Gopi | Surface Reconstruction based on Lower Dimensional Localized Delaunay Triangulation[END_REF].

Organisation of the paper

In Section 2, we introduce the basic concepts used in this paper. We recall the notion of weighted Voronoi (or power) diagrams and Delaunay triangulations in Section 2.1 and define sampling conditions in Section 2.2. We introduce various quantities to measure the shape of simplices in Section 2.3 and, in particular, the central notion of fatness. In Section 2.4, we define the two main notions of this paper: the tangential complex and inconsistent configurations.

The algorithmic part of the paper is given in Section 3.

The main structural results are given in Section 4. Under some sampling condition, we bound the shape measure of the simplices of the tangential complex in Section 4.2 and of inconsistent configurations in Section 4.3. A crucial fact is that inconsistent configurations cannot be fat. We also bound the number of simplices and inconsistent configurations that can be incident on a point in Section 4.4. In Sections 4.5 and 4.6, we prove the correctness of the algorithm, and space and time complexity respectively. In Section 5, we prove that the simplicial complex output by the algorithm is indeed a good approximation of the sampled manifold. Finally, in Section 6, we conclude with some possible extensions.

The list of main notations has been added at the end of the paper as a reference for the readers.

Definitions and preliminaries

The standard Euclidean distance between two points x, y ∈ R d will be denoted by

x-y . Distance between a point x ∈ R d and a set X ⊂ R d , is defined as dist(x, X) = inf x∈X x -p .
We will denote the topological boundary and interior of a set X ⊆ R d by ∂X and int X respectively.

For a map f : X → Y and X 1 ⊆ X, f | X 1 denotes the restriction of map f to the subset X 1 . A j-simplex is the convex hull of j + 1 affinely independent points. If τ is a j-simples with vertices {x 0 , . . . , x j }, we will also denote τ as [x 0 , . . . , x j ]. For convenience, we often identify a simplex and the set of its vertices. Hence, if τ is a simplex, p ∈ τ means that p is a vertex of τ . If τ is a j-simplex, aff(τ ) denotes the j-dimensional affine hull of τ and N τ denotes the (d-j)-dimensional normal space of aff(τ ).

For a given simplicial complex K and a simplex σ in K, the star and link of a simplex σ in K σ ∈ K denotes the subcomplexes st(σ, K) = {τ : for some

τ 1 ∈ K, σ, τ ⊆ τ 1 } and lk(σ, K) = {τ ∈ st(σ, K) : τ ∩ σ = ∅} respectively. The open star of σ in K, st(σ, K) = st(σ, K) \ lk(σ, K).
For all x, y in R d , [x y] denotes the line segment joining the points x and y.

In this paper, M denotes a differentiable manifold of dimension k embedded in R d and P = {p 1 , . . . , p n } a finite sample of points from M. We will further assume that M has positive reach (see Section 2.2). We denote by T p and N p the k-dimensional tangent space and (d-k)-dimensional normal space at point p ∈ M respectively. xp .

If U and V are vector subspaces of R d , with dim(U ) ≤ dim(V ), the angle between them is defined by

∠(U, V ) = max u∈U min v∈V ∠(u, v),
where u and v are vectors in U and V respectively. This is the largest principal angle between U and V . In case of affine spaces, the angle between two of them is defined as the angle between the corresponding parallel vector subspaces.

The following lemma is a direct consequence of the above definition. We have included the proof for completeness.

Lemma 2.1 Let U and V be affine spaces of R d with dim(U ) ≤ dim(V ), and let U ⊥ and V ⊥ are affine spaces of

R d with dim(U ⊥ ) = d -dim(U ) and dim(V ⊥ ) = d -dim(V ). 1. If U ⊥ and V ⊥ are the orthogonal complements of U and V in R d , then ∠(U, V ) = ∠(V ⊥ , U ⊥ ). 2. If dim(U ) = dim(V ) then ∠(U, V ) = ∠(V, U ).
Proof. Without loss of generality we assume that the affine spaces U, V, U ⊥ and V ⊥ are vector subspaces of R d , i.e. they all passes through the origin.

1. Suppose ∠(U, V ) = α. Let v * ∈ V ⊥ be a unit vector. There are unit vectors u ∈ U , and u * ∈ U ⊥ such that v * = au + bu * . We will show that ∠(v * , u * ) ≤ α. First note that this angle is complementary to ∠(v * , u), i.e.,

∠(v * , u * ) = π 2 -∠(v * , u). (1)
There is a unit vector v ∈ V such that ∠(u, v) = α 0 ≤ α. Viewing angles between unit vectors as distances on the unit sphere, we exploit the triangle inequality:

∠(v * , v) ≤ ∠(v * , u) + ∠(u, v), from whence ∠(v * , u) ≥ π 2 -α 0 .
Using this expression in Eq. (1), we find

∠(v * , u * ) ≤ α 0 ≤ α, which implies, since v * was chosen arbitrarily, that ∠(V ⊥ , U ⊥ ) ≤ ∠(U, V ). Since dim V ⊥ ≤ dim U ⊥ ,
and the orthogonal complement is a symmetric relation on subspaces, the same argument yields the reverse inequality.

2. Let ∠(U, V ) = α, and let P : U → V denotes the projection map of the vector space U on V .

Case a. α = π/2. Since α = π/2 and dim(U ) = dim(V ), the projection map P is an isomorphism between vector spaces U and V . Therefore, for any unit vector v ∈ V there exist a vector u ∈ U such that P (u) = v. From the definition of angle between affine space and the linear map P , we have ∠(u, v) (= ∠(v, u)). This implies, ∠(V, U ) ≤ ∠(U, V ) = α. Using the same arguments, we can show that ∠(U,

V ) ≤ ∠(V, U ) hence ∠(U, V ) = ∠(V, U ).
Case b. α = π/2. We have ∠(V, U ) = π/2, otherwise using the same arguments as in Case 1 we can show that α = ∠(U, V ) ≤ ∠(V, U ) < π/2.

Weighted Delaunay triangulation

Weighted points

A weighted point is a pair consisting of a point p of R d , called the center of the weighted point, and a non-negative real number ω(p), called the weight of the weighted point. It might be convenient to identify a weighted point (p, ω(p)) and the hypersphere (we will simply say sphere in the sequel) centered at p of radius ω(p).

Two weighted points (or spheres) (p, ω(p)) and (q, ω(q)) are called orthogonal when pq 2 = ω(p) 2 + ω(q) 2 , further than orthogonal when pq 2 > ω(p) 2 + ω(q) 2 , and closer than orthogonal when pq 2 < ω(p) 2 + ω(q) 2 .

Given a point set P = {p 1 , . . . , p n } ⊆ R d , a weight function on P is a function ω that assigns to each point p i ∈ P a non-negative real weight ω(p i ): ω(P) = (ω(p 1 ), ..., ω(p n )). We write p ω i = (p i , ω(p i )) and P ω = {p ω 1 , . . . , p ω n }. We define the relative amplitude of ω as ω = max p∈P, q∈P\{p}

ω(p) ||p -q|| . (2)
In the paper, we will assume that ω ≤ ω 0 , for some constant ω 0 ∈ [0, 1/2) (see Hypothesis 3.2). Under this hypothesis, all the balls bounded by weighted spheres are disjoint.

Given a subset τ of d + 1 weighted points whose centers are affinely independent, there exists a unique sphere orthogonal to the weighted points of τ . The sphere is called the orthosphere of τ and its center o τ and radius Φ τ are called the orthocenter and the orthoradius of τ . If the weights of the vertices of τ are 0 (or all equal), then the orthosphere is simply the circumscribing sphere of τ whose center and radius are respectively called circumcenter and circumradius. If τ is a j-simplex, j < d, the orthosphere of τ is the smallest sphere that is orthogonal to the (weighted) vertices of τ . Its center o τ lies in aff(τ ). Note that a simplex τ may have an orthoradius which is imaginary, i.e. Φ 2 τ < 0. This situation however cannot happen if the relative amplitude of the weight function is < 1/2.

A finite set of weighted points P ω is said to be in general position if there exists no sphere orthogonal to d + 2 weighted points of P ω .

Weighted Voronoi diagram and Delaunay triangulation

Let ω be a weight function defined on P. We define the weighted Voronoi cell of p ∈ P as

Vor ω (p) = {x ∈ R d : ||p -x|| 2 -ω(p) 2 ≤ ||q -x|| 2 -ω(q) 2 , ∀q ∈ P}.
(3)

The weighted Voronoi cells and their k-dimensional faces, 0 ≤ k ≤ d, form a cell complex that decomposes R d into convex polyhedral cells. This cell complex is called the weighted Voronoi diagram or power diagram of P [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF].

Let τ be a subset of points of P and write

Vor ω (τ ) =

x ∈ τ

Vor ω (x).

We say P ω is in general position1 if Vor ω (τ ) = ∅ when |τ | > d + 1. The collection of all simplices conv(τ ) such that Vor ω (τ ) = ∅ constitutes a triangulation of conv(P) if dim aff(P) = d, and this triangulation is called the weighted Delaunay triangulation Del ω (P). The mapping that associates to the face Vor ω (τ ) of Vor ω (P) the face conv(τ ) of Del ω (P) is a duality, i.e. a bijection that reverses the inclusion relation. See [START_REF] Delaunay | Sur la sphère vide[END_REF][START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF][START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]. Alternatively, a d-simplex τ is in Del ω (P) if the orthosphere o τ of τ is further than orthogonal from all weighted points in P ω \ τ .

Observe that the definition of weighted Voronoi diagrams makes sense if, for some p ∈ P, ω(p) 2 < 0, i.e. some of the weights are imaginary. In fact, since adding a same positive quantity to all ω(p) 2 does not change the diagram, handling imaginary weights is as easy as handling real weights. In the sequel, we will only consider positive weight functions with relative amplitude < 1/2.

The weighted Delaunay triangulation of a set of weighted points can be computed efficiently in small dimensions and has found many applications, see e.g. [START_REF] Aurenhammer | Power diagrams: properties, algorithms and applications[END_REF]. In this paper, we use weighted Delaunay triangulations for two main reasons. The first one is that the restriction of a d-dimensional weighted Voronoi diagram to an affine space of dimension k is a k-dimensional weighted Voronoi diagram that can be computed without computing the d-dimensional diagram (see Lemma 2.2). The other main reason is that some flat simplices named slivers can be removed from a Delaunay triangulation by weighting the vertices (see [BGO09, CDE + 00, CDR05] and Section 3). Lemma 2.2 Let H be a k-dimensional affine space of R d . The restriction of the weighted Voronoi diagram of P = {p 0 , . . . , p m } to H is identical to the k-dimensional weighted Voronoi diagram of P = {p 0 , . . . , p m } in H, where p i denotes the orthogonal projection of p i onto H and the squared weight of

p i is ω(p i ) 2 -p i -p i 2 + λ 2
where λ = max p j ∈P p jp j is used to have all weights non-negative. In other words,

Vor ω (p i ) ∩ H = Vor ξ (p i ) where Vor ξ (p i ) = {x ∈ H : x -p i 2 -ξ(p i ) 2 ≤ x -p j 2 -ξ(p j ) 2 , ∀p j ∈ P }.

Proof.

By Pythagoras theorem, we have

∀ x ∈ H ∩ Vor ω (p i ), x -p i 2 -ω(p i ) 2 ≤ x - p j 2 -ω(p j ) 2 ⇔ x -p i 2 + p i -p i 2 -ω(p i ) 2 ≤ x -p j 2 + p j -p j 2 -ω(p j ) 2
, where p i denotes the orthogonal projection of p i ∈ P onto H. See Figure 1. Hence the restriction of Vor ω (P) to H is the weighted Voronoi diagram Vor ξ (P ) of the points P with the weight function:

ξ : P → [0, ∞), with ξ(p i ) 2 = -p i -p i 2 + ω(p i ) 2 + λ 2 where λ = max p j ∈P p j -p j . x H p i p i p j p j
Vor ω (p i p j )

Sampling conditions

Local feature size

The medial axis of M is the closure of the set of points of R d that have more than one nearest neighbor on M. The local feature size of x ∈ M, lfs(x), is the distance of x to the medial axis of M [AB99]. As is well known and can be easily proved, lfs is 1-Lipschitz, i.e., lfs(x) ≤ lfs(y) + xy .

The infimum of lfs over M is called the reach of M. In this paper, we assume that M is a compact submanifold of R d of (strictly) positive reach. Note that the class of submanifolds with positive reach includes C 2 -smooth submanifolds [Fed59].

Sampling parameters

The point sample P is said to be an (ε, δ)-sample (where 0 < δ < ε < 1) if (1) for any point x ∈ M, there exists a point p ∈ P such that ||x -p|| ≤ ε lfs(x), and (2) for any two distinct points p, q ∈ P, ||p -q|| ≥ δ lfs(p). The parameter ε is called the sampling rate, δ the sparsity, and ε/δ the sampling ratio of the sample P.

The following lemma, proved in [GW04], states basic properties of manifold samples. As before, we write nn(p) for the distance between p ∈ P and its nearest neighbor in P \ {p}.

Lemma 2.3 Given an (ε, δ)-sample P of M, we have 1. For all p ∈ P, we have

δ lfs(p) ≤ nn(p) ≤ 2ε 1 -ε lfs(p).
2. For any two points p, q ∈ M such that ||p -q|| = t lfs(p), 0 < t < 1, sin ∠(pq, T p ) ≤ t/2.

3. Let p be a point in M. Let x be a point in T p such that ||p-x|| ≤ t lfs(p) for some 0 < t ≤ 1/4. Let x be the point on M closest to x. Then ||xx || ≤ 2t 2 lfs(p).

Fat simplices

Consider a j-simplex τ , where 1 ≤ j ≤ k + 1. We denote by R τ , ∆ τ , L τ , V τ and Γ τ = ∆ τ /L τ the circumradius, the longest edge length, the shortest edge length, the j-dimensional volume, and spread i.e. the longest edge to shortest edge ratio of τ respectively. We define the fatness of a j-dimensional simplex τ as

Θ τ = 1 if j = 0 V 1/j τ /∆ τ if j > 0 (4)
The following important lemma is due to Whitney [Whi57, Chap. II].

Lemma 2.4 Let τ = [p 0 , . . . , p j ] be a j-dimensional simplex and let H be an 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)! Θ j τ L τ .
We deduce from the above lemma the following corollary. See also Lemma 1 in [START_REF] Fu | Convergence of curvature in secant approximations[END_REF] and Lemma 16 in [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF].

Corollary 2.5 (Tangent space approximation) Let τ be a j-simplex, j ≤ k, with vertices on M, and let p be vertex of τ . Assuming that ∆ τ < lfs(p), we have

sin ∠(aff(τ ), T p ) ≤ ∆ 2 τ Θ k τ L τ lfs(p) = Γ τ ∆ τ Θ k τ lfs(p)
.

Proof.

It suffices to take H = T p and to use η = ∆ 2 τ /2 lfs(p) (from Lemma 2.3 (2)) and

Γ τ = ∆ τ /L τ . Hence sin ∠(aff(τ ), T p ) ≤ 2η (j -1)! Θ j τ L τ ≤ 2η Θ j τ L τ ≤ ∆ 2 τ Θ j τ L τ lfs(p) = Γ τ ∆ τ Θ j τ lfs(p)
A sliver is a special type of flat simplex. The property of being a sliver is defined in terms of a parameter Θ 0 , to be fixed later in Section 4.

The following definition is a variant of a definition given in [START_REF] Li | Generating Well-Shaped d-dimensional Delaunay Meshes[END_REF].

Definition 2.6 (Θ 0 -fat simplices and Θ 0 -slivers) Given a positive parameter Θ 0 , a simplex τ is said to be Θ 0 -fat if the fatness of τ and of all its subsimplices is at least Θ 0 . A simplex of dimension at least 2 which is not Θ 0 -fat but whose subsimplices are all Θ 0 -fat is called a Θ 0 -sliver.

Tangential Delaunay complex and inconsistent configurations

Let P be a point sample of M and u be a point of P. We denote by Del ω u (P) the weighted Delaunay triangulation of P restricted to the tangent space T u . Equivalently, the simplices of Del ω u (P) are the simplices of Del ω (P) whose Voronoi dual faces intersect T u , i.e., τ ∈ Del ω u (P) iff Vor ω (τ ) ∩ T u = ∅. Observe that, from Lemma 2.2, Del ω u (P) is in general isomorphic to a weighted Delaunay triangulation of a k-dimensional convex hull. Since this situation can always be ensured by applying some infinitesimal perturbation on P or of the weights, in addition to the general position assumption of P ω we will also assume, in the rest of the paper, that all Del ω u (P) are isomorphic to a weighted Delaunay triangulation of a k-dimensional convex hull. We will discuss this in more details at the beginning of Section 5. Finally, write st(u) for the star of u in Del ω u (P), i.e., st(u) = st(u, Del ω u (P)) (see Figure 2). We denote by tangential Delaunay complex or tangential complex for short, the simplicial complex {τ : τ ∈ st(u), u ∈ P}. We denote it by Del ω T M (P). By our assumption above, Del ω T M (P) is a k-dimensional subcomplex of Del ω (P). Note that st(u) = st(u, Del ω u (P)) and st(u, Del ω T M (P)), the star of u in the complex Del ω T M (P), are in general different. By duality, computing st(u) is equivalent to computing the restriction to T u of the (weighted) Voronoi cell of u, which, by Lemma 2.2, reduces to computing a cell in a k-dimensional weighted Voronoi diagram embedded in T u . To compute such a cell, we need to compute the intersection of |P| -1 halfspaces of T u where |P| is the cardinality of P. Each halfspace is bounded by the bisector consisting of the points of T u that are at equal weighted distance from u ω and some other point in P ω . This can be done in optimal time [Cha93, CS89]. It follows that the tangential complex can be computed without constructing any data structure of dimension higher than k, the intrinsic dimension of M.

The tangential Delaunay complex is not in general a triangulated manifold and therefore not a good approximation of M. This is due to the presence of so-called inconsistencies. Consider a k-simplex τ of Del ω T M (P) with two vertices u and v such that τ is in st(u) but not in st(v) (refer to Figure 3). We write B u (τ ) (and B v (τ )) for the open ball centered on T u (and T v ) that is orthogonal to the (weighted) vertices of τ , and denote by m u (τ ) (and m v (τ )), or m u (m v ) for short, its center.

u v τ B v (τ ) B u (τ ) w T u Vor(τ ) aff(Vor(τ )) m u T v m v M i φ Figure 3: An inconsistent configuration in the unweighted case. Edge [u v] is in Del u (P) but not in Del v (P) since Vor(uv) intersects T u but not T v . This happens because [m u m v ] penetrates (at i φ ) the Voronoi cell of a point w = u, v, therefore creating an inconsistent configuration φ = [u, v, w].
Also note that i φ is the center of an empty sphere circumscribing simplex φ.

According to our definition, τ is inconsistent iff B u (τ ) is further than orthogonal from all weighted points in P ω \ τ while there exists a weighted point in P ω \ τ that is closer than orthogonal from B v (τ ). We deduce from the above discussion that the line segment [m u m v ] has to penetrate the interior of Vor ω (w), where w ω is a weighted point in P ω \ τ .

We now formally define an inconsistent configuration as follows.

Definition 2.7 (Inconsistent configuration) Let φ = [p 0 , . . . , p k+1 ] be a (k + 1)-simplex, and let u, v, and w be three vertices of φ. We say that φ is a Θ 0 -inconsistent (or inconsistent for short) configuration of Del ω T M (P) witnessed by the triplet • τ is a Θ 0 -fat simplex.

(u, v, w) if • The k-simplex τ = φ \ {w} is in st(u) but not in st(v). • Vor ω (w)
Note that i φ is the center of a sphere that is orthogonal to the weighted vertices of τ and also to w ω , and further than orthogonal from all the other weighted points of P ω . Equivalently, i φ is the point on [m u m v ] that belongs to Vor ω (φ). An inconsistent configuration is therefore a (k+1)-simplex of Del ω (P). However, an inconsistent configuration does not belong to Del ω T M (P) since Del ω T M (P) has no (k+1)-simplex under our general position assumption. Moreover, the lower dimensional faces of an inconsistent configuration do not necessarily belong to Del ω T M (P). Since the inconsistent configurations are k + 1-dimensional simplices, 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 write Inc ω (P) for the subcomplex of Del ω (P) consisting of all the Θ 0 -inconsistent configurations of Del ω T M (P) and their subfaces. We also define the completed complex C ω (P) as

C ω (P) = Del ω T M (P) ∪ Inc ω (P).
Refer to Figure 4. An important observation, stated as Lemma 4.9 in Section 4.3, is that, if ε is sufficiently small with respect to Θ 0 , then the fatness of φ is less than Θ 0 . Hence, if the subfaces of φ are Θ 0 -fat simplices, φ will be a Θ 0 -sliver. This observation is at the core of our reconstruction algorithm.

Manifold reconstruction

The algorithm removes all Θ 0 -slivers from C ω (P) by weighting the points of P. By the observation mentioned above, all inconsistencies in the tangential complex will then also be removed. All simplices being consistent, the resulting weighted tangential Delaunay complex M output by the algorithm will be a simplicial k-manifold that approximates M well, as will be shown in Section 5.

In this section, we describe the algorithm. Its analysis is deferred to Section 4. We will make the two following hypotheses.

Hypothesis 3.1 M is a compact smooth submanifold of R d without boundary, and P is an (ε, δ)sample of M of sampling ratio ε/δ ≤ η 0 for some positive constant η 0 .

As shown in [CWW08, GW04], we can estimate the tangent space T p at each sample point p and also the dimension k of the manifold from P and η 0 . We assume now that T p , for any point p ∈ P, k and η 0 are known.

The algorithm fixes a bound ω 0 on the relative amplitude of the weight assignment :

Hypothesis 3.2 ω ≤ ω 0 , for some constant ω 0 ∈ [0, 1/2)
Observe that, under this hypothesis, all the balls bounded by weighted spheres are disjoint.

The algorithm also fixes Θ 0 to a constant defined in Theorem 4.15, that depends on k, ω 0 and η 0 .

We define the local neighborhood of p ∈ P as

LN (p) = {q ∈ P : |B(p, p -q ) ∩ P| ≤ N }. ( 5
)
where N is a constant that depends on k, ω 0 and η 0 to be defined in Section 4.4. We will show in Lemma 4.12, that LN (p) includes all the points of P that can form an edge with p in C ω (P). In fact, the algorithm can use instead of LN (p) any subset of P that contains LN (p). This will only affect the complexity of the algorithm, not the output.

Outline of the algorithm. Initially, all the sample points in P are assigned zero weights, and the completed complex C ω (P) is built for this zero weight assignment. Then the algorithm processes each point p i ∈ P = {p 1 , . . . , p n } in turn, and assigns a new weight to p i . The new weight is chosen so that all the simplices of all dimensions in C ω (P) are Θ 0 -fat. See Algorithm 1. We now give the details of the functions used in the manifold reconstruction algorithm. The function update completed complex(Q, ω) is described as Algorithm 2. It makes use of two functions, build star(p) and build inconsistent configurations(p, τ ).

The function build star(p) calculates the weighted Voronoi cell of p, which reduces to computing the intersection of the halfspaces of T p bounded by the (weighted) bisectors between p and other points in LN (p).

The function build inconsistent configurations(u, τ ) adds to C ω (P) all the inconsistent configurations of the form φ = τ ∪ {w} where τ is an inconsistent simplex of st(u). More precisely, for each vertex v = u of τ such that τ ∈ st(v), we calculate the points w ∈ LN (p), such that (u, v, w) witnesses the inconsistent configuration φ = τ ∪ {w}. Specifically, we compute the restriction of the Voronoi diagram of the points in LN (u) to the line segment [m u m v ], where m u = T u ∩ aff(Vor ω (τ )) and m v = T v ∩ aff(Vor ω (τ )). According to the definition of an inconsistent configuration, w is one of the sites whose (restricted) Voronoi cell is the first to be intersected by the line segment [m u m v ], oriented from m u to m v . We add inconsistent configuration φ = τ ∪ {w} to the completed complex.

We now give the details of function weight(p, ω) that computes ω(p), keeping the other weights fixed (see Algorithm 3). This function extends a similar subroutine introduced in [CDE + 00] for Algorithm 1 Manifold reconstruction(P = {p 0 , . . . , p n }, η 0 ) // Initialization for i = 1 to n do calculate the local neighborhood LN (p i ) end for for i = 1 to n do ω(p i ) ← 0 end for // Build the full unweighted completed complex C ω (P) C ω (P) ← update completed complex(P, ω) // Weight assignment to remove inconsistencies

for i = 1 to n do ω(p i ) ← weight(p i , ω) update completed complex(LN (p i ), ω) end for // Output output : M ← Del ω T M (P)
Algorithm 2 Function update completed complex(Q, ω)

for each point q ∈ Q do build star(q) end for for each q ∈ Q do for each k-simplex τ in st(q) do if τ is Θ 0 -fat and ∃v ∈ τ, τ ∈ st(v) then // τ is inconsistent build inconsistent configurations(q, τ ) end if end for end for removing slivers in R 3 . We need first to define candidate simplices. A candidate simplex of p is defined as a simplex of C ω (P) that becomes incident to p when the weight of p is varied from 0 to ω 0 nn(p), keeping the weights of all the points in P \ {p} fixed. Note that a candidate simplex of p is incident to p for some weight ω(p) but does not necessarily belong to st(p).

Let τ be a candidate simplex of p that is a Θ 0 -sliver. We associate to τ a forbidden interval W (p, τ ) that consists of all squared weights ω(p) 2 for which τ appears as a simplex in C ω (P) (the weights of the other points remaining fixed).

Algorithm 3 Function weight(p, ω) S(p) ← candidate slivers(p, ω) // J(p) is the set of squared weights of p such that C ω (P) contains // no Θ 0 -sliver incident to p

J(p) ← [0, ω 2 0 nn(p) 2 ] \ τ ∈S(p) W (p, τ ) ω(p) 2 ← a squared weight from J(p) return ω(p)
The function candidate slivers(p, ω) varies the weight of p and computes all the candidate slivers of p and their corresponding weight intervals W (p, τ ). More precisely, this function follows the following steps.

1. We first detect all candidate j-simplices for all 2 ≤ j ≤ k + 1. This is done in the following way. We vary the weight of p from 0 to ω 0 nn(p), keeping the weights of the other points fixed. For each new weight assignment to p, we modify the stars and inconsistent configurations of the points in LN (p) and detect the new j-simplices incident to p that have not been detected so far. The weight of point p changes only in a finite number of instances 0 = P 0 < P 1 < • • • < P n-1 < P n = ω 0 nn(p).

2. We determine the next weight assignment of p in the following way. For each new simplex τ currently incident to p, we keep it in a priority queue ordered by the weight of p at which τ will disappear for the first time. Hence the minimum weight in the priority queue gives the next weight assignment for p. Since the number of points in LN (p) is bounded, the number of simplices incident to p is also bounded, as well as the number of times we have to change the weight of p.

3. For each candidate sliver τ of p which is detected, we compute W (p, τ ) on the fly. The candidate slivers are stored in a list with their corresponding forbidden intervals.

Analysis of the algorithm

The analysis of the algorithm relies on structural results that will be proved in Sections 4.2, 4.3 and 4.4. We will then prove that the algorithm is correct and analyze its complexity in Sections 4.5 and 4.6. In Section 5, we will show that the output M of the reconstruction algorithm is a good approximation of M.

The bounds to be given in the lemmas of this section will depend on the dimension k of M, the bound η 0 on the sampling ratio, and on a positive scalar Θ 0 that bounds the fatness and will be used to define slivers, fat simplices and inconsistent configurations.

First lemmas

Circumradius and orthoradius

The following lemma states some basic facts about weighted Voronoi diagrams when the relative amplitude of the weighting function is bounded. Similar results were proved in [CDE + 00].

Lemma 4.1 Assume that Hypothesis 3.2 is satisfied. If τ is a simplex of Del ω (P) and p and q are any two vertices of τ , then Proof. Refer to Figure 5.

1. ∀z ∈ aff(Vor ω (τ )), q -z ≤ p-z √ 1-4ω 2 0 . 2. ∀z ∈ aff(Vor ω (τ )), z -p 2 -ω 2 (p) ≥ Φ τ . 3. ∀ σ ⊆ τ , Φ σ ≤ Φ τ . o τ o τ x p q x z aff(Vor ω (τ )) τ
1. If zq ≤ zp , then the lemma is proved since 0 < 1 -4 ω 2 0 ≤ 1. Hence assume that zq > zp . Since z ∈ aff(Vor ω (τ ))

z -p 2 = z -q 2 + ω(p) 2 -ω(q) 2 ≥ z -q 2 -ω(q) 2 ≥ z -q 2 -ω 2 0 p -q 2 ≥ z -q 2 -ω 2 0 ( z -p + z -q ) 2 > z -q 2 -4ω 2 0 z -q 2 = (1 -4ω 2 0 ) z -q 2 . 2.
We know that o τ = aff(Vor ω (τ )) ∩ aff(τ ). Therefore, using Pythagoras theorem,

z -p 2 -ω(p) 2 = p -o τ 2 + o τ -z 2 -ω(p) 2 = Φ 2 τ + o τ -z 2 ≥ Φ 2 τ .
3. The result directly follows from part 2 and the fact that aff(Vor ω (τ )) ⊆ aff(Vor ω (σ)) (since σ ⊆ τ ).

Altitude and fatness

If p ∈ τ , we define τ p = τ \ {p} to be the (j -1)-face of τ opposite to p. We also write D τ (p) for the distance from p to the affine hull of τ p . D τ (p) will be called the altitude of p in τ .

From the definition of fatness, we easily derive the following lemma.

Lemma 4.2 Let τ = [p 0 , . . . , p j ] be a j-dimensional simplex and p be a vertex of τ .

1.

Θ j τ ≤ 1 j! 2. j! Θ j τ ≤ Dτ (p) ∆τ ≤ j Γ j-1 τ × Θ j τ Θ j-1 τp
Proof. 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 τ Θ j τ = V τ = | det(p 1 -p 0 . . . p j -p 0 )| j! ≤ ∆ j τ j! .
2. Using the above inequality and the definition of fatness, we get

D τ (p) = j V τ V τp ≥ j Θ j τ ∆ j τ ∆ j-1 τ (j-1)! ≥ j! Θ j τ ∆ τ . Moreover, using ∆ τp ≥ L τp ≥ L τ and Γ τ = ∆τ Lτ , D τ (p) ∆ τ = j V τ ∆ τ V τp = j Θ j τ ∆ j-1 τ Θ j-1 τp ∆ j-1 τp ≤ j Θ j τ ∆ j-1 τ Θ j-1 τp L j-1 τ = jΓ j-1 τ × Θ j τ Θ j-1 τp .

Excentricity

Let τ be a simplex and p be a vertex of τ . We define the excentricity H τ (p, ω(p)) of τ with respect to p as the signed distance from o τ to aff(τ p ). Hence, H τ (p, ω(p)) is positive if o τ and p lie on the same side of aff(τ p ) and negative if they lie on different sides of aff(τ p ).

The following lemma is a generalization of Claim 13 from [CDE + 00]. It bounds the excentricity of a simplex τ as a function of ω(p) where p is a vertex of τ . The proof is included for completeness even though it is exactly the same as the one given in [CDE + 00].

Lemma 4.3 Let τ be a simplex of Del ω (P) and let p be any vertex of τ . We have

H τ (p, ω(p)) = H τ (p, 0) - ω(p) 2 2 D τ (p) .
Proof.

Refer to Figure 6. For convenience, we write R = Φ τ for the orthoradius of τ and r = Φ τp for the orthoradius of τ p . The orthocenter o τp of τ p is the projection of o τ onto τ p . When the weight ω(p) varies while the weights of other points remain fixed, o τ moves on a (fixed) line L that passes through o τp . Now, let p and p be the projections of p onto L and aff(τ p ) respectively. We have

R 2 + ω(p) 2 = (H τ (p, ω(p)) -D τ (p)) 2 + δ 2 . We also have R 2 = H τ (p, ω(p)) 2 + r 2 and therefore H τ (p, ω(p)) 2 = (H τ (p, ω(p)) -D τ (p)) 2 + δ 2 -ω(p) 2 -r 2 . We deduce that H τ (p, ω(p)) = D τ (p) 2 + δ 2 -r 2 2D τ (p) - ω(p) 2 2D τ (p) .
The first term on the right side is H τ (p, 0) and the second is the displacement of o τ when we change the weight of p to ω 2 (p).

Properties of the tangential Delaunay complex

The following two lemmas are slight variants of results of [START_REF] Cheng | Manifold Reconstruction from Point Samples[END_REF]. The first lemma states that the restriction of a (weighted) Voronoi cell to a tangent space is small. Lemma 4.4 Assume that Hypotheses 3.1 and 3.2 are satisfied. There exists a positive constant C 1 such that for all T p ∩ Vor ω (p), ||x -p|| ≤ C 1 εlfs(p).

Proof.

Assume for a contradiction that there exists a point

x ∈ Vor ω (p) ∩ T p s.t. ||x -p|| > C 1 εlfs(p) with C 1 (1 -C 1 ε) > 2 + C 1 ε(1 + C 1 ε) (6)
Let q be a point on the line segment [px] such that ||p -q|| = C 1 εlfs(p)/2. Let q be the point closest to q on M. From Lemma 2.3, we have ||q - Hence,

q || ≤ C 2 1 ε 2 lfs(p)/2. L/2 L/2 < L/2 L/2 > π/2 p q a t x x
||p -q || ≤ ||p -q|| + ||q -q || < C 1 2 ε(1 + C 1 ε) lfs(p)
Since P is an ε-sample, there exists a point t ∈ P s.t. ||q -t|| ≤ εlfs(q ). Using the fact that lfs is 1-Lipschitz and Eq. (6), lfs(q ) ≤ lfs(p) + ||p -

q || < (1 + C 1 2 ε(1 + C 1 ε)) lfs(p) < C 1 2 (1 -C 1 ε) lfs(p), which yields ||q -t|| < C 1 2 ε(1 -C 1 ε) lfs(p). We thus have ||q -t|| ≤ ||q -q || + ||q -t|| < C 1 2 εlfs(p).
It follows (see Figure 7), that ∠ptx > π/2, which implies that

||x -p|| 2 -||x -t|| 2 -||p -t|| 2 > 0.
Hence,

||x -p|| 2 -||x -t|| 2 -ω 2 (p) + ω 2 (t) ≥ ||p -t|| 2 -ω 2 (p) ≥ ||p -t|| 2 -ω 2 0 ||p -t|| 2 > 0 (since ω 0 < 1 2 )
This implies x ∈ Vor ω (p), which contradicts our initial assumption. We conclude that Vor ω (p)∩T p ⊆ B(p, C 1 εlfs(p)) if Eq. ( 6) is satisfied, which is true for C 1 def = 3 + √ 2 ≈ 4.41 and ε ≤ ε 0 < 0.09. The following lemma states that, under Hypotheses 3.2 and 3.1, the simplices of Del ω p (P) are small, have a good radius-edge ratio and a small excentricity. Lemma 4.5 Assume that Hypotheses 3.1 and 3.2 are satisfied. There exists positive constants C 2 , C 3 and C 4 that depend on ω 0 and η 0 such that, if ε < 1 2C 2 , the following holds.

1. If pq is an edge of Del ω T M (P), then ||p -q|| < C 2 εlfs(p).

||q -x|| ≤ p -x 1 -4ω 2 0 ≤ C 1 εlfs(p) 1 -4ω 2 0 . Hence, p -q ≤ C 1 εlfs(p) where C 1 def = C 1 (1 + 1/ 1 -4ω 2 0 ).
(b) From the definition of Del ω T M (P), there exists a vertex r of τ such that [pq] ∈ st(r). From (a), rp and rq are at most C 1 εlfs(r). Using the fact that lfs is 1-Lipschitz, lfs(p

) ≥ lfs(r) -p -r ≥ (1 -C 1 ε)lfs(r) (from part 1 (a) ), which yields lfs(r) ≤ lfs(p) 1-C 1 ε . It follows that p -q ≤ p -r + r -q ≤ 2C 1 εlfs(p) 1 -C 1 ε .
The first part of the lemma is proved by taking C 2 def = 5C 1 2 and using 2C 2 ε < 1.

2. Without loss of generality, let τ 1 ∈ st(p) be a simplex incident to p and τ ⊆ τ 1 . Let z ∈ Vor ω (τ 1 ) ∩ T p , and r z = ||z -p|| 2ω 2 (p). The ball centered at z with radius r z is orthogonal to the weighted vertices of τ 1 . From Lemma 4.1 (2) and (3), we have r z ≥ Φ τ 1 and Φ τ ≤ Φ τ 1 . Hence it suffices to prove that there exists a constant C 3 such that r z ≤ C 3 L τ . Since z ∈ Vor ω (τ ) ∩ T p , we deduce from Lemma 4.4 that ||z -p|| ≤ C 1 εlfs(p). Therefore

r z = ||z -p|| 2 -ω 2 (p) ≤ ||z -p|| ≤ C 1 εlfs(p).
For any vertex q of τ , we have pq ≤ C 2 εlfs(p) (By part 1). Using 2C 2 ε < 1 and the fact that lfs is 1-Lipschitz, lfs(p) ≤ 2lfs(q). Therefore, taking for q a vertex of the shortest edge of τ , we have, using Lemma 2.3 and Hypothesis 1,

r z ≤ C 1 ε lfs(p) ≤ C 1 ε δ × δ × 2lfs(q) ≤ 2C 1 η 0 L τ .
From part (1) of the lemma we have ∆ τ 1 ≤ 2C 2 εlfs(p). Therefore

Γ τ = ∆ τ L τ ≤ ∆ τ 1 L τ ≤ 2C 2 εlfs(p) δlfs(q) ≤ 4C 2 η 0 def = C 3 .
The last inequality follows from the fact that lfs(p) ≤ 2lfs(q) and ε/δ ≤ η 0 .

3. From the definition of H τ (p, ω(p)), we have for all vertices q ∈ τ p

|H τ (p, ω(p))| = dist(o τ , aff(τ p )) = o τ -o τp ≤ o τ -q
Using the facts that Φ τ ≤ C 3 L τ (from part 2), L τ ≤ pq , ω(q) ≤ ω 0 pq , and pq ≤ C 2 εlfs(p) (from part 1), we have

o τ -q = Φ 2 τ + ω(q) 2 ≤ C 2 3 L 2 τ + ω 2 0 p -q 2 ≤ p -q C 2 3 + ω 2 0 ≤ C 2 C 2 3 + ω 2 0 × ε lfs(p) def = C 4 ε lfs(p)

Properties of inconsistent configurations

We now give lemmas on inconsistent configurations which are central to the proof of correctness of the reconstruction algorithm given later in the paper. The first lemma is the analog of Lemma 4.5 applied to inconsistent configurations. Differently from Lemma 4.5, we need to use Corollary 2.5 to control the orientation of the facets of Del ω T M (P) and require the following additional hypothesis relating the sampling rate ε and the fatness bound Θ 0 . Lemma 4.7 (Tangent approximation) Let τ be a k-simplex of Del ω T M (P) that is Θ 0 -fat. For all vertices p of τ , we have sin ∠(aff(τ ), T p ) = sin ∠(T p , aff(τ )) ≤ Aε. The constant A, defined in Hypothesis 4.6, depends on k, ω 0 , η 0 and Θ 0 .

Proof.

Since dim(τ ) = dim(T p ) (= k), we have from Lemma 2.1, ∠(aff(τ ), T p ) = ∠(T p , aff(τ )). From Lemma 4.5 (1) and the Triangle inequality, we get ∆ τ ≤ 2C 2 ε lfs(p), and, from Lemma 4.5 (2), we get Γ τ ≤ C 3 . Using Corollary 2.5, we obtain sin ∠(aff(τ

), T p ) ≤ Γ τ ∆ τ Θ k τ lfs(p) ≤ 2C 2 C 3 ε Θ k 0 = Aε
Lemma 4.8 Assume that Hypotheses 3.1, 3.2 and 4.6 are satisfied. Let φ ∈ Inc ω (P) be an inconsistent configuration witnessed by (u, v, w). There exist positive constants

C 2 > C 2 , C 3 > C 3 and C 4 > C 4 that depend on ω 0 and η 0 s.t., if ε < 1/C 2 , then 1. p -i φ ≤ C 2
2 εlfs(p) for all vertices p of φ.

2. If pq is an edge of φ then pq ≤ C 2 ε lfs(p).

If

σ ⊆ φ, then Φ σ ≤ C 3 L σ and Γ σ = ∆ σ /L σ ≤ C 3 . 4. If σ ⊆ φ and p is any vertex of σ, |H σ (p, ω(p))| of σ is at most C 4 ε lfs(p).
Proof. From the definition of inconsistent configurations, the k-dimensional simplex τ = φ \ {w} belongs to Del ω u (P). We first bound dist(i φ , aff(τ )) = o τi φ where o τ is the orthocenter of τ . Let m u (∈ Vor ω (τ ) ∩ T u ) denote, as in Section 2.4, the point of T u that is the center of the ball orthogonal to the weighted vertices of τ . By definition, m u is further than orthogonal to all other weighted points of P \τ . Observe that ||u-o τ || ≤ ||u-m u ||, since o τ belongs to aff(τ ) and therefore is the closest point to u in aff(Vor ω (τ )). Moreover, by Lemma 4.4, um u ≤ C 1 εlfs(u). Then, by Lemma 4.1, we have for all vertices p ∈ τ

p -o τ ≤ u -o τ 1 -4ω 2 0 ≤ u -m u 1 -4ω 2 0 ≤ C 1 ε lfs(u) 1 -4ω 2 0 . ( 7
)
Since τ is a Θ 0 -fat simplex (by definition of an inconsistent configuration), we get from Lemma 4.7 that sin ∠(aff(τ

), T p ) ≤ A ε (8)
for all vertices p of τ . This implies, together with 2A ε < 1 (Hypothesis 4.6),

tan 2 ∠(aff(τ ), T u ) ≤ A 2 ε 2 1 -A 2 ε 2 < 4A 2 ε 2 , (9)
Observing again that ||um u || ≤ C 1 εlfs(u) (from Lemma 4.4) and Eq. 8, we deduce

m u -o τ ≤ m u -u sin ∠(aff(τ ), T u ) ≤ A C 1 ε 2 lfs(u). (10)
We also have, ||vo τ || ≤ ||vm u || as o τ is the closest point to v in aff(Vor ω (τ )). Hence we have, using Eq.s (7) and (9),

m v -o τ ≤ v -o τ tan ∠(aff(τ ), T v ) < 2 A C 1 ε 2 1 -4ω 2 0 lfs(u) (11)
Let i φ denote, as in Section 2.4, the first point of the line segment [m u m v ] that is in Vor ω (φ). We get from Eq.s (10) and (11) that

o τ -i φ ≤ 2A C 1 ε 2 1 -4ω 2 0 lfs(u).
1. Using Lemma 4.4, and the facts that 2A ε < 1 and uo τ ≤ um u , we get

u -i φ ≤ u -o τ + o τ -i φ ≤ u -m u + o τ -i φ ≤ C 1 ε + 2A C 1 ε 2 1 -4ω 2 0 lfs(u) ≤ C 1 ε + C 1 ε 1 -4ω 2 0 lfs(u) ≤ C 2 4 εlfs(u) (12)
where C 2 is the constant introduced in Lemma 4.5. Eq. ( 12), together with Lemma 4.1 and

C 2 def = C 2 / 1 -4ω 2 0 > C 2 , yields p -i φ ≤ u -i φ 1 -4ω 2 0 ≤ C 2 4 εlfs(u) (13)
for all vertices p of φ. We deduce

p -u ≤ p -i φ + u -i φ ≤ C 2
2 εlfs(u). We now express lfs(u) in terms of lfs(p) using the fact that lfs is 1-Lipschitz and using C 2 ε < 1:

|lfs(p) -lfs(u)| ≤ p -u ≤ C 2 2 εlfs(u) ≤ 1 2 lfs(u). (14)
We deduce that lfs(u) ≤ 2lfs(p) and pi φ ≤

C 2 2 εlfs(p) 2. Using Eq. ( 13) and ( 14), from part 1 of this lemma, we have

p -q ≤ p -i φ + q -i φ ≤ C 2 2 εlfs(u) ≤ C 2 εlfs(p).
3. If σ belongs to Del ω T M (P), the result has been proved in Lemma 4.5 (3) with

C 3 = C 3 . Let r φ = i φ -u 2 -ω(u) 2 . Since i φ ∈
Vor ω (σ), the sphere centered at i φ with radius r φ is orthogonal to the weighted vertices of σ. From Lemma 4.1 (2) and (3), we have r φ ≥ Φ φ and Φ φ ≥ Φ σ respectively. Hence it suffices to show that there exists a constant C 3 such that r φ ≤ C 3 L σ . Using Eq. ( 12), we get

r φ = i φ -u -ω(u) 2 ≤ i φ -u ≤ C 2 4 εlfs(u).
Let q be a vertex of a shortest edge of σ. We have, from part 2, of this lemma, u-q ≤ C 2 εlfs(q) < lfs(q). From which we deduce that lfs(u) ≤ 2lfs(q). Therefore, using Hypothesis 3.1,

Φ σ ≤ r φ ≤ C 2 2 εlfs(q) = C 2 2 ε δ × δ × lfs(q) ≤ C 2 2 η 0 L σ .
From part 2 of the lemma, we have ∆ φ ≤ 2C 2 εlfs(q), and using the facts that L σ ≥ δlfs(q) and ε/δ ≤ η 0 , this implies

Γ σ = ∆ σ L σ ≤ ∆ φ L σ ≤ 2C 2 η 0
We get the result by taking C 3 def = max{4C 2 η 0 , 2C 2 η 0 }. 4. Using the same arguments as in the proof of Lemma 4.5 (3) we have for all vertices q of σ p

|H σ (p, ω(p))| ≤ o σ -p
Using the facts that Φ σ ≤ C 3 L σ (from part 3), L σ ≤ pq , ω(q) ≤ ω 0 pq , and pq ≤ 2C 2 εlfs(p) (from part 2), we have

o σ -q = Φ 2 σ + ω(q) 2 ≤ C 3 2 L 2 σ + ω 2 0 p -q 2 ≤ p -q C 3 2 + ω 2 0 ≤ C 2 C 3 2 + ω 2 0 × ε lfs(p) def = C 4 ε lfs(p)
The next crucial lemma bounds the fatness of inconsistent configurations.

Lemma 4.9 Assume Hypotheses 3.1, 3.2 and 4.6, and ε < 1/2C 2 . The fatness Θ φ of an inconsistent configuration φ is at most

C 2 ε (k + 1)! 1 + 2C 3 Θ k 0 1 k+1
Proof. Let φ be witnessed by (u, v, w). From the definition of inconsistent configurations, the kdimensional simplex τ = φ \ {w} belongs to st(u) and τ is a Θ 0 -fat simplex. Using Lemma 4.8 (2), we have ∆ φ ≤ 2C 2 ε lfs(u). As in the proof of Lemma 4.8, we have sin

∠(T u , aff(τ )) ≤ Γτ ∆τ Θ k τ lfs(u)
(refer to Eq. ( 8)) by Corollary 2.5 and the fact that ∆ φ < lfs(u). Also, using the fact that uw < C 2 ε lfs(u) < lfs(u) (from Lemma 4.5 (2) and ε < 1/2C 2 ) and Lemma 2.3 (2), we have

sin ∠(uw, T u ) ≤ u -w 2 lfs(u) ≤ ∆ φ 2 lfs(u)
.

We can bound the altitude

D φ (w) of w in φ D φ (w) = dist(w, aff(τ )) = sin ∠(uw, aff(τ )) × u -w ≤ (sin ∠(uw, T u ) + sin ∠(aff(τ ), T u )) × ∆ φ ≤ ∆ φ 2 lfs(u) + Γ τ ∆ τ Θ k τ lfs(u) ∆ φ ≤ ∆ 2 φ 2lfs(u) 1 + 2Γ τ Θ k τ . (15)
From the definition of fatness of a simplex and Lemma 4.2 (1), we get

V τ = Θ k τ ∆ k τ ≤ ∆ k τ k! . (16)
We deduce

Θ k+1 φ = V φ ∆ k+1 φ = D φ (w) V τ (k + 1) × 1 ∆ k+1 φ ≤ ∆ 2 φ 2 lfs(u) 1 + 2Γ τ Θ k τ × ∆ k τ (k + 1)! ∆ k+1 φ
using Eq. ( 15) and ( 16)

≤ C 2 ε (k + 1)! 1 + 2C 3 Θ k 0 .
The last inequality comes from the facts that σ is Θ 0 -fat, ∆ φ ≤ 2C 2 ε lfs(u) (from Lemma 4.8 (2)) and Γ τ ≤ C 3 (from Lemma 4.5 (3)).

A consequence of the lemma is that, if the subfaces of φ are Θ 0 -fat simplices and if the following hypothesis Hypothesis 4.10

C 2 ε (k+1)! 1 + 2C 3 Θ k 0 < Θ k+1 0
is satisfied, then φ is a Θ 0 -sliver. Hence, techniques to remove slivers can be used to remove inconsistent configurations.

In the above lemmas, we assumed that ε is small enough. Specifically in addition to Hypotheses 3.1, 3.2 and 4.6, we assumed that 2C 2 ε < 1 in Lemma 4.5, C 2 ε < 1 in Lemma 4.8 and 2C 2 ε < 1 in Lemma 4.9. We will make another hypothesis that subsumes these two previous conditions.

Hypothesis 4.11 C 2 (1 + C 2 η 0 ) ε < 1/2.
Observe that this hypothesis implies C 2 (1 + C 2 ) ε < 1/2 since η 0 > 1. Observe also that Hypotheses 4.6, 4.10 and 4.11 are satisfied for fixed ω 0 , η 0 , Θ 0 and a sufficiently small ε.

Number of local neighbors

We will use the result from this section for the analysis of the algorithm, and also for calculating its time and space complexity.

Let N def = (4C 2 η 0 + 6) k , where the constant C 2 is defined in Lemma 4.8.

Lemma 4.12 Assume Hypotheses 3.1, 3.2, 4.6 and 4.11. The set

LN (p) = {q ∈ P : |B(p, p -q ) ∩ P| ≤ N } ,
where N = 2 O(k) and the constant in big-O depends on ω 0 and η 0 , includes all the points of P that can form an edge with p in C ω (P).

Proof. Lemmas 4.5 and 4.8 show that, in order to construct st(p) and search for inconsistencies involving p, it is enough to consider the points of P that lie in ball B p = B(p, C 2 εlfs(p)). Therefore it is enough to count the number of points in B p ∩ P. Let x and y be two points of B p ∩ P. Since lfs() is a 1-Lipschitz function, we have

lfs(p) (1 -C 2 ε) ≤ lfs(x), lfs(y) ≤ lfs(p) (1 + C 2 ε). (17)
By definition of an (ε, δ)-sample of M, the two balls B x = B(x, r x ) and B y = B(y, r y ), where r x = δ lfs(x)/2 and r y = δ lfs(y)/2, are disjoint. Moreover, both balls are contained in the ball B + p = B(p, r + ), where

r + = C 2 ε lfs(p) + (1 + C 2 ε)δ lfs(p). Let B x = B x ∩ T p , B y = B y ∩ T p and B + p = B + p ∩ T p . From Lemma 2.3(2), the distance from x to T p is dist(x, T p ) = p -x × sin(px, T p ) ≤ C 2 2 ε 2 lfs(p)/2. ( 18
)
Using Eq.s (17), (18) and the fact that ε/δ ≤ η 0 , we see that B x is a k-dimensional ball of squared radius

δ 2 lfs 2 (x)/4 -dist(x, T p ) 2 ≥ δ 2 lfs 2 (p)(1 -C 2 ε) 2 /4 -C 4 2 ε 4 lfs(p)/4 ≥ δ 2 lfs 2 (p)/4 × (1 -C 2 ε) 2 -C 4 2 η 2 0 ε 2 def = (r -) 2 .
We can now use a packing argument. Since the balls B x , x in B p ∩ P, are disjoint and all contained in B + p , the number of points of B p ∩ P is at most

r + r - k = (C 2 ε + (1 + C 2 ε)δ) 2 δ 2 /4 × (1 -C 2 ε) 2 -C 4 2 η 2 0 ε 2 k/2 ≤ 4 (C 2 η 0 + (1 + C 2 ε)) 2 (1 -C 2 ε) 2 -C 4 2 η 2 0 ε 2 k/2 = 4 (C 2 η 0 + (1 + C 2 ε)) 2 (1 -C 2 ε -C 2 2 η 0 ε)(1 -C 2 ε + C 2 2 η 0 ε) k/2 ≤ 4C 2 η 0 + 6 k def
= N (using Hypothesis 4.11)

And the result follows.

Correctness of the algorithm

Definition 4.13 (Sliverity range) Let ω be a weight assignment satisfying Hypothesis 3.2. The weight of all the points in P \ {p} are fixed and the weight ω(p) of p is varying. The sliverity range Σ(p) of a point p ∈ P is the measure of the set of all squared weights ω(p) 2 for which p is a vertex of a Θ 0 -sliver in C ω (P).

Lemma 4.14 Under Hypotheses 3.1, 3.2, 4.6, 4.10 and 4.11, the sliverity range satisfies

Σ(p) < 2N k+1 C 5 Θ 0 nn(p) 2
for some constant C 5 that depends on k, ω 0 and η 0 but not on Θ 0 .

Proof. Let τ be a j-dimensional simplex of C ω (P) incident on p (with 2 ≤ j ≤ k + 1). assume that τ is a Θ 0 -sliver. If ω(p) is the weight of p, we write H τ (p, ω(p)) for the excentricity of τ with respect to p and parameterized by ω(p). From Lemma 4.8(4), we have

|H τ (p, ω(p))| ≤ C 4 εlfs(p) def = D (19)
Using Lemma 4.2 (2), we have

D p (τ ) ≤ j Γ j-1 τ ∆ τ × Θ j τ Θ j-1 τp ≤ (k + 1) C k 3 Θ 0 ∆ τ def = E (20)
The last inequality follows from the facts that j ≤ k + 1, Γ τ ≤ C 3 (from Lemmas 4.5 (2) and 4.8 (3)) and τ is a Θ 0 -sliver. Moreover, from Lemma 4.3,

H τ (p, ω(p)) = H τ (p, 0) - ω(p) 2 2D p (τ ) . ( 21
)
It then follows from Eq. ( 19), ( 20) and (21) that the set of squared weights of p for which τ belongs to C ω (P) is a subset of the following interval

[2D p (τ )H τ (p, 0) -β, 2D p (τ )H τ (p, 0) + β],
where β = 2DE. Therefore, from Eq. ( 20) and ( 21), the measure of the set of weights for which τ belongs to C ω (P) is at most

2β = 4DE = 4(k + 1)C 3 k Θ 0 ∆ τ C 4 εlfs(p).
Let q 1 and q 2 be two vertices of τ such that ∆ τ = q 1q 2 . Using Lemmas 4.5 (1) and 4.8 (2), we get ∆ τ ≤ pq 1 + pq 2 ≤ 2C 2 εlfs(p).

Using this inequality, lfs(p) ≤ nn(p)/δ (Lemma 2.3) and ε/δ ≤ η 0 (Hypothesis 3.1), the sliverity

range of τ is at most 8(k + 1)C 3 k C 2 C 4 Θ 0 η 2 0 nn(p) 2 = C 5 Θ 0 nn(p) 2 with C 5 def = 8(k + 1)C 3 k C 2 C 4 η 2 0
. By Lemma 4.12, the number of j-simplices that are incident to p is at most N j . Hence, the sliverity range of p is at most

Σ(p) ≤ k+1 j=3 N j C 5 Θ 0 nn(p) 2 < 2N k+1 C 5 Θ 0 nn(p) 2 .
The last inequality follows from the fact that k+1 j=3 N j < 2N k+1 (as N = (4C 2 η 0 + 6) k > 2).

Theorem 4.15 Take Θ 0 =

ω 2 0 2N k+1 C 5
and assume that Hypotheses 3.1, 3.2, 4.6, 4.10 and 4.11 are satisfied. Then, the simplicial complex M output by Algorithm 1 has no inconsistencies and its simplices are all Θ 0 -fat.

Proof.

The sliverity range Σ(p) of p is at most 2N k+1 C 5 Θ 0 nn(p) 2 from Lemma 4.14. Since Θ 0 = ω 2 0 2N k+1 C 5 , Σ(p) is less than the total range of possible squared weights ω 2 0 nn(p) 2 . Hence, Function weight (p, ω) will always find a weight for any point p ∈ P and any weight assignment of relative amplitude at most ω 0 for the points of P \ {p}.

Since the algorithm removes all the simplices of C ω (P) that are not Θ 0 -fat, all the simplices of M are Θ 0 -fat.

By Lemma 4.9 and Hypothesis 4.10, all the inconsistent configurations in C ω (P) are either Θ 0 -slivers or contain a subfaces which are Θ 0 -slivers. It follows that M has no inconsistency since, when the algorithm terminates, all simplices of C ω (P) are Θ 0 -fat.

Time and space complexity

Theorem 4.16 Assume that Hypotheses 3.1, 3.2, 4.6, 4.10 and 4.11 are satisfied. Then the space complexity of the algorithm is

(O(d) + 2 O(k 2 ) )|P| and the time complexity is O(d) |P| 2 + d 2 O(k 2 ) |P|,
where |P| denotes the cardinality of P.

Proof. Space Complexity: For each point p ∈ P we maintain LN (p). The total space complexity for storing LN (p) for each point p ∈ P is thus O(N |P|) by definition of LN (p).

Rather than maintaining C ω (P) as one data structure, for all p ∈ P, we maintain st(p) and cosph(p), where cosph(p) = inconsistent configuration φ φ is witnessed by (p, q, r), for some q, r ∈ P \ p Note that C ω (P) = p∈P st(p) ∪ cosph(p) .

By Lemma 2.2, each st(p), p ∈ P, has the same combinatorial complexity as a Voronoi cell in the k-dimensional flat T p . Since the sites needed to compute this Voronoi cell all belong to LN (p), there number is at most N by Lemma 4.12. From the Upper Bound Theorem of convex geometry, see e.g. [BY98], the combinatorial complexity of each star is therefore O(N k/2 ). Hence the total space complexity of the tangential Delaunay complex is O(kN k/2 )|P|.

For a given inconsistent Θ 0 -fat k-simplex in st(p), we can have from Lemmas 4.8 (2) and 4.12, at most k|LN (p)| ≤ kN different inconsistent configurations. Hence, the number of inconsistent configurations to be stored in the completed complex C ω (P) is at most O(k 2 N k/2 +1 )|P|.

The main loop of the algorithm processes each point p ∈ P in turn, and finds a weight assignment ω(p) using the function weight(, ). In the function weight, we maintain a list of candidate simplices of p with their corresponding forbidden intervals and a priority queue of candidate simplices of p. By Lemmas 4.5 (1), 4.8 (2) and 4.12, the space complexities of both these data structures are bounded by O(kN k+1 ).

With N = 2 O(k) (refer to Lemma 4.12), we conclude that the total space complexity of the algorithm is Let τ = [p 0 , . . . , p k ] be a Θ 0 -fat k-simplex in st(u). For each vertex v ( = u) of τ with τ ∈ st(v), we need to compute the inconsistent configurations of the form φ = [p 0 , . . . , p k , w] witnessed by (u, v, w) where w ∈ LN (p) \ τ . The number of such inconsistent configurations is therefore less than |LN (p)| ≤ N . The time complexity to compute all the inconsistent configurations of the form φ = [p 0 , . . . , p k , w] witnessed by the triplet (u, v, w) is O(dN ). Since the number of choices of v is at most k, hence the time complexity for building all the inconsistent configurations of the form φ = [p 0 , . . . , p k , w] witnessed by (u, v, w) with v ( = u) being a vertex of τ and w a point in

O(k 2 N k/2 +1 |P| + d|P| + kN k+1 ) = (O(d) + 2 O(k 2 ) )|P|.
LN (u) \ τ is O(dkN ) (22)
The time complexity to build all the inconsistent configurations corresponding to st(u) is O(dkN k/2 +1 ) since the number of k-simplices in the star of a point p is O(N k/2 ).

Hence, the time complexity for building the inconsistent configurations of C ω (P) is O(dN + kN k/2 +1 ) |P|. Therefore the total time complexity of the initialization phase is

O(dkN + k 3 N log N + (dk + k 3 )N k/2 +1 ) |P|
Consider now the main loop of the algorithm. The time complexity of function weight(p, ω) is O((d + k 3 )N k+1 ) since we need to sweep over at most all (k + 1)-simplices incident on p with vertices in LN (p). The number of such simplices is at most N k+1 . We easily deduce from the above discussion that the time complexity of Function update complete complex(LN (p), ω) is

O(dkN + k 3 N log N + (dk + k 3 )N k/2 +1 ) N.
Since functions weight(p, Θ 0 , ω) and update complete complex(C ω (P), p, ω) are called |P| times, we conclude that the time complexity of the main loop of the algorithm is O(dk N 2 + (k 3 + dk)N k+1 ) |P|.

Combining the time complexities for all the steps of the algorithm and using N = O2 O(k) (refer to Lemma 4.12), we get the total time complexity of the algorithm

O(d)|P| 2 + O(dk N 2 + (dk + k 3 )N k+1 ) |P| = O(d) |P| 2 + d 2 O(k 2 ) |P| Observe that, since P is an (ε, δ)-sample of M with ε/δ ≤ η 0 , |P| = O(ε k ).

Topological and geometric guarantees

In this section we give conditions under which Del ω T M (P) is isotopic to and a close geometric approximation of the manifold M. The results in this section will be used to prove the theoretical guarantees on the quality of the simplicial complex M output by the Algorithm 1 in Section 3.

General position assumption

We assume that the weighted points in P ω are in general position in R d and this will imply Del ω (P) is a simplicial complex, see Section 2.1.2. Therefore Del ω T M (P) is a simplicial complex. We know from Lemma 2.2 that Vor ω (P) ∩ T p is identical to the weighted Voronoi diagram Vor ξp (P p ) of P ξp in T p , where P p is obtained by projecting P orthogonally onto T p and ξ p is the weight assignment defined in Lemma 2.2. The following lemma is a generalization of Proposition 6.2 in [NSW08] which bounds the variation of the angle between tangent spaces between two points on the manifold M.

Basic lemmas

Lemma 5.8 (Tangent variation) Let p, q ∈ M and pq = t lfs(p) with t ≤ 1. Then, sin ∠(T q , T p ) = sin ∠(T p , T q ) < O(t).

Using Lemma 5.8, we can show that the derivative of the restriction πp of π p to a neighborhood of p is nonsingular and that πp is injective. Using the Inverse function theorem this implies that πp is a diffeomorphism. Note that M p ⊆ π p ( T p ) follows from Lemma 2.3 (3).

The next lemma follows directly from Lemmas 2.3 4.7, 5.8 and 5.9.

Lemma 5.10 Under Conditions C1 and C2 of Theorem 5.3, for any x ∈ st(p, M). max{ π p (x)x , π(x)x , π * p (x)x } = O(ε 2 ) lfs(p).

We will need the following technical lemma bounding the fatness of a perturbed simplex. A similar inequality can be found in the proof of Lemma 14c in [Whi57, Chapter 4].

Lemma 5.11 Let σ = [p 0 , . . . , p j ] be a j-dimensional simplex of R d , and let σ = [p 0 , . . . , p j ] where p ip i ≤ ρ∆ σ with ρ < 1/2, for 0 ≤ i ≤ j. Then,

Θ 2j σ ≥ Θ 2j σ (1 + 2ρ) 2j - 2jρ 1 + 2ρ
This directly implies that, if ρ ≤ Θ 2j σ j2 j+1 1 -1 2 j , then Θ σ ≥ Θ σ /2.

Proof. By Lemma 4.5 (1) and the definition of σ t , we have σ t ⊂ B(p, O(ε) lfs(p)). Therefore, for a sufficiently small ε, σ t cannot intersect the medial axis of M, which implies, from Lemma 5.2 (1), that the restriction of the map π to σ t is a C 1 -function. The rest of the proof is devoted to showing that the restriction of π to σ t is injective and that dπ is nonsingular. Let ρ def = max i∈{1, ..., k} ( q ip i ) /∆ σ . Since ∆ σ ≥ δlfs(p) and ε/δ ≤ η 0 , ρ = O(ε 2 ) lfs(p)/∆ σ = O( ). It then follows from Lemma 5.11 and the fact that Θ σ ≥ Θ 0 , that Θ σt ≥ Θ 0 /2 for ε sufficiently small. Using Γ σ ≤ C 3 (Lemma 4.5 (2)), we get for sufficiently small ε

L σt ≥ L σ -2ρ∆ σ = L σ (1 -2ρΓ σ ) ≥ C 3 δlfs(p)(1 -O(ε)) ≥ C 3 δlfs(p)/2
We then deduce from Corollary 2.5 and the bounds on Θ σt and L σt

Let x be a point of σ. We have, from Lemmas 4.5 (1) and 5.10, pπ(x) = O(ε) lfs(p). Therefore, from Lemma 5.8 and Eq. ( 24), sin ∠(T π(x) , aff(σ t )) ≤ sin ∠(T π(x) , T p ) + sin ∠(T p , aff(σ t )) = O(ε).

Eq. ( 25) implies that π restricted to σ t is injective for ε sufficiently small. Indeed otherwise, there would exist x 1 , x 2 ∈ σ t such that π(x 1 ) = π(x 2 ), and the line segment in aff(σ t ) joining the points x 1 and x 2 would be parallel to N π(x 1 ) , contradicting Eq. (25).

We deduce also from Eq. (25) and Lemma 5.2 (2) that the derivative dπ is non-singular when π is restricted to σ t . This completes the proof of the lemma. Lemma 5.16 (π * C 1 -embeds σ t ) For ε sufficiently small, π * p C 1 -embeds σ t in T p for all t ∈ [0, 1].

Proof.

From the definition of σ t , for all x ∈ σ t , we have dist(x, σ) = O(ε 2 )lfs(p). From Lemma 5.10, we have xπ(x) = O(ε 2 )lfs(p) for all x ∈ σ. Combining these facts we get xπ(x) = O(ε 2 )lfs(p) for all x ∈ σ t . This implies, from Lemma 4. Since for all x ∈ π(σ t ), sin ∠(T x , T p ) < 1, we get from Lemma 5.9,

N x ∩ T p = π -1 p (x), ∀ x ∈ π(σ t ) c c ij q i z ij y ij x ij q j z ji y ji x ji B(c, R) B(c ij , r ij )
B(q i , λ(q i )) B(q j , λ(q j )) Let λ(q i ) = max{ω(q i ), bδ lfs(q i )/4} for q i ∈ {q 0 , . . . , q k }. Observe that λ(q i ) < nn(q i )/2. For ε ≤ 1/2 and Lemma 2.3 (1), λ(q i ) < nn(q i ) 2 ≤ ε 1ε lfs(q i ) ≤ 2ε lfs(q i ) .

Let x ij = [c ij q i ]∩∂B(c ij , r ij ) and y ij = [c ij q i ]∩∂B(q i , λ(q i )). Note that x ij = ∂B(c, R)∩[c ij q i ], see Figure 9. Therefore x ijy ij = r ij + λ(q i )r 2 ij + ω(q i ) 2 = 2r ij λ(q i ) + λ(q i ) 2ω(q i ) 2 r ij + λ(q i ) + r 2 ij + ω(q i ) 2 ≥ 2r ij λ(q i ) R + λ(q i ) + R 2 + λ(q i ) 2 since λ(q i ) ≥ ω(q i ) and r ij ≤ R ≥ 2r ij λ(q i ) 2C 1 ε + 2ε + 4C 2 1 ε 2 + 4ε 2 lfs(q i ) from Eq. ( 32) and (34) ≥ b 2 δ 2 lfs(q i ) 4ε(C 1 + 1 + C 2 1 + 1) using λ(q i ) ≥ bδlfs(q i ) 4 and Eq. (33) ≥ b 2 δ lfs(q i ) 4η 0 (C 1 + 1 + C 2 1 + 1) since ε/δ ≤ η 0 from Hypothesis 3.1

Using the fact that r ij ≥ x ijy ij , we have

A 2 ij def = (2r ij -x ij -y ij ) × x ij -y ij ≥ r ij × x ij -y ij ≥ C 2
6 δ 2 lfs(q i ) 2 where C 6 is a constant that depends on ω 0 and η 0 .

For a given

  p ∈ R d and r ≥ 0, B(p, r) ( B(p, r)) denotes the d-dimensional Euclidean open (close) ball centered at p of radius r, and B M (p, r) ( BM (p, r)) denotes B(p, r) ∩ M ( B(p, r) ∩ M). For a given p ∈ P, nn(p) denotes the distance of p to its nearest neighbor in P \ {p}, i.e., nn(p) = min x∈P,x =p

Figure 1 :

 1 Figure 1: Refer to Lemma 2.2. The grey line denotes the k-dimensional plane H and the black line denotes Vor ω (p i p j ).

Figure 2 :

 2 Figure2: M is the black curve. The sample P is the set of small circles. The tangent space at u is denoted by T u . 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, st(u) = {uv, uw}.

  is one of the first weighted Voronoi cells of Vor ω (P), other than the weighted Voronoi cells of the vertices of τ , that is intersected by the line segment [m u m v ] oriented from m u to m v . Here m u = T u ∩ Vor ω (τ ) and m v = T v ∩ aff(Vor ω (τ )). Let i φ denote the point where the oriented segment [m u m v ] first intersects Vor ω (w).

Figure 4 :

 4 Figure 4: In Figure (a), M is the black curve, the sample P is the set of small circles, the tangent space at a point x ∈ P is denoted by T x and the 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). In Figure (b), the grey triangles denote the inconsistent configurations corresponding to the inconsistent simplices in Figure (a).

Figure 5 :

 5 Figure 5: For the proofs of Lemma 4.1.

Figure 6 :

 6 Figure 6: For the proof of Lemma 4.3.

Figure 7 :

 7 Figure 7: Refer to Lemma 4.4. x is a point on the line segment such that ||px || = C 1 ε lfs(p), L = C 1 ε lfs(p), ∠pax = π/2 and ∠ptx ≥ ∠ptx > π/2.

2.

  If τ is a simplex of Del ω T M (P), then Φ τ ≤ C 3 L τ and Γ τ = ∆ τ /L τ ≤ C 3 . 3. If τ is a simplex of Del ω T M (P) and p a vertex of τ , the excentricity |H τ (p, ω(p))| is at most C 4 ε lfs(p). Proof. 1. We will consider the following two cases: (a) Consider first the case where pq is an edge of Del ω p (P). Then T p ∩ Vor ω (pq) = ∅. Let x ∈ T p ∩ Vor ω (pq). From Lemma 4.4, we have ||p -x|| ≤ C 1 εlfs(p). By Lemma 4.1,

Hypothesis 4. 6

 6 2A ε < 1 where A def = 2C 2 C 3 /Θ k 0 , and C 2 and C 3 are the constants defined in Lemma 4.5.

Time complexity:

  In the initialization phase, the algorithm computes LN (p) for all p ∈ P and initializes the weights to 0. This can easily be done in time O(d)|P| 2 . Then the algorithm builds C ω (P) for the zero weight assignment. The time to compute st(p) is dominated by the time to compute the cell of p in the weighted k-dimensional Voronoi diagram of the projected points of LN (p) onto T p . Since, by definition, |LN (p)| ≤ N , the time for building the star of p is the same as the time to compute the intersection of N halfspaces in R k , O(kdN + k 3 (N log N + N k/2 )), see e.g. [Cha93, BY98]. The factor O(kd) appears in the first term because calculating the projection of a point in R d on a k-flat requires to compute k inner products. TheO(k 3) factor comes from the fact that the basic operation we need to perform is to decide whether a point lies in the ball orthogonal to a k-simplex. This operation reduces to the evaluation of the sign of the determinant of a (k + 2) × (k + 2) matrix. The N k/2 term bounds the combinatorial complexity of a cell in the Voronoi diagram of N sites in a k-flat. Therefore the time needed to build the stars of all the points p in P is O(kdN + k 3 (N log N + N k/2 )) |P|.

For a point

  p ∈ P, let π p : R d → T p denotes the orthogonal projection of R d onto T p , and π p : T p \ O → M denotes the restriction of π to T p \ O. We define in addition π * p : R d \ O -→ T p as the map that maps each point x ∈ R d \ O to the point of intersection of T p and N π(x) , the normal space of M at π(x).

  Lemma 5.9 ( π p is a C 1 -diffeomorphism ) There exists absolute constants r 1 and r 2 such that π p restricted to T p def = B(p, r 1 lfs(p)) ∩ T p is a diffeomorphism onto π p ( T p ) and π p ( T p ) ⊇ M p def = B M (p, r 2 lfs(p)).

  sin ∠(aff(σ), aff(σ t )) = O(ε 2 ) lfs(p) Θ k σt L σt = O(ε)and from Lemma 4.7 sin ∠(T p , aff(σ t )) ≤ sin ∠(T p , aff(σ)) + sin ∠(aff(σ), aff(σ t )) = O(ε).

  5 (1), π(σ t ) ⊂ B M (p, O(ε) lfs(p)) ⊂ M p From Lemma 5.8 and the fact that π(σ t ) ⊂ B M (p, O(ε)lfs(p)), we have sin ∠(T x , T p ) < 1, ∀ x ∈ π(σ t )(26)

Figure 9 :

 9 Figure 9: Refer to the proof of Lemma A.1, Case 2.

This is an extension of the general position notion introduced by Delaunay [Del34] for Delaunay triagulations. See also[START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF].

Acknowledgements

The authors thank Mariette Yvinec and Ramsay Dyer for their careful reading of earlier drafts and insightful comments. The authors also gratefully acknowledge the work of the reviewers and thank them for their detailed comments that helped to improve the manuscript.

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

† This work was done when the author was a PhD student at Geometrica, INRIA Sophia Antipolis.

Under Hypotheses 3.1 and 3.2, Vor ξp (p) is bounded (Lemma 4.4) and therefore dim aff(P p ) = k. In addition to the general position assumption of P ω , we will also assume that the weighted points P ξp p are in general position on T p , meaning that there exists no orthosphere centered on T p that is orthogonal to k + 2 weighted points of P ξp p . See Sections 2.1.1 and 2.1.2. As we have already mentioned in Section 2.4, the general position assumption can be fulfilled by applying an infinitesimal perturbation to the point sample or to the weight assignment ω.

Under the general position assumption above, Del ξp (P p ) is a triangulation of conv(P p ) [Aur87, [START_REF] Aurenhammer | An optimal algorithm for constructing the weighted voronoi diagram in the plane[END_REF][START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]. Moreover, since Vor ξp (p) is bounded (from Lemma 4.4), p is an interior vertex of Del ξp (P p ). This implies the following result.

Lemma 5.1 For all p ∈ P, st(p) is isomorphic to the star of an interior vertex of a compact triangulated k-dimensional polytope.

Medial axis and the projection map

Let O denote the medial axis of M, and let π : R d \ O -→ M denote the projection map that maps each point of R d \ O to its closest point on M. The following lemma is a standard result from Federer [Fed69].

Lemma 5.2 Let M be a smooth submanifold of R d without boundary. Then, the projection map

1. The map π is a C 1 -function.

2. For all x ∈ R d \ O, the kernel of the linear map dπ(x) : R d → T π(x) , where dπ(x) denotes the derivative of π at x, is parallel to N π(x) and has dimension dk.

Main result

The main result of this section is the following theorem.

Theorem 5.3 Asume that Hypotheses 3.1 and 3.2 are satisfied as well as the following conditions:

C1. Del ω T M (P) has no k-dimensional inconsistent simplices.

C2. All the simplices in Del ω T M (P) are Θ 0 -fat.

Then there exists ε 0 that depends only on k, ω 0 , η 0 , and Θ 0 such that for all ε ≤ ε 0 , the simplicial complex Del ω T M (P), also noted M for convenience, and the map π | M satisfy the following properties:

P1. Tangent space approximation: Let τ be a k-simplex in M. For all vertices p of τ , we have sin ∠(T p , aff(τ )) = sin ∠(aff(τ ), T p ) = O(ε). Note that this property directly follows from Lemma 4.7.

P2. PL k-manifold without boundary: M is a piecewise linear k-manifold without boundary;

P3. Homeomorphism: The map π | M provides a homeomorphism between M and M;

P5. Isotopy: There exists an isotopy F : R d × [0, 1] → R d such that the map F (•, 0) restricted to M is the identity map on M and F ( M, 1) = M.

The constants in the big-O notations depend on k, ω 0 , η 0 and Θ 0 .

Since Conditions C1 and C2 are satisfied when the conditions of Theorem 4.15 are satisfied, it follows that the simplicial complex M output by Algorithm 1 satisfies the above properties P1-P5 provided that the additional condition on ε in the above theorem is satisfied.

The rest of the section will be devoted to the proof of Theorem 5.3.

Piecewise-linear k-manifold

The following lemma is a direct consequence of condition C1.

Lemma 5.4 Under condition C1 of Theorem 5.3 and for a sufficiently small ε, st(p) = st(p, M) for all p ∈ P.

We can now prove that M is a PL k-manifold without boundary (Property P2).

Lemma 5.5 (PL k-manifold without boundary) Assume condition C1 of Theorem 5.3. If ε is sufficiently small, then M is a PL k-manifold without boundary.

Proof.

Since Del ξp (P p ) is a triangulation of a k-dimensional polytope, Del ξp (P p ) is a PL kmanifold [Zee66, Lem. 8 of Chap. 3] and, since p is an interior vertex of Del ξp (P p ), lk(p, Del ξp (P p )) is a PL (k -1)-sphere [Zee66, Proof of Lem. 9 of Chap. 3].

From Lemma 5.4, we have st(p) = st(p, M). By Lemma 5.1, st(p, M) is isomorphic to the star of an interior vertex of a k-dimensional triangulated convex domain. Hence, lk(p, M) is a PL (k -1)-sphere. This directly implies that M is a PL k-manifold with no boundaries.

Property P2 follows since the output M of the Algorithm 1 has no inconsistency and therefore satisfies Condition C1.

Homeomorphism

Following lemma establishes Property P3 of Theorem 5.3, and the rest of the section is devoted to its proof.

Lemma 5.6 (Homeomorphism) Assume conditions C1 and C2 of Theorem 5.3. For ε sufficiently small, the map π restricted to M gives a homeomorphism between M and M.

In the broad outline, the proof of Lemma 5.6 is similar to the proof of homeomorphism given in [ACDL02] but the technical details are quite different in most places. The difficulties that arose in getting results analogous to the ones in [ACDL02] were handled using ideas from [Whi57, Chap. II].

In the rest of this paper, we assume that ω 0 , η 0 , Θ 0 are fixed constants, and that ε is small enough. We use the asymptotic notations O() and Ω() hiding the constants.

Remark 5.7 For the details of the missing proofs of Lemmas 5.8, 5.9, 5.10, and 5.11 refer to [Gho12, Chap. 4].

We will also need the following standard lemma from convex geometry which bounds the distance between an interior point and a point on the boundary of a simplex. The proof can be found in [Whi57, Lem. 14b of Chap. 4].

Lemma 5.12 Let σ = [p 0 , . . . , p j ] be a j-simplex and let p be a point of sigma. Writing p = j i=0 µ i p i with j i=0 µ i = 1 and µ i ≥ 0, we have dist(p, ∂σ) ≥ j!Θ j σ ∆ σ × min{µ 0 , . . . , µ j }

We will use the standard notion of oriented simplices, positively oriented simplices and oriented PL manifolds (see [START_REF] Zeeman | Seminar on Combinatorial Topology[END_REF][START_REF] Rourke | Introduction to Piecewise-Linear Topology[END_REF] or [Whi57, Appendix II]).

Recall that the restriction of π p to st(p) provides a simplicial isomorphism between st(p) and st(p, Del ξp (P p )) (Lemma 2.2), and that st(p) = st(p, M) (by Lemma 5.4). We orient T p and orient positively the k-dimensional simplices of st(p, Del ξp (P p )). We then orient the k-dimensional simplices of st(p, M) isomorphically, using the orientation of the corresponding k-dimensional simplices of st(p, Del ξp (P p )).

Since lk(p, M) is a PL (j -1)-sphere (see the proof of Lemma 5.5), we immediately get that st(p, M) is an oriented PL k-ball. Definition 5.13 (Simplexwise positive map) Let σ be a positively oriented i-simplex of R i , and let f : σ → R i be a C 1 -function. The map f is called simplexwise positive if det(J(f )) > 0 for all x ∈ σ, where J(f ) and det(J(f )) denote the Jacobian and the determinant of the Jacobian of the map f respectively.

The following result is a special case of a standard result from PL-topology (see, e.g., Appendix II of [Whi57, Lemma 15a]).

Lemma 5.14 Let K be a j-dimensional PL-manifold whose j-simplices are oriented, and let the continuous map f : K → R j be a simplexwise positive map. Then, for any connected open subset W of R j \ f (∂K), any two points of W not in f (K -1 ) are covered the same number of times. If this number is 1, then f , restricted to the open subset W = f -1 (W ) of K, is injective.

Properties of π *

p Let σ = [p, p 1 , . . . , p k ] be a k-simplex of st(p, M). Now, for all points q ∈ σ, let π p,t (q) = (1t)q + t π p (q), and σ t = π p,t (σ). Since π p is affine on each simplex, π p,t is also affine. Therefore σ t is a simplex with vertices p, p 1t , . . . , p kt with p jt = π p,t (p j ), j ∈ {1, . . . , k}.

It directly follows from Lemma 5.10 that ∀j ∈ {1, . . . , k}, p j -

Let τ be an i-simplex, and let f : τ → M be a C 1 -function. The simplex τ is said to be C 1 -embedded by f in M if f is an injective mapping and for all x ∈ τ , the rank of the linear map df (x) : Therefore, we deduce π * p σt = π -1 p • π. Refer to Fig. 8. Then, by Lemma 5.9, π -1 p is a diffeomorphism from π p ( T p ) ⊇ M p to T p , and, by Lemma 5.15,

The proof of the following lemma is similar to Lemma 23a in Chapter 4 from [Whi57] and it will be used for proving that π *

p restricted to open star st(p, M) is injective in Lemma 5.18. We give a proof for completeness.

The proof of the following lemma is similar to Lemma 23a in Chapter 4 from [Whi57]. We give a proof for completeness.

Lemma 5.17 Assume conditions C1 and C2 of Theorem 5.3. For ε sufficiently small, π * p is a simplexwise positive mapping of st(p, M) into T p .

Proof.

We use the same notation σ t as in the beginning of this section. We first observe that, since the simplex σ 1 is in T p , π * p is the identity on σ 1 . Therefore, det(J(π * p )) > 0 on σ 1 . Moreover, since from Lemma 5.16, det(J(π * p )) = 0 in σ t for all t ∈ [0, 1], we also have det(J(π * p )) > 0 in σ 0 . Which proves the lemma.

We are now ready to prove that π * p restricted to the open star st(p, M) is injective.

Lemma 5.18 (Injectivity of π * p restricted to st(p, M)) Assume conditions C1 and C2 of Theorem 5.3 and let ε be sufficiently small. For each point p ∈ P, the map π * p is injective on the open star st(p, M).

Proof. Note that, since M) is a PL-manifold without boundary (Lemma 5.5), lk(p, M) is a PL (k -1)-sphere, st(p, M) is a PL k-ball, and ∂ st(p, M) = lk(p, M).

For convenience, rewrite f = π * p and S = st(p, M) and orient S as indicated in Section 5.2.1. S is thus an oriented PL k-ball, and, by Lemma 5.17, f = π * p is a simplexwise positive mapping from S to T p . Let f (S k-1 ) be the image by f of the (k -1)-skeleton of S (i.e. the set of faces of S of dimension at most k -1) and let W be any connected open subset of T p \ f (∂S). From Lemma 5.14, any two points of W \ f (S k-1 ) are covered the same number of times. If this number is 1, then f , restricted to the open subset f -1 (W ) of st(p), is injective onto W .

To prove the lemma, it is therefore sufficient to prove that there exists a point g in S \ S k-1 whose image f (g) is not covered by any other point x of S \ S k-1 , i.e. f (g) = f (x) for all x = g and x ∈ S \ S k-1 . Let σ be a k-simplex [q 0 , . . . , q k ] of S and let

Using Lemma 5.12 and the facts that Θ σ ≥ Θ 0 (since the simplices of M are Θ 0 -fat), ∆ σ ≥ δ lfs(p), and ε/δ ≤ η 0 (Hypothesis 3.1), we have

From the above equation and using the fact that π p (x)x = O(ε 2 ) lfs(p) for all x ∈ σ (Lemma 5.10), we have for ε sufficiently small

Since π p embeds S into T p (from Lemma 2.2) with π p (S) = st(p, Del ξp (P p )) and since Del ξp (P p ) is a triangulation of a k-dimensional convex hull (see the discussion on the general position assumption at the beginning of Section 5), Eq. ( 27) implies that, for all x ∈ S \ int σ

From Eq. (28), we have for all

which is positive for a sufficiently small ε. Lemma 5.16 together with Eq. (29) shows that π * p (x) = π * p (g) for all x ∈ S \ {g}. Hence, g is not covered by any other point of S, and the lemma follows.

Properties of π

Lemma 5.19 (Injectivity of the restriction of π to st(p, M)) Assume conditions C1 and C2 of Theorem 5.3. Let ε be sufficiently small. For all p in P, the map π restricted to st(p, M) is injective.

Proof. To reach a contradiction, assume that there exist

Which contradicts the fact that π * p is injective when restricted to st(p, M) (Lemma 5.18). A standard result in topology [Mas67, Lem. 3.4 of Chap. 5] asserts that π covers all the points of a connected component of its image M the same number of times. Let us show that this number is exactly one. It cannot be zero by definition. Neither can it be more than one since, otherwise, some sample point p would be covered more than once, which would contradict Lemma A.1 from Appendix A. We have therefore proved that the restriction of π to M is injective. Since M is compact and M is a Hausdorff space, Lemmas 5.20 and 5.21, and a standard result from topology [Bre94, Thm. 7.8 of Chap. 1] implies that π : M → M is a homeomorphism.

This completes the proof of Lemma 5.6.

Pointwise approximation

Following lemma is a direct consequence is Lemmas 5.10 and 5.6, and the fact that lfs is 1-Lipschitz function. For the details of the proof, refer to [Gho12, Chap. 4].

Lemma 5.22 (Pointwise approximation) Under conditions C1 and C2 of Theorem 5.3, and

Isotopy

We need the following lemma to show the isotopy property P5 of Theorem 5.3. The proof of the lemma follows directly from [GW04, Lem. 5] and one can check the details in [Gho12, Chap. 4].

Lemma 5.23 Assume conditions C1 and C2 of Theorem 5.3. If ε is sufficiently small, for all points x, y

Consider the following map:

We will denote by F t = F (•, t) and M t = F t ( M) for all t ∈ [0, 1]. It is easy to see that F (x, t) (and F t ) is continuous, M 0 = M, and since π restricted to M gives a homeomorphism between M and M, M 1 = M. Using the facts that M is compact, M t is a Hausdorff space (as M t is a subspace with the subspace topology of R d ), and F t is continuous, we get [Bre94, Thm. 7.8 of Chap. 1], the following result:

Lemma 5.24 Assume conditions C1 and C2 of Theorem 5.3. If ε is sufficiently small, then for all t ∈ [0, 1], F t : M → M t is a homeomorphism.

This implies:

Lemma 5.25 (Isotopy) Assume conditions C1 and C2 of Theorem 5.3. If ε is sufficiently small then

is an isotopy.

Conclusion

We have given the first algorithm that is able to reconstruct a smooth closed manifold in a time that depends only linearly on the dimension of the ambient space. We believe that our algorithm is of practical interest when the dimension of the manifold is small, even if it is embedded in a space of high dimension. This situation is quite common in practical applications in machine learning. Unlike most surface reconstruction algorithms in R 3 , our algorithm does not need to orient normals (a critical issue in practical applications) and, in fact, works for non orientable manifolds.

The algorithm is simple. The basic ingredients we need are data structures for constructing weighted Delaunay triangulations in R k . We have assumed that the dimension of M is known. If not, we can use algorithms given in [CWW08, GW04] to estimate the dimension of M and the tangent space at each sample point. Moreover, our algorithm is easy to parallelize. One interesting feature of our approach is that it is robust and still works if we only have approximate tangent spaces at the sample points. We will report on experimental results in a forthcoming paper.

We have assumed that we know an upper bound on the sampling ratio η 0 of the input sample P. Ideas from [BGO09, FR02] may be useful to convert a sample to a subsample with a bounded sampling ratio.

We foresee other applications of the tangential complex and of our construction each time computations in the tangent space of a manifold are required, e.g. for dimensionality reduction and approximating the Laplace Beltrami operator [START_REF] Belkin | Discrete laplace operator on meshed surfaces[END_REF]. It easily follows from [BNN10] that our reconstruction algorithm can also be used in Bregman spaces where the Euclidean distance is replaced by any Bregman divergence, e.g. Kullback-Leibler divergence. This is of particular interest when considering statistical manifolds like, for example, spaces of images [START_REF] Carlsson | On the local behavior of spaces of natural images[END_REF].

A Appendix to Section 5.2

This section is devoted to the proof of the following lemma.

Lemma A.1 Let ε be sufficiently small. For all p in P and restricting the map π to M, we have π -1 (p) = {p}.

Proof. To reach a contradiction, we assume that there exists a k-simplex τ = [q 0 , . . . , q k ] in M such that there exists x ∈ τ with x = p and π(x) = p.

We will have to consider the following two cases.

Case 1. p is a vertex of τ . This implies that the unit vector u ∈ aff(τ ) along the line joining the points p and x lies in N p . But from Lemma 4.7 and ε sufficiently small, we have sin ∠(aff(τ ), T p ) ≤ Aε < 1.

Case 2. p is not a vertex of τ . The outline of the proof for case 2 is the following: we will divide the k-simplex τ into a union of k + 2 convex cells and show that for each convex cell the distance of any point of the cell to p is Ω(ε lfs(q)), where q is a vertex of τ . For ε sufficiently small, we will reach a contradiction from Lemma 5.10. Since M has no inconsistent configuration, τ is in st(q 0 ). Let c = Vor ω (τ ) ∩ T q 0 . From Lemma 4.4, we have

Using the facts that for all vertices q i , q j of τ , q iq j ≤ C 2 ε lfs(q i ) (from Lemma 4.5 (1)), lfs is 1-Lipschitz and ε sufficiently small, we have lfs(q j ) ≥ lfs(q i )q iq j ≥ (1 -C 2 ε) lfs(q i ) ≥ lfs(q i ) 2 (31) Therefore using Eq. ( 30) and (31), we have for all vertices q i of τ

Consider the edge q i q j of τ and let c ij be the projection of c onto the line segment [q i q j ]. Observe that the ball of radius r ij = c ijq i 2ω(q i ) 2 ≤ R centered at c ij is orthogonal to the balls B(q i , ω(q i )) and B(q j , ω(q j )). Using the fact that P is a (ε, δ)-sample of M and Lemma 4.1 (2), we have

1 -4ω 2 0 q iq j 2 ≥ bδ max{lfs(q i ), lfs(q j )} ( 33)

Let z ij denote the point closest to y ij on ∂B(c, R).

From the Intersecting Chords Theorem (see Figure 10) for circles, we have :

By the definition of z ij , the solution to the above quadratic equation in z ijy ij is the smaller root:

The last inequality follows from the facts that R ≤ 2C 1 ε lfs(q i ) (from Eq. ( 32)) and ε/δ ≤ η 0 . Using the fact that ε is sufficiently small and Eq. (31), we have for all vertices q of τ z

Let conv(S) denote the convex hull of the points y ij :

Using convexity, we have dist(conv(S), ∂B(c, R)) = min i, j ( =i) ∈ {0, ..., k}

for all vertices q ∈ {q 0 , . . . , q k }.

Recall that we have assumed in the beginning of the proof that x ∈ τ and π(x) = p. Eq. (35) implies that x ∈ conv(S). Indeed, if x ∈ conv(S), then from Eq. ( 35) and the fact that the ball B(c, R) is empty, we will have

xp ≥ C 7 ε lfs(q) for all vertices q of τ . But from Lemma 5.10, we have

for all vertices q of τ . We have reached a contradiction for ε sufficiently small. Hence x ∈ conv(S).

Let S i = {y ij | j ∈ {0, . . . , k} \ {i}} ∪ {q i }. The convex hulls of S i , conv(S i ), satisfies the following properties : conv(S i ) ⊂ B(q i , λ(q i ))

If x ∈ conv(S i ), then from Lemma 5.10, we have for ε sufficiently small

The last inequality holds for ε < 1 2Cη 0 . We have reached a contradiction. Since x ∈ conv(S i ) and x ∈ conv(S), it follows that x ∈ τ . We have therefore reached a contradiction, i.e. there does not exit x ∈ τ such that π(x) = p.

Main notations

General notations (Section 2)

line segment connecting points x and y xy standard Euclidean distance between the points x and y ∠ (u, v) angle between the vectors u and v ∠(U, V ) angle between the affine spaces U and V Submanifold (Section 2)