Jean-Daniel Boissonnat

Ramsay Dyer

Arijit Ghosh

Boissonnat

Dyer

Ghosh

Constructing Intrinsic Delaunay Triangulations of Submanifolds

Keywords: Delaunay triangulation, submanifolds, Riemannian geometry triangulation de Delaunay, sous-variétés, géométrie riemannienne Constructing intrinsic DTs

We describe an algorithm to construct an intrinsic Delaunay triangulation of a smooth closed submanifold of Euclidean space. Using results established in a companion paper on the stability of Delaunay triangulations on δ-generic point sets, we establish sampling criteria which ensure that the intrinsic Delaunay complex coincides with the restricted Delaunay complex and also with the recently introduced tangential Delaunay complex. The algorithm generates a point set that meets the required criteria while the tangential complex is being constructed. In this way the computation of geodesic distances is avoided, the runtime is only linearly dependent on the ambient dimension, and the Delaunay complexes are guaranteed to be triangulations of the manifold.

Introduction

This paper addresses the problem of constructing an intrinsic Delaunay triangulation of a smooth closed submanifold M ⊂ R N . We present an algorithm which generates a point set P ⊂ M and a simplicial complex on P that is homeomorphic to M and has a connectivity determined by the Delaunay triangulation of P with respect to the intrinsic metric of M.

For a submanifold of Euclidean space, the restricted Delaunay complex [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF], which is defined by the ambient metric restricted to the submanifold, was employed by Cheng et al. [START_REF] Cheng | Manifold reconstruction from point samples[END_REF] as the basis for a triangulation. However, it was found that sampling density alone was insufficient to ensure a triangulation, and manipulations of the complex were employed.

In an earlier work, Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF] announced sampling density conditions which would ensure that the Delaunay complex defined by the intrinsic metric of the manifold was a triangulation. In fact, as shown in Section 2.4.3 and Appendix A, the stated result is incorrect: sampling density alone is insufficient to guarantee an intrinsic Delaunay triangulation (see Theorem A.3). Topological defects can arise when the vertices lie too close to a degenerate or "quasi-cospherical" configuration.

Our interest in the intrinsic Delaunay complex stems from its close relationship with other Delaunay-like structures that have been proposed in the context of non-homogeneous metrics. For example, anisotropic Voronoi diagrams [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF] and anisotropic Delaunay triangulations emerge as natural structures when we want to mesh a domain of R m while respecting a given metric tensor field.

This paper builds over preliminary results on anisotropic Delaunay meshes [START_REF] Boissonnat | Anisotropic Delaunay mesh generation[END_REF] and manifold reconstruction using the tangential Delaunay complex [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. The central idea in both cases is to define Euclidean Delaunay triangulations locally and to glue these local triangulations together by removing inconsistencies between them. We view the inconsistencies as arising from instability in the Delaunay triangulations, and exploit the results of a companion paper [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] to define sampling conditions under which these inconsistencies cannot arise.

The algorithm is based on the tangential Delaunay complex [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF], and is an adaptation of a Delaunay refinement algorithm designed to avoid poorly shaped "sliver" simplices [START_REF] Li | Generating well-shaped d-dimensional Delaunay meshes[END_REF][START_REF] Boissonnat | Triangulating smooth submanifolds with light scaffolding[END_REF]. The tangential Delaunay complex is defined with respect to local Delaunay triangulations restricted to the tangent spaces at sample points. We demonstrate that the algorithm produces sampling conditions such that the tangential Delaunay complex coincides with the restricted Delaunay complex and the intrinsic Delaunay complex. The refinement algorithm avoids the problem of slivers without the need to resort to a point weighting strategy [CDE + 00, [START_REF] Cheng | Manifold reconstruction from point samples[END_REF][START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF], which alters the definition of the restricted Delaunay complex.

We present background and foundational material in Section 2. Then, in Section 3, we exploit results established in [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] to demonstrate sampling conditions under which the intrinsic Delaunay complex, the restricted Delaunay complex, and the tangential Delaunay complex coincide and are manifold. The algorithm itself is presented in Section 4, and the analysis of the algorithm is presented in Section 5.

Background

Within the context of the standard m-dimensional Euclidean space R m , when distances are determined by the standard norm, • , we use the following conventions. The distance between a point p and a set X ⊂ R m , is the infimum of the distances between p and the points of X, and is denoted d R m (p, X). We refer to the distance between two points a and b as ba or d R m (a, b) as convenient. A ball B R m (c, r) = {x | xc < r} is open, and B R m (c, r) is its topological closure. Generally, we denote the topological closure of a set X by X, the interior by int(X), and the boundary by ∂X. The convex hull is denoted conv(X), and the affine hull is aff(X).

We will make use of other metrics besides the Euclidean one. A generic metric is denoted d, and the associated open and closed balls are B(c, r), and B(c, r). If a specific metric is intended, it will be indicated by a subscript, for example in Section 3 we introduce d M , the intrinsic metric on a manifold M, which has associated balls B M (c, r).

If A is a k × j matrix, we denote its i th singular value by s i (A). We use the operator norm A = s 1 (A) = sup x =1 Ax .

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

sin ∠(U, V ) = sup u∈U u -π V u ,
where π V is the orthogonal projection onto V . This is the largest principal angle between U and V . The angle between affine subspaces K and H is defined as the angle between the corresponding parallel vector subspaces.

Sampling parameters and perturbations

The structures of interest will be built from a finite set P ⊂ R m , which we consider to be a set of sample points. If D ⊂ R m is a bounded set, then P is an -sample set for D if d R m (x, P) < for all x ∈ D. We say that is a sampling radius for D satisfied by P. If no domain D is specified, we say P is an -sample set if d(x, P ∪ ∂conv(P)) < for all x ∈ conv(P). Equivalently, P is an -sample set if it satisfies a sampling radius for

D (P) = {x ∈ conv(P) | d R m (x, ∂conv(P)) ≥ }.
In particular, if P is an -sample set for U , and P = U ∩ P, and conv(P) ⊂ U , then P is an -sample set.

A set P is λ-sparse if d R m (p, q) > λ for all p, q ∈ P. We usually assume that the sparsity of a -sample set is proportional to , thus: λ = µ 0 .

We consider a perturbation of the points P ⊂ R m given by a function ζ : P → R m . If ζ is such that d R m (p, ζ(p)) ≤ ρ, we say that ζ is a ρ-perturbation. As a notational convenience, we frequently define P = ζ(P), and let p represent ζ(p) ∈ P. We will only be considering ρ-perturbations where ρ is less than half the sparsity of P, so ζ : P → P is a bijection.

Points in P which are not on the boundary of conv(P) are interior points of P.

Simplices

Given a set of j + 1 points {p 0 , . . . , p j } ⊂ P ⊂ R m , a (geometric) j-simplex σ = [p 0 , . . . , p j ] is defined by the convex hull: σ = conv({p 0 , . . . , p j }). The points p i are the vertices of σ. Any subset {p i0 , . . . , p i k } of {p 0 , . . . , p j } defines a k-simplex τ which we call a face of σ. We write τ ≤ σ if τ is a face of σ, and τ < σ if τ is a proper face of σ, i.e., if the vertices of τ are a proper subset of the vertices of σ.

The boundary of σ, is the union of its proper faces: ∂σ = τ <σ τ . In general this is distinct from the topological boundary defined above, but we denote it with the same symbol. The interior of σ is int(σ) = σ \ ∂σ. Again this is generally different from the topological interior. Other geometric properties of σ include its diameter (the length of its longest edge), ∆(σ), and the length of its shortest edge, L(σ). If σ is a 0-simplex, we define L(σ) = ∆(σ) = 0.

Inria

For any vertex p ∈ σ, the face oppposite p is the face determined by the other vertices of σ, and is denoted σ p . If τ is a j-simplex, and p is not a vertex of τ , we may construct a (j + 1)simplex σ = p * τ , called the join of p and τ . It is the simplex defined by p and the vertices of τ , i.e., τ = σ p .

Our definition of a simplex has made an important departure from standard convention: we do not demand that the vertices of a simplex be affinely independent. A j-simplex σ is a degenerate simplex if dim aff(σ) < j. If we wish to emphasise that a simplex is a j-simplex, we write j as a superscript: σ j ; but this always refers to the combinatorial dimension of the simplex.

If σ is non-degenerate, then it has a circumcentre, C(σ), which is the centre of the smallest circumscribing ball for σ. The radius of this ball is the circumradius of σ, denoted R(σ). A degenerate simplex may or may not have a circumcentre and circumradius. We write R(σ) to indicate that a simplex has a circumcentre. We will make use of the affine space N (σ) composed of the centres of the balls that circumscribe σ. We sometimes refer to a point c ∈ N (σ) as a centre for σ. The space N (σ) is orthogonal to aff(σ) and intersects it at the circumcentre of σ.

Its dimension is m -dim aff(σ). The altitude of p in σ is D(p, σ) = d R m (p, aff(σ p )).
A poorly-shaped simplex can be characterized by the existence of a relatively small altitude. The thickness of a j-simplex σ is the dimensionless quantity

Υ(σ) = 1 if j = 0 min p∈σ D(p,σ) j∆(σ) otherwise. We say that σ is Υ 0 -thick, if Υ(σ) ≥ Υ 0 . If σ is Υ 0 -thick,
then so are all of its faces. Indeed if τ ≤ σ, then the smallest altitude in τ cannot be smaller than that of σ, and also ∆(τ ) ≤ ∆(σ).

Although he worked with volumes rather than altitudes, Whitney [Whi57, p. 127] proved that the affine hull of a thick simplex makes a small angle with any hyperplane which lies near all the vertices of the simplex. We can state this [BDG12, Lemma 2.5] as: Lemma 2.1 (Whitney angle bound) Suppose σ is a j-simplex whose vertices all lie within a distance η from a k-dimensional affine space, H ⊂ R m , with k ≥ j. Then sin ∠(aff(σ), H) ≤ 2η Υ(σ)∆(σ) .

Simplex perturbation

We will make use of two results displaying the robustness of thick simplices with respect to small perturbations of their vertices. The first observation bounds the change in thickness itself under small perturbations:

Lemma 2.2 (Thickness under perturbation) Let σ = [p 0 , . . . , p j ] and σ = [p 0 , . . . , pj ] be j-simplices such that pip i ≤ ρ for all i ∈ {0, . . . , j}. For any positive η ≤ 1, if

ρ ≤ (1 -η)Υ(σ) 2 L(σ) 14 , (1) then
D(p i , σ) ≥ ηD(p i , σ),
for all i ∈ {0, . . . , j}. It follows that

Υ(σ)∆(σ) ≥ ηΥ(σ)∆(σ), RR n°8273
Boissonnat, Dyer, & Ghosh and

Υ(σ) ≥ 1 - 2ρ ∆(σ) ηΥ(σ) ≥ 6 7 ηΥ(σ).

Proof

Let p, q ∈ σ with p, q the corresponding vertices of σ. Let v = pq and ṽ = pq. Define θ = ∠(v, aff(σ p )) and θ = ∠(ṽ, aff(σ p)). Since Υ(σ) ≤ Υ(σ p ), Whitney's Lemma 2.1 lets us bound ∠(aff(σ p ), aff(σ p)) by the angle α defined by

sin α = 2ρ Υ(σ)∆(σ)
.

Also, by an elementary geometric argument,

sin γ = 2ρ v
defines γ as an upper bound on the angle between the lines generated by v and ṽ.

Thus we have

D(p, σ) = ṽ sin θ ≥ ( v -2ρ) sin(θ -α -γ).
Using the addition formula for sine together with the facts that for x, y ∈ [0, π 2 ], (1x) ≤ cos x; 2 sin x ≥ x; and sin x + sin y ≥ sin(x + y), we get

D(p, σ) ≥ ( v -2ρ) 1 -2 2ρ Υ(σ)∆(σ) + 2ρ v D(p, σ) v - 2ρ Υ(σ)∆(σ) + 2ρ v . For convenience, define µ = 2ρ L(σ) ≥ 2ρ v ≥ 2ρ ∆(σ) . Then D(p, σ) ≥ v (1 -µ) 1 -2 1 + 1 Υ(σ) µ D(p, σ) v -1 + 1 Υ(σ) µ ≥ (1 -µ) 1 - 4µ Υ(σ) D(p, σ) - 2µ v Υ(σ) ≥ (1 -µ) 1 - 4µ Υ(σ) D(p, σ) - 2µ v Υ(σ) 2 ∆(σ) D(p, σ) ≥ (1 -µ) 1 - 4µ Υ(σ) - 2µ Υ(σ) 2 D(p, σ) ≥ 1 - 7µ Υ(σ) 2 D(p, σ) ≥ KD(p, σ) when µ ≤ (1 -K)Υ(σ) 2 7 .
The condition on µ is satisfied when ρ satisfies Inequality (1). The bound on Υ(σ)∆(σ) follows immediately from the bounds on the D(p, σ), and the bound on Υ(σ) itself follows from the observation that

∆(σ) ∆(σ) ≥ ∆(σ) ∆(σ) + 2ρ ≥ 1 - 2ρ ∆(σ) ≥ 1 - Υ(σ) 2 7 ≥ 6 7 ,
when ρ satisfies Inequality (1).

We will also make use of a bound relating circumscribing balls of a simplex that undergoes a perturbation: Inria Lemma 2.3 (Circumscribing balls under perturbation) Let σ = [p 0 , . . . , p j ] and σ = [p 0 , . . . , pj ] be j-simplices such that pip i ≤ ρ for all i ∈ {0, . . . , j}. Suppose B = B R m (c, r), with r < , is a circumscribing ball for σ. If

ρ ≤ Υ(σ) 2 L(σ) 28 ,
then there is a circumscribing ball

B = B R d (c, r) for σ with c -c < 8 Υ(σ)∆(σ) ρ (2)
and

|r -r| < 9 Υ(σ)∆(σ) ρ.
If, in addition, we have that p0 = p 0 , then |r -r| ≤ cc , and (2) serves also as a bound on |r -r|.

Proof

By the perturbation bounds, the distances between c and the vertices of σ differ by no more than 2ρ. Also, cp i < ˜ = + ρ, and so by [BDG12, Lemma 4.3] we have

d R m (c, N (σ)) < 2˜ ρ Υ(σ)∆(σ)
.

The bound on ρ allows us to apply Lemma 2.2 with K = 1 2 , so Υ(σ)∆(σ) ≥ 1 2 Υ(σ)∆(σ), and we obtain the bound on cc with the observation that ˜ ≤ 2 . Indeed, ρ ≤ because L(σ) ≤ 2 .

By the triangle inequality |r -r| ≤ p0p 0 + cc , and the stated bound on |r -r| follows from the observation that Υ(σ)∆(σ) ≥ 1 if j > 1. Under the assumption that p0 = p 0 , the bound on cc also serves as a bound on |r -r|.

Flakes

For algorithmic reasons, it is convenient to have a more structured constraint on simplex geometry than that provided by a simple thickness bound. A simplex that is not thick has a relatively small altitude, but we wish to exploit a family of bad simplices for which all the altitudes are relatively small. As shown by Lemma 2.7 below, the Γ 0 -flakes form such a family. The flake parameter Γ 0 is a positive real number smaller than one. Definition 2.4 (Γ 0 -good simplices and Γ 0 -flakes) A simplex σ is Γ 0 -good if Υ(σ j ) ≥ Γ j 0 for all j-simplices σ j ≤ σ. A simplex is Γ 0 -bad if it is not Γ 0 -good. A Γ 0 -flake is a Γ 0 -bad simplex in which all the proper faces are Γ 0 -good.

Observe that a flake must have dimension at least 2, since Υ(σ j ) = 1 for j < 2. A flake that has an upper bound on the ratio of its circumradius to its shortest edge is called a sliver . The flakes we will be considering have no upper bound on their circumradius, and in fact they may be degenerate and not even have a circumradius.

Ensuring that all simplices are Γ 0 -good is the same as ensuring that there are no flakes. Indeed, if σ is Γ 0 -bad, then it has a j-face σ j ≤ σ that is not Γ j 0 -thick. By considering such a face with minimal dimension we arrive at the following important observation: Lemma 2.5 A simplex is Γ 0 -bad if and only if it has a face that is a Γ 0 -flake.

We obtain an upper bound on the altitudes of a Γ 0 -flake through a consideration of dihedral angles. In particular, we observe the following general relationship between simplex altitudes: Boissonnat, Dyer, & Ghosh Lemma 2.6 If σ is a j-simplex with j ≥ 2, then for any two vertices p, q ∈ σ, the dihedral angle between σ p and σ q defines an equality between ratios of altitudes: sin ∠(aff(σ p ), aff(σ q )) = D(p, σ) D(p, σ q ) = D(q, σ) D(q, σ p ) .

Proof

Let σ pq = σ p ∩σ q , and let p * be the projection of p into aff(σ pq ). Taking p * as the origin, we see that p-p * D(p,σq) has the maximal distance to aff(σ p ) out of all the unit vectors in aff(σ q ), and this distance is D(p,σ) D(p,σq) . By definition this is the sine of the angle between aff(σ p ) and aff(σ q ). A symmetric argument is carried out with q to obtain the result.

We arrive at the following important observation about flake simplices: Lemma 2.7 (Flakes have small altitudes) If a k-simplex σ is a Γ 0 -flake, then for every vertex p ∈ σ, the altitude satisfies the bound

D(p, σ) < k∆(σ) 2 Γ 0 (k -1)L(σ) .

Proof

Recalling Lemma 2.6 we have

D(p, σ) = D(q, σ)D(p, σ q ) D(q, σ p ) ,
and taking q to be a vertex with minimal altitude, we have

D(q, σ) = kΥ(σ)∆(σ) < kΓ k 0 ∆(σ), and
D(q, σ p ) ≥ (k -1)Υ(σ p )∆(σ p ) ≥ (k -1)Γ k-1 0 L(σ), and
D(p, σ q ) ≤ ∆(σ q ) ≤ ∆(σ),
and the bound is obtained.

Complexes

Given a finite set P, an abstract simplicial complex is a set of subsets K ⊂ 2 P such that if σ ∈ K, then every subset of σ is also in K. The Delaunay complexes we study are abstract simplicial complexes, but their simplices carry a canonical geometry induced from the inclusion map ι : P → R m . (We assume ι is injective on P, and so do not distinguish between P and ι(P).) To each abstract simplex σ ∈ K, we have an associated geometric simplex conv(ι(σ)), and normally when we write σ ∈ K, we are referring to this geometric object. Occasionally, when it is convenient to emphasise a distinction, we will write ι(σ) instead of σ.

Thus we view such a K as a set of simplices in R m , and we refer to it as a complex , but it is not generally a (geometric) simplicial complex. A geometric simplicial complex is a finite collection G of non-degenerate simplices in R N such that if σ ∈ G, then all of the faces of σ also belong to G, and if σ, σ ∈ G and τ = σ ∩ σ = ∅, then τ ≤ σ and τ ≤ σ. An abstract simplicial complex is defined from a geometric simplicial complex in an obvious way. A geometric realization of an abstract simplicial complex K is a geometric simplicial complex whose associated abstract simplicial complex may be identified with K. A geometric realization always exists for any Inria complex. Details can be found in algebraic topology textbooks; the book by Munkres [START_REF] Munkres | Elements of Algebraic Topology[END_REF] for example.

The carrier of an abstract complex K is the underlying topological space |K|, associated with a geometric realization of K. Thus if G is a geometric realization of K, then |K| = σ∈G σ. For our complexes, the inclusion map ι induces a continous map ι : |K| → R m , defined by barycentric interpolation on each simplex. If this map is injective, we say that K is embedded . In this case ι also defines a geometric realization of K, and we may identify the carrier of K with the image of ι.

A subset K ⊂ K is a subcomplex of K if it is also a complex. The star of a subcomplex K ⊆ K is the subcomplex generated by the simplices incident to K . I.e., it is all the simplices that share a face with a simplex of K , plus all the faces of such simplices. This is a departure from a common usage of this same term in the topology literature. The star of K is denoted star(K ) when there is no risk of ambiguity, otherwise we also specify the parent complex, as in star(K ; K).

A triangulation of P ⊂ R m is an embedded complex K with vertices P such that |K| = conv(P).

Definition 2.8 (Triangulation at a point) A complex K is a triangulation at p ∈ R m if:

• p is a vertex of K.

• star(p) is embedded.

• p lies in int(|star(p)|).

• For all τ ∈ K, and σ ∈ star(p), if int(τ ) ∩ σ = ∅, then τ ∈ star(p).

A complex K is a j-manifold complex if the star of every vertex is isomorphic to the star of a triangulation of R j . If σ is a simplex with vertices in P, then any map ζ : P → P ⊂ R m defines a simplex ζ(σ) whose vertices in P are the images of vertices of σ. If K is a complex on P, and K is a complex on P, then ζ induces a simplicial map K → K if ζ(σ) ∈ K for every σ ∈ K. We denote this map by the same symbol, ζ. We are interested in the case when ζ is an isomorphism, which means it establishes a bijection between K and K. We then say that K and K are isomorphic, and write

K ∼ = K, or K ζ ∼ = K if
we wish to emphasise that the correspondence is given by ζ.

We use the following local version of a standard result [BDG12, Lemma 2.7]: Lemma 2.9 Suppose K is a complex with vertices P ⊂ R m , and K a complex with vertices P ⊂ R m . Suppose also that K is a triangulation at p ∈ P, and that ζ : P → P induces an injective simplicial map star(p

) → star(ζ(p)). If K is a triangulation at ζ(p), then ζ(star(p)) = star(ζ(p)).

The Delaunay complex

An empty ball is one that contains no point from P. Definition 2.10 (Delaunay complex) A Delaunay ball is a maximal empty ball. Specifically, B = B R m (x, r) is a Delaunay ball if any empty ball centred at x is contained in B. A simplex σ is a Delaunay simplex , if there exists some Delaunay ball B such that the vertices of σ belong to ∂B ∩ P. The Delaunay complex is the set of Delaunay simplices, and is denoted Del(P).

The Delaunay complex has the combinatorial structure of an abstract simplicial complex, but Del(P) is embedded only when P satisfies appropriate genericity requirements [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF].

Protection

A Delaunay simplex σ is δ-protected if it has a Delaunay ball B such that d R m (q, ∂B) > δ for all q ∈ P\σ. We say that B is a δ-protected Delaunay ball for σ. If τ < σ, then B is also a Delaunay ball for τ , but it cannot be a δ-protected Delaunay ball for τ . We say that σ is protected to mean that it is δ-protected for some unspecified δ > 0.

Definition 2.11 (δ-generic) A finite set of points P ⊂ R m is δ-generic if all the Delaunay m-simplices are δ-protected. The set P is simply generic if it is δ-generic for some unspecified δ > 0.

We have previously demonstrated [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] that δ-generic point sets impart a quantifiable stability on the Delaunay complex. In Section 3 we review the main stability result and develop it to define the sampling conditions that will be met by the algorithm that we introduce in Section 4.

The Delaunay complex in other metrics

We will also consider the Delaunay complex defined with respect to a metric d on R m which differs from the Euclidean one. Specifically, if P ⊂ U ⊂ R m and d : U × U → R is a metric, then we define the Delaunay complex Del d (P) with respect to the metric d.

The definitions are exactly analogous to the Euclidean case: A Delaunay ball is a maximal empty ball B(x, r) in the metric d. The resulting Delaunay complex Del d (P) consists of all the simplices which are circumscribed by a Delaunay ball with respect to the metric d. The simplices of Del d (P) are, possibly degenerate, geometric simplices in R m . As for Del(P), Del d (P) has the combinatorial structure of an abstract simplicial complex, but unlike Del(P), Del d (P) may fail to be embedded even when there are no degenerate simplices.

Obtaining Delaunay triangulations in other metrics

Delaunay [START_REF] Delaunay | Sur la sphère vide[END_REF] showed that if P ⊂ R m is generic, then Del(P) is a triangulation. Point sets that are not generic are often dismissed in theoretical work, because an arbitrarily small perturbation of the points can be made which will yield a generic point set. Thus in the sense of the standard measure in the configuration space R m×|P| , almost all point sets will yield a Delaunay triangulation. However, when the metric is no longer Euclidean, this is no longer true.

In contrast to the purely Euclidean case, topological problems arise in point sets that are "near degenerate" , i.e., point sets that are not δ-generic for a sufficiently large δ. How large δ needs to be depends on how much the metric differs from the Euclidean one. Indeed, this was the initial motivation for the introduction of δ-generic point sets [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF], which are central to the results presented in this paper.

As we show here with a qualitative argument, the problem can be viewed as arising from the fact that when m is greater than two, the intersection of two metric spheres is not uniquely specified by m points. We demonstrate the issue in the context of Delaunay balls. The problem is developed quantitatively in terms of the Voronoi diagram in Appendix A.

We work exclusively on a three dimensional domain, and we are not concerned with "boundary conditions"; we are looking at a coordinate patch on a densely sampled compact 3-manifold.

One core ingredient in Delaunay's triangulation result [START_REF] Delaunay | Sur la sphère vide[END_REF] is that any triangle τ is the face of exactly two tetrahedra. This follows from the observation that a triangle has a unique circumcircle, and that any circumscribing sphere for τ must include this circle. The affine hull of τ cuts space into two components, and if τ ∈ Del(P), then it will have an empty circumsphere centred at a point c on the line through the circumcentre and orthogonal to aff(τ ). The point c is contained on an interval on this line which contains all the empty spheres for τ . The endpoints of the interval are the circumcentres of the two tetrahedra that share τ as a face.

The argument hinges on the assumption that the points are in general position, and the uniqueness of the circumcircle for τ . If there were a fourth vertex lying on that circumcircle, then there would be three tetrahedra that have τ as a face, but this configuration would violate the assumption of general position. Now if we allow the metric to deviate from the Euclidean one, no matter how slightly, the guarantee of a well defined unique circumcircle for τ is lost. In particular, If three spheres S 1 , S 2 and S 3 all circumscribe τ , their pairwise intersections will be different in general. I.e.,

S 1 ∩ S 3 = S 2 ∩ S 3 .
Although these intersections may be topological circles that are "arbitrarily close" assuming the deviation of the metric from the Euclidean one is small enough, "arbitrarily close" is not good enough when the only genericity assumption allows configurations that are arbitrarily bad.

An attempt to illustrate the problem is given in Figure 1, where τ = {u, v, p}. Here, the sphere S 3 would be contained inside the spheres S 1 and S 2 if the metric were Euclidean, but any aberration in the metric may leave a part of S 3 exposed to the outside. This means that in principle another sample point w could lie on S 3 , while S 1 and S 2 remain empty. Thus there are three tetrahedra that share τ as a face.

The essential difference between dimension 2 and the higher dimensions can be observed by examining the topological intersection properties of spheres. Specifically, two (m -1)-spheres intersect transversely in an (m -2)-sphere. For a non-Euclidean metric, even if this property holds for sufficently small geodesic spheres, only in dimension two is the sphere of intersection of the Delaunay spheres of two adjacent m-simplices uniquely determined by the vertices of the shared (m -1)-simplex. See Figure 1.

Boissonnat, Dyer, & Ghosh

The Voronoi diagram

We will occasionally make reference to the Voronoi diagram, which is a structure dual to the Delaunay complex. It offers an alternative way to interpret observations made with respect to the Delaunay complex.

The Voronoi cell associated with p ∈ P with respect to the metric d : U × U → R is given by

V d (p) = {x ∈ U | d(x, p) ≤ d(x, q) for all q ∈ P}.
More generally, a Voronoi face is the intersection of a set of Voronoi cells: given {p 0 , . . . , p k } ⊂ P, let σ denote the corresponding abstract simplex. We define the associated Voronoi face as

V d (σ) = k i=0 V d (p i ).
It follows that σ is a Delaunay simplex if and only if V d (σ) = ∅. In this case, every point in

V d (σ)

Equating Delaunay structures

We now turn to the task of triangulating M, a smooth, compact m-manifold, without boundaries embedded in R N . In this section we demonstrate our main structural result, Theorem 3.5, which is stated at the end of Section 3.1. It says that the complex constructed by the algorithm we describe in Section 4 is in fact an intrinsic Delaunay triangulation of the manifold, which we introduce next.

Delaunay structures on manifolds

The restricted Delaunay complex is the Delaunay complex Del R N | M (P) obtained when distances on the manifold are measured with the metric d R N | M . This is the Euclidean metric of the ambient space, restricted to the submanifold M. In other words,

d R N | M (x, y) = d R N (x, y).
We use this notation to avoid ambiguities in conjunction with the local Euclidean metrics discussed below. The Delaunay complex Del R N | M (P) is a substructure of Del R N (P).

Alternatively, distances on the manifold may be measured with d M , the intrinsic metric of the manifold. This metric defines the distance between x and y as the infimum of the lengths of the paths on M which connect x and y. Since the length of a path on M is defined as its length as a curve in R N , this metric is also induced from d R N . The intrinsic Delaunay complex is the Delaunay structure Del M (P) associated with this metric.

Although neither of these metrics are Euclidean, the idea is that locally, in a small neighbourhood of any point, these metrics may be well approximated by d R m . Then, if the sampling satisfies appropriate δ-generic and -dense criteria in these local Euclidean metrics, the global Delaunay complex in the metric of the manifold will coincide locally with a Euclidean Delaunay triangulation, and we can thus guarantee a manifold complex.

Local Euclidean metrics

A local coordinate chart at a point p ∈ M, is a pair (W, φ p ), where W ⊂ M is an open neighbourhood of p, and φ p : W → U = φ p (W ) ⊂ R m is a homeomorphism onto its image, with Inria φ p (p) = 0. A local coordinate chart allows us to pull back the Euclidean metric to W . For all x, y ∈ W , the metric d φp (x, y) = d R m (φ p (x), φ p (y)) is a a local Euclidean metric for p on W . This metric depends upon the choice of φ p ; there are different ways to impose a Euclidean metric on W .

It is convenient to take the reciprocal point of view, and work with a local parameterization at a point p ∈ M. This is a pair (U, ψ p ), such that U ⊂ R m , and (W, ψ -1 p ) is a local coordinate chart for p, where W = ψ p (U ). We can then use ψ p to pull back the metric of the manifold to U , and to simplify the notation we write d M (x, y) for x, y ∈ U , where it is to be understood that this means d M (ψ p (x), ψ p (y)), and likewise for d R N | M (x, y). Indeed, once W and U have been coupled together by a homeomorphism, we can transfer the metrics between them and the distinction becomes only one of perspective; the standard metric d R m on U is a local Euclidean metric for p.

We wish to generate a sample set P ⊂ M that will allow us to exploit the stability results for Delaunay triangulations [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]. We consider the stability of a Delaunay triangulation in a local Euclidean metric. The following definition is convenient when stating the stability results:

Definition 3.1 (Secure simplex) A simplex σ ∈ Del(P) is secure if it is a δ-protected m- simplex that is Υ 0 -thick and satisfies R(σ) < and L(σ) ≥ µ 0 .
We will make reference to the following result [BDG12, Theorem 4.17]:

Theorem 3.2 (Metric stability assuming thickness) Suppose conv(P) ⊆ U ⊂ R m and the metric d : U × U → R is such that |d(x, y) -d R m (x, y)| ≤ ρ for all x, y ∈ U . Suppose also that P J ⊆ P is such that every m-simplex σ ∈ star(P J ; Del(P)) is secure and satisfies d R m (p, ∂U ) ≥ 2 for every vertex p ∈ σ. If ρ ≤ Υ 0 µ 0 36 δ,
then star(P J ; Del d (P)) = star(P J ; Del(P)).

In our context the point set P used in Theorem 3.2 will come from a larger point set P, such that P = W ∩ P. We will write P W in order to emphasise this dependence on W . We want to ensure that star(P J ; Del d (P W )) = star(P J ; Del d (P)).

(

) 3
This requirement is attained by demanding that P satisfy a sampling radius of with respect to the metric d M . Since d R m (x, y) ≤ d M (x, y) for all x, y ∈ U ∼ = W , by our particular choice of ψ p , we will have that P W is an -sample set with respect to the metric d R m . We ensure that U is large enough so that d R m (p, ∂U ) ≥ 4 for all p ∈ P J . It then follows that R(σ) < for any simplex σ ∈ star(P J ; Del(P W )), because P W is an -sample set [BDG12, Lemma 3.6], and thus d R m (q, ∂U ) ≥ 2 for any q ∈ σ. It follows that d M (q, ∂U ) ≥ 2 as well, and thus the sampling radius on P ensures that Equation (3) is satisfied. For our purposes P J will consist of a single point p, and the sampling radius is constrained by the requirement that U be small enough that the metric distortion introduced by ψ p meets the requirements of Theorem 3.2.

The tangential Delaunay complex

The algorithm we describe in Section 4 is a variation of the algorithm described by Boissonnat and Ghosh [START_REF] Boissonnat | Triangulating smooth submanifolds with light scaffolding[END_REF]. This algorithm builds the tangential Delaunay complex , which we denote by Del T M (P). This is not a Delaunay complex as we have defined them, since it cannot be defined by the Delaunay empty ball criteria with respect to any single metric. However, it is a Delaunaytype structure, and as with Del R N | M (P), the tangential Delaunay complex is a substructure of

Boissonnat, Dyer, & Ghosh

Del R N (P). We will demonstrate sampling conditions which ensure that Del T M (P) = Del M (P) = Del R N | M (P).

Definition 3.3 (Tangential Delaunay complex) The tangential Delaunay complex for P ⊂ M ⊂ R N is defined by the criterion that σ ∈ Del T M (P) if it has an empty circumscribing ball B R N (c, r) such that c ∈ T p M for some vertex p ∈ σ.

We define some local complexes to facilitate discussions of the tangential Delaunay complex. For all p ∈ P, let

K(p) = {σ | V R N (σ) ∩ T p M = ∅},
and define star(p) = star(p; K(p)).

(4)

Then the tangential Delaunay complex is the union of the complexes star(p) for all p ∈ P.

Boissonnat et al. [BG11, Lemma 2.3] showed that V R N (P)∩T p M is equal to the m-dimensional weighted Voronoi diagram of P ⊂ T p M
, where P is the orthogonal projection of P onto T p M and the squared weight of a point p i ∈ P isp ip i 2 . Therefore, K(p) is isomorphic to a dual complex (the nerve) of the k-dimensional weighted Voronoi diagram of P .

Power protection

The algorithm introduced in Section 4.2 will ensure that for every simplex σ in the tangential Delaunay complex, and every vertex p ∈ σ, there is a Delaunay ball for σ that is centred on T p M and is protected in the following sense:

Definition 3.4 (Power protection) A simplex σ with Delaunay ball B R N (C, R) is δ2 -power- protected if d R N (C, q) 2 -R 2 > δ2 for all q ∈ P \ σ.
Observe that, if C ∈ M, the ball B R N (C, R) is not an object that can be described by the metric d R N | M . In the context of the tangential Delaunay complex we use power-protection rather than the protection described in Section 2.4.1 because working with squared distances is convenient when we consider the Delaunay complex restricted to an affine subspace.

Main structural result

The rest of Section 3 is devoted to the proof of Theorem 3.5 below. It says that for the point set generated by our algorithm, the tangential Delaunay complex is isomorphic with the intrinsic Delaunay complex of M. It then follows, from a previously published result [BG10, Theorem 5.1], that the intrinsic Delaunay complex is in fact homeormorphic to M; it is an intrinsic Delaunay triangulation.

Thus we obtain a partial recovery of the kind of results attempted by Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF]. Our sampling conditions, and our algorithm (existence proof) rely on the embedding of M in R N ; we leave purely intrinsic sampling conditions for future work.

Theorem 3.5 (Intrinsic Delaunay triangulation) Suppose P ⊂ M is (μ 0 )-sparse with respect to d R N , and every m-simplex σ ∈ Del T M (P) is Υ0 -thick, and has, for every vertex p ∈ σ, a δ2 -power-protected empty ball of radius less than centred on T p M, with δ ≥ δ 0 μ0 . If δ 2 0 μ2 0 ≤ 1 7 , and

≤ Υ2 0 μ3 0 δ 2 0 rch(M) 1.5 × 10 6 , then Del T M (P) = Del R N | M (P) = Del M (P),
Inria and for sufficiently small, these will be homeomorphic to M:

|Del M (P)| ∼ = M.

Choice of local Euclidean metric

A local parameterization at p ∈ M will be constructed with the aid of the orthogonal projection

π p : R N → T p M, (5)
restricted to M. As shown in Lemma B.4, Niyogi et al. [NSW08, Lemma 5.4] demonstrated that if r < rch(M) 2 , then π p is a diffeomorphism from W = B R N | M (p, r) onto its image U ⊂ T p M.
We will identify T p M with R m , and define the homeomorphism

ψ p = π p | -1 W : U -→ W. (6)
Using ψ p to pull back the metrics d M and d R N | M to R m , we can view them as perturbations of d R m . The magnitude of the perturbation is governed by the radius of the ball used to define W .

Definition 3.6 We call a neighbourhood W of p ∈ M admissible if W ⊆ B R N | M (p, r), with r ≤ rch(M)
100 . In all that follows, any mention of a local Euclidean metric refers to the one defined by π p restricted to an admissible neighbourhood. The requirement r ≤ rch(M) 100 is simply a convenient bound that yields a small integer constant in the perturbation bound of the following lemma, and does not constrain subsequent results. The bound could be relaxed to r ≤ rch(M) 4 at the expense of a weaker bound on the perturbation.

Lemma 3.7 (Metric distortion) Suppose (U, ψ p ) is a local parameterisation at p ∈ W ⊂ M with W = ψ p (U ). If W ⊆ B R N | M (p, r), with r ≤ rch(M)
100 , then for all x, y ∈ U ,

d R N | M (x, y) -d R m (x, y) ≤ |d M (x, y) -d R m (x, y)| ≤ 23r 2 rch(M) . Proof Let u, v ∈ W ⊂ B R N | M (p, r)
, and let θ be the angle between the line segments [u, v] and [π p (u), π p (v)], θ 1 the angle between [u, v] and T u M, and θ 2 the angle between T p M and T u M.

Thus θ ≤ θ 1 + θ 2 , and d R m (π p (u), π p (v)) = d R N (u, v)cos θ. Defining η = r rch(M) , Lemma B.5 yields d M (u, v) ≤ d R N (u, v) (1 + 4η) , (7)
and so

d R m (π p (u), π p (v)) ≥ d M (u, v) cos θ 1 + 4η .
Using Lemma B.1, we find sin θ 1 ≤ η, and Lemma B.3, yields sin θ 2 ≤ 6η. Therefore, since sin θ ≤ sin θ 1 + sin θ 2 , we have

cos θ = (1 -sin 2 θ) 1/2 ≥ 1 -sin θ ≥ 1 -7η and we get d R m (π p (u), π p (v)) ≥ d M (u, v) 1 -7η 1 + 4η ≥ d M (u, v) (1 -7η) (1 -4η) ≥ d M (u, v)(1 -11η). Using Equation (7) we find d M (u, v) ≤ 208r 100 , so d R m (π p (u), π p (v)) ≥ d M (u, v) -23 r 2 rch(M)
, and the result follows since

d M (u, v) ≥ d R N | M (u, v) ≥ d R m (π p (u), π p (v)).

Boissonnat, Dyer, & Ghosh

Our sampling radius is constrained by the size of a Euclidean ball that can be contained in an admissible neighbourhood. The following lemma gives a convenient expression for this:

Lemma 3.8 If 1 < a ≤ 10 4 and a ≤ rch(M)
100 , and

U = B R m (p, (a -1) ), then ψ p (U ) = W ⊆ B R N | M (p, a ). Proof Using Lemma B.1, we have that B R m (p, r) ⊆ π p (B R N | M (p, a )) if r 2 ≤ a 2 2 - a 2 2 2rch(M) 2 = a 2 2 1 - a 2rch(M) 2 ≤ a 2 2 1 - 1 200 2 .
Thus we require r ≤ 200 2 -1 200 2 a , which is satisfied by r = (a -1) if a ≤ 79999.

Lemmas 3.7 and 3.8 lead to a sampling radius which allows us to employ Theorem 3.2, and so obtain an equivalence between Delaunay structures:

Proof

As usual, let U = π p (W ). Then by Lemma 3.8 B R m (p, 4 ) ⊆ U , and thus d R m (q, ∂U ) ≥ 2 for any vertex q of a simplex in star(p; Del(P W )). Thus Lemma 3.7 allows us to apply Theorem 3.2 provided 23a 2 2 rch(M) ≤ Υ 0 µ 0 ν 0 36 , when a = 5, and we obtain the required bound on . Thus the star of every vertex in Del M (P) is equal to the star of that point in the local Euclidean metric, and likewise for Del R N | M (P). The claim follows since σ ∈ Del M (P) if and only if it is in the local Euclidean Delaunay triangulation of every one of its vertices, and likewise for the simplices in Del R N | M (P).

The protected tangential complex

We obtain Theorem 3.5 by means of Theorem 3.2 via the observation that power protection of the ambient Delaunay balls translates into protection in the local Euclidean metrics. We must distinguish between the geometry of a simplex defined with respect to the Euclidean metric d R N of the ambient space, as opposed to a local Euclidean metric d R m . In general, we use a tilda to indicate simplices in the ambient space, and their properties.

Inria

Lemma 3.10 (Protection under projection) Suppose P ⊂ M and that σ ∈ Del R N (P) is an Υ0 -thick m-simplex, with L(σ) ≥ μ0 and B R N (C, R) is a δ2 -power-protected empty ball for σ, with respect to the metric d R N , where δ2 ≥ δ 2 0 μ2 0 2 . Suppose also that C ∈ T p M, for some vertex p ∈ σ.

If R < , with

≤ Υ2 0 μ3 0 δ 2 0 rch(M) 512 , (9)
ν 0 = δ 2 0 μ2 0 8 . ( 10
)

Proof

We first find a bound for d R m (C, c) and r. Let σ = [p 0 , . . . , pm ], and σ = [p 0 , . . . , p m ] so that π p (p i ) = p i , and p = p 0 = p0 . We will first show that, near C, there is a circumcentre c for σ in the metric d R m . For any p i ∈ σ, d R N (p, p i ) < 2R, and so by Lemma B.1 we have

d R N (p i , p i ) ≤ 2R 2 rch(M) < 2 2 rch(M)
.

In order to apply Lemma 2.3 we require

2 2 rch(M) ≤ Υ2 0 μ0
28 , or

≤ Υ2 0 μ0 rch(M) 56 ,
which is satisfied by Equation (9). Since aff(σ) = T p M, the circumcentre c ∈ T p M is the closest point in N (σ) to C, Lemma 2.3 yields

|R -r| ≤ d R m (C, c) = d R N (C, c) < 16 2 Υ0 μ0 rch(M)
.

Now we seek a lower bound on the protection of B R m (c, r). Suppose q ∈ P \ σ. We wish to establish a lower bound on d R m (c, q)r, where q = π p (q). We may assume that d R N (C, q) < 3 , since otherwise q will lie outside of our region of interest.

Let z = (3 ) 2 2rch(M) be the upper bound on d R N (q, q) given by Lemma B.1.

Then d R m (C, q) 2 ≥ d R N (C, q) 2 -z 2 > R 2 + δ2 -z 2 . Thus d R m (C, q) -R > δ2 -z 2 d R m (C, q) + R > δ2 -z 2 4 , since R < . Then d R m (c, q) -r ≥ (d R m (C, q) -d R m (C, c)) -(R + |R -r|) > δ2 -z 2 4 -2d R m (C, c). Putting this together, using δ2 ≥ δ 2 0 μ2 0 2 , we get d R m (c, q) -r > 1 4 δ 2 0 μ2 0 - 81 2 16rch(M) 2 -

Υ0 μ0 rch(M)

.

In order to simplify away the final term, we demand

32 Υ0 μ0 rch(M) ≤ 1 16 δ 2 0 μ2 0 ,
which is satisfied by Equation (9). Under this constraint, the central term is also seen to be less than 1 16 δ 2 0 μ2 0 , and we obtain δ ≥ 1 8 δ 2 0 μ2 0 .

Proposition 3.9 requires a thickness Υ 0 and shortest edge bound µ 0 for the simplex σ ⊂ R m , but Lemma 3.10 is expressed in terms of the corresponding quantities Υ0 and μ0 for the corresponding simplex σ ⊂ R N .

Lemma 3.11 (Simplex distortion under projection) Let σ ∈ Del R N (P) be an m-simplex as described in Lemma 3.10, and let σ = π p (σ) be its projection in the local Euclidean metric for p on any admissible neighbourhood that contains B R N | M (p, 2 ), where p is a vertex of σ. If satisfies Equation (9), and δ 2 0 μ2 0 ≤ 1 7 , then L(σ) > µ 0 , where

µ 0 = 20 21 μ0 ,
and Υ(σ) > Υ 0 , where

Υ 0 = 6 49 Υ0 . Proof Since π p (B R N | M (p, 2 )) ⊆ B R m (p, 2 )
, it is sufficient to apply the Metric distortion lemma 3.8 with a = 3.

For the shortest edge length, we find

L(σ) ≥ L(σ) - 3 2 × 23 2 rch(M) = μ0 1 - 207 Υ2 0 δ 2 0 μ2 0 512 > μ0 1 - Υ2 0 δ 2 0 μ2 0 3 > 20 21 μ0 .
For the thickness bound, in order to apply Lemma 2.2, using η = (1η), we require

207 Υ2 0 δ 2 0 μ3 0 512 ≤ η Υ2 0 μ0 14 , which is satisfied if we choose η > 6δ 2 0 μ2 0 . Then Lemma 2.2 yields Υ(σ) ≥ 6 7 1 -6δ 2 0 μ2 0 Υ(σ) > 6 49 Υ(σ).
We can now express the sampling conditions in terms of the output parameters of the tangential complex algorithm, and this allows us to apply Proposition 3.9 and obtain our main structural result: Inria Proof of Theorem 3.5 We first translate the sampling requirements of Proposition 3.9 in terms of properties of simplices in the ambient metric d R N . Using Lemma 3.11, together with Equation (10), the upper bound on the sampling radius demanded by Proposition 3.9 becomes

≤ 20 × 6 Υ0 δ 2 0 μ3 0 rch(M) 21 × 49 × 8 × 20700
.

We obtain the stated sampling radius bound after multiplying by Υ0 in order to ensure that the demand of Equation ( 9) is also met. Thus the stated sampling radius satisfies the requirements of both Lemma 3.10 and Proposition 3.9. The fact that the structures are isomorphic follows from the fact that they are all locally isomorphic to the Delaunay triangulation in the local Euclidean metric. To see that star(p; Del T M (P)) ∼ = star(p; Del(P W )), observe that Lemma 3.10 implies that there is an injective simplicial map star(p; Del T M (P)) → star(p; Del(P W )). The isomorphism is established by Lemma 2.9, once it is established that star(p; Del T M (P)) is a triangulation at p. In fact star(p; Del T M (P)) is isomorphic to the star of p in a regular triangulation of the projected points P W ; it is a weighted Delaunay triangulation [BG10], and with our choice of W , the point p is an interior point in this triangulation [BG10, Lemma 2.7(1)]. Thus star(p; Del T M (P)) is a triangulation at p, and it follows that star(p; Del T M (P)) ∼ = star(p; Del(P W )).

The equality of the Delaunay complexes now follows from Proposition 3.9, Equation (8).

The homeomorphism assertion follows from previous work [BG10, Theorem 5.1].

Algorithm

In this section we introduce a Delaunay refinement algorithm which, while constructing a tangential Delaunay complex, will transform the input sample set into one which meets the requirements of Theorem 3.5. In particular we wish to construct a tangential Delaunay complex in which every m-simplex σ is Υ0 -thick and for every p ∈ σ, there is a δ2 -power-protected Delaunay ball for σ centred on T p M. We demand δ ≥ δ 0 μ0 , where provides a strict upper bound on the radius of these Delaunay balls, and μ0 provides a lower bound on the shortest edge length of any simplex in Del T M (P). The constants δ 0 and μ0 are both positive and smaller than one. The algorithm is in the same vein as that of Boissonnat and Ghosh [START_REF] Boissonnat | Triangulating smooth submanifolds with light scaffolding[END_REF], which is in turn an adaptation of the algorithm introduced by Li [START_REF] Li | Generating well-shaped d-dimensional Delaunay meshes[END_REF]. It is described in Section 4.2, after we introduce terminology and constructs which are used in the algorithm in Section 4.1.

Components of the algorithm

We now introduce the primary concepts that are used as building blocks of the algorithm.

Elementary weight functions

Elementary weight functions are a convenient device to facilitate the identification of simplices σ that are not δ2 -power-protected for δ = δ 0 L(σ).

In order to emphasise that we are considering a function defined only on the set of vertices of a simplex, we denote by σ the set {p 0 , . . . , p k } of vertices of σ = [p 0 . . . p k ]. We will call ω σ : σ → [0, ∞) an elementary weight function if it satisfies the following conditions:

RR n°8273

1. There exists p i ∈ σ such that ω σ (p i ) ∈ [0, δ 0 L(σ)], and 2. for all p j ∈ σ \ p i , ω σ (p j ) = 0.

For a given σ = [p 0 , . . . , p k ] and elementary weight function ω σ , we define N (σ, ω σ ) as the set of solutions to the following system of k equations:

x -p i 2 -x -p 0 2 = ω σ (p i ) 2 -ω σ (p 0 ) 2 .
In direct analogy with the space N (σ) of centres of σ, the set N (σ, ω σ ) is an affine space of dimension mdim aff(σ) that is orthogonal to aff(σ). We denote by C(σ, ω σ ) the unique point in N (σ, ω σ ) ∩ aff(σ), and we define

R(σ, ω σ ) 2 = p 0 -C(σ, ω σ ) 2 -ω σ (p 0 ) 2 ,
where the notation is chosen to emphasise the close relationship with the circumcentre C(σ) and circumradius R(σ). The following lemma exposes some properties of R(σ, ω σ ) in this spirit:

Lemma 4.1 For a given σ = [p 0 , . . . , p k ], with k ≥ 1, and elementary weight function ω σ , we have:

1. If σ 1 ≤ σ then ω σ1 = ω σ | σ1
is an elementary weight function, and

R(σ 1 , ω σ1 ) ≤ R(σ, ω σ ). 2. ∆(σ) ≤ 2 1-δ 2 0 R(σ, ω σ ). 3. If Υ(σ) > 0, then 1 -η ≤ R(σ, ω σ ) R(σ) ≤ 1 + η, with η = δ 2 0 Υ(σ) .
Proof 1. That ω σ1 is an elementary weight function follows from the observation that L(σ 1 ) ≥ L(σ). Since N (σ, ω σ ) ⊆ N (σ 1 , ω σ1 ), the projection of C(σ, ω σ ) into aff(σ 1 ) is C(σ 1 , ω σ1 ). The result then follows from the Pythagorean theorem.

2. Let e = [p 0 , p 1 ] be the longest edge of σ, and let c denote the projection of C(σ, ω σ ) onto aff(e). Without loss of generality we assume that ω(p 0 ) = 0.

We have

p 0 -c 2 = p 1 -c 2 -ω σ (p 1 ) 2 = (p 1 -p 0 ) -(c -p 0 ) 2 -ω σ (p 1 ) 2 = ∆(σ) 2 -2(p 1 -p 0 ) • (c -p 0 ) + p 0 -c 2 -ω σ (p 1 ) 2 .
Since p 0 , p 1 , and c are colinear, we have 2(p 1p 0 ) • (cp 0 ) = 2∆(σ) p 0c , and using the fact that ω σ (p 1 ) ≤ δ 0 L(σ), we get

p 0 -c = ∆(σ) 2 1 - ω σ (p 1 ) 2 ∆(σ) 2 ≥ (1 -δ 2 0 )∆(σ) 2 .

Inria

The result follows from the fact that R(σ, ω σ ) ≥ p 0c . 3. Using the fact that ω σ (p) = 0 for all vertices p ∈ σ, except at most one, we get p i ∈ ∂B R k (C(σ, ω σ ), R(σ, ω σ )) for all p i ∈ σ except at most one.

Let η = C(σ, ω σ ) -C(σ) , and assume, without loss of generality, that the vertex p 0 ∈ ∂B(C(σ, ω σ ), R(σ, ω σ )). Therefore,

C(σ) -p 0 -η ≤ C(σ, ω σ ) -p 0 ≤ C(σ) -p 0 + η R(σ) + η ≤ R(σ, ω σ ) ≤ R(σ) + η. (11)
Since the point in N (σ) that is closest to C(σ, ω σ ) is C(σ), and Υ(σ) > 0, we obtain the following bound using Lemma 4.1 from [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]:

η ≤ δ 2 0 L(σ) 2 2Υ(σ)∆(σ) ≤ δ 2 0 Υ(σ) R(σ), since L(σ) ≤ ∆(σ) ≤ 2R(σ). (12)
The result now follows from Eq. ( 11) and ( 12).

If σ = p * σ p , and ω σ is an elementary weight function that vanishes on σp , then N (σ, ω σ ) ⊆ N (σ p ), but no point in N (σ, ω σ ) can be the centre of a δ2 -power-protected Delaunay ball for σ p for any δ ≥ δ 0 L(σ). In other words, σ and ω σ define a quasi-cospherical configuration that is an obstruction to the power protection of σ p at all points in N (σ, ω σ ).

Quasicospherical configurations

We now define the family of simplices that our algorithm must eliminate in order to ensure that the final point set has the desired protection properties.

Recalling the definition (4) of star(p), we have the following [BG10, Lemma 2.7 (1)]:

Lemma 4.2 Let P ⊂ M satisfy a sampling radius of with respect to d R N such that ≤ rch(M)/16. Then for all x ∈ V R N (p) ∩ T p M, we have px ≤ 4 . In particular, for all p ∈ P, and every m-simplex σ ∈ star(p), we have R p (σ) ≤ 4 .

Since by Lemma 4.2, the Voronoi cell of p restricted to T p M is bounded, we get:

Lemma 4.3 If ≤ rch(M)
16 , then the combinatorial dimension of the maximal simplices in star(p) is at least m.

We will always assume that P satisfies a sampling radius of ≤ rch(M)

16 . If σ is a maximal simplex in star(p), then V R N (σ) intersects T p M at a single point. Indeed, since V R N (σ) ⊂ V R N (p), by Lemma 4.2 the convex set V R N (σ) ∩ T p M is bounded, and if it had a nonempty interior, then σ would not be maximal. Let σ be a maximal simplex in star(p). Then, for all σ m ≤ σ, the unique point in V R N (σ) ∩ T p M . will be denoted by c p (σ m ). We denote the radius of the circumscribing ball centred at c p (σ m ) by R p (σ m ), i.e., R p (σ m ) = pc p (σ m ) .

In our algorithm we will use the following complex, whose definition employs a particular elementary weight function:

cosph δ0 (p) = σ m+1 = p m+1 * σ m σ m ∈ star(p), R p (σ m ) < , σ m is Γ 0 -good, and ∃ ω σ m+1 with ω σ m+1 | σm = 0 and c p (σ m ) ∈ N (σ m+1 , ω σ m+1 ) . (13)
The (m + 1)-dimensional simplices in cosph δ0 (p) are analogous to inconsistent configurations defined in [START_REF] Boissonnat | Triangulating smooth submanifolds with light scaffolding[END_REF][START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF]. Unless otherwise stated, whenever σ m+1 = p m+1 * σ m ∈ cosph δ0 (p), with σ m ∈ star(p), the mention of ω σ m+1 will refer to the elementary weight function identified in Equation (13). In particular, ω σ m+1 (p i ) = 0 for all p i ∈ σm+1 \ p m+1 , and

ω σ m+1 (p m+1 ) ∈ [0, δ 0 L(σ m+1 )] satisfies c p (σ m ) -p 2 = c p (σ m ) -p m+1 2 -ω σ m+1 (p m+1 ) 2 .
We will exploit the following observations:

Lemma 4.4 If σ m+1 = p m+1 * σ m ∈ cosph δ0 (p) with σ m ∈ star(p), then R(σ m+1 , ω σ m+1 ) ≤ R p (σ m ) and ∆(σ m+1 ) ≤ 2 1 -δ 2 0 R p (σ m ) Proof Since c p (σ m ) ∈ N (σ m+1 , ω σ m+1 ), it follows that C(σ m+1 , ω σ m+1
) is the projection of c p (σ m ) into aff(σ m+1 ), and therefore R p (σ m ) ≥ R(σ m+1 , ω σ m+1 ). The bound on ∆(σ m+1 ) now follows directly from Lemma 4.1.

Boissonnat et al. [BG11]

, using Lemma 4.2, showed that we can compute star(p) by computing a weighted Delaunay triangulation on T p M of the points obtained by projecting P onto T p M. Once star(p) has been computed, we can compute cosph δ0 (p) by a simple distance computation.

The importance of cosph δ0 (p) lies in the observation that if an m-simplex σ m ∈ star(p) is not sufficiently power-protected, then there will be a simplex in cosph δ0 (p) that is a witness to this. It is a direct consequence of the definitions, but we state it explicitly for reference: Lemma 4.5 If P is μ0 -sparse, and cosph δ0 (p) = ∅, then every σ m ∈ star(p) is δ 2 0 μ2 0 2 -power protected on T p M.

Unfit configurations and the picking region

The refinement algorithm, at each step, kills an unfit configuration by inserting a new point x = ψ p (x ) where x belongs to the so-called picking region of the unfit configuration, and ψ p is the inverse projection defined in Equation (6). We use the term unfit configuration to distinguish the elements under consideration from other simplices. An unfit configuration φ may be one of two types:

Big configuration: An m-simplex φ = σ m in star(p) is a big configuration if R p (σ m ) ≥ .
Bad configuration: A simplex φ is a bad configuration if it is Γ 0 -bad and it is either an msimplex φ = σ m ∈ star(p) that is not a big configuration, or it is an (m + 1)-simplex φ = σ m+1 ∈ cosph δ0 (p).

We will show in Section 5.2, Lemma 5.13, that in fact every (m + 1)-simplex in cosph δ0 (p) is a bad configuration. The size of the picking region is governed by a positive parameter α < 1 called the picking ratio.

Inria

Definition 4.6 (Picking region) The picking region of a bad configuration, σ m ∈ star(p) or p m+1 * σ m ∈ cosph δ0 (p) with σ m ∈ star(p), denoted by P (σ m , p) and P (σ m+1 , p) respectively, is defined to be the m-dimensional ball

B R N (c p (σ m ), αR p (σ m )) ∩ T p M.
We choose a point in the picking region so as to minimize the introduction of new unfit configurations. We are able to avoid creating new bad configurations provided that the radius of the potential configuration is not too large. To this end, we introduce the parameter β > 1.

Definition 4.7 (Hitting sets and good points) Let φ = σ m ∈ star(p) or φ = q * σ m ∈ cosph δ0 (p) with σ m ∈ star(p), and x = ψ p (y) where y ∈ P (φ, p). A set σ ⊂ P of size k, with k ≤ m + 1, is called a hitting set of x if a. τ = x * σ is a k-dimensional Γ 0 -flake and there exists an elementary weight function ω τ satisfying the following condition: b. R(τ, ω τ ) < βR p (σ m ) The elementary weight function ω τ is called a hitting map, and we sometimes say σ hits x. A point x = ψ p (y), where y ∈ P (φ, p), is said to be a good point if it is not hit by any set σ ⊂ P with |σ| ≤ m + 1.

A simplex σ which defines a hitting set of x, is necessarily Γ 0 -good. This follows from the requirement that x * σ be a Γ 0 -flake.

The refinement algorithm

In this section, we show that we can refine an -net of M so that the simplices of the Delaunay tangential complex of the refined sample Del T M (P) are power-protected. An -net is a point sample P ⊂ M that is an -sparse -sample set of M for the metric d R N . One can obtain an -net by using a farthest point strategy to select a subset of a sufficiently dense sample set. We will assume that we know the dimension m of the submanifold M and the tangent space T p M at any point p in M.

The algorithm takes as input P 0 , an -net of M, and the positive input parameters , Γ 0 , α < 1 2 , β > 1 and δ 0 < 1 4 . The algorithm refines the input point sample such that: (1) The output sample P ⊇ P 0 is an μ0 -sparse -sample set of M with respect to d R N , where µ 0 = 1 9 .

(2) For all p ∈ P, every m-simplex σ m ∈ star(p; Del T M (P)), σ m is Γ 0 -good and δ 2 0 μ2 0 2 -power protected on T p M.

The algorithm, described in Algorithm 1, applies two rules with a priority order: Rule (2) is applied only if Rule (1) cannot be applied. The algorithm ends when no rule applies any more. Each rule inserts a new point to kill an unfit configuration: either a big configuration or a bad configuration.

A crucial procedure, that selects the location of the point to be inserted, is Pick_valid, given in Algorithm 2. Pick_valid(φ, p) returns a good point ψ p (y) where y ∈ P (φ, p).

The refinement algorithm will also use the procedure Insert(p), given in Algorithm 3.

Analysis of the algorithm

We now turn to the demonstration of the correctness of Algorithm 1. In Section 5.1 we show that the algorithm must terminate, and in Section 5.2 we show that the output of the algorithm

RR n°8273

Algorithm 1 Refinement algorithm Input -net P 0 of M, and input parameters Γ 0 , α and δ 0 ; Initalize P ← P 0 , and calculate Del T M (P); Rule (1) Big configuration ( -big radius): Step 1. Add p to P;

if ∃ p ∈ P such that ∃ σ m ∈ star(p) with R p (σ m ) ≥ , then Insert(ψ p (c p (σ m ))); Rule (2) Bad configuration (Γ 0 -bad): if ∃ p ∈ P and ∃ σ m ∈ star(p) s.t. σ m is Γ 0 -bad, then Insert(Pick_valid(σ m , p)); if ∃ p ∈ P and ∃ σ m+1 ∈ cosph δ0 (p) s.t. σ m+1 is Γ 0 -bad
Step 2. Compute star(p) and cosph δ0 (p);

Step 3. For all x ∈ P \ {p}, update star(x) and cosph δ0 (x); Inria meets the requirements of Theorem 3.5. In order to complete the demonstrations we impose a number of requirements on the input parameters, listed as Hypotheses H0 to H5 below.

Recall that our input parameters are the following positive numbers: , which is the sampling radius and sparsity bound satisfied by P 0 , the input -net sample set; δ 0 , which is used to describe the amount of power-protection enjoyed by the m-simplices in the final complex; Γ 0 , which is used to quantify the quality of the output simplices; β, which is used to describe an upper bound on the radius of the bad configurations that we will avoid; and α, which governs the relative size of the picking region.

It is often convenient to represent the sampling radius by a dimension-free parameter that has the reach of the manifold factored out. We define

˜ = rch(M)
.

The volume of the m-dimensional Euclidean unit-ball is denoted V m . In order to state the hypotheses on the input parameters, we use some additional symbols:

˜ 0 = 1 2 4 (2 4 + 1) 2 , B = 4 + 2(1 + 2 7 3 2 β 2 ) 2 , β = β 1 -2 4 ˜ 0 ,
as well as ξ, E, and D. The term ξ is introduced in Lemma 5.5, and depends on m and rch(M), and the term E, defined in Equation ( 17), depends on ξ and β. The symbol D is introduced in Lemma 5.8, where it is said to depend on m and β.

In order to guarantee termination, we demand the following hypotheses on the input parameters:

H0. α < 1/2 H1. β ≥ 2 (1-δ 2 0 )(1-α-4.5 ˜ 0) H2. Γ 0 < min Vmα m E m+1 β m D , 1 B+1 H3. δ 2 0 ≤ Γ m+1 0 H4. ˜ ≤ min ξ 2(β+β )rch(M) , Γ m+1 0 8β
To meet the quality requirements of Theorem 3.5 we demand an additional constraint on the sampling radius:

H5. ˜ ≤ δ 2 0 Γ 2m 0 1.1×10 9
The make use of the following observation:

Lemma 5.1 From hypotheses H0 to H4 we have ˜ < ˜ 0 and δ 2 0 < 2 4 ˜ 0 , and

(1 -δ 2 0 )(1 -α -4.5˜ 0 ) 4 > 9 def = μ0 . ( 14
)
Boissonnat, Dyer, & Ghosh Proof From H1 we have β > 2 and using the fact that B > β 4 and H2 we have Γ 0 < 1 2 4 +1 . And using the fact, from H4, that

˜ ≤ Γ m+1 0 8β ≤ Γ 2 0 8β < 1 2 4 (2 4 + 1) 2 = ˜ 0 .
Similarly the bound on δ 2 0 follows from H3. Inequality (14) follows from H0 and the definition of ˜ 0 .

From Equation ( 14) we can see that we require β ≥ 4.5. Given α satisfying H0, and a valid choice for β, the hypotheses H2 to H4 sequentially yield upper bounds on the parameters Γ 0 , δ 0 , and ˜ ; we are able to choose parameters that satisfy all of the hypotheses.

The main result of this section can now be summarised:

Theorem 5.2 (Algorithm guarantee) If the input parameters satisfy hypotheses H0 to H5, then Algorithm 1 terminates after producing an intrinsic Delaunay complex Del M (P) that triangulates M.

Termination of the algorithm

This subsection is devoted to the proof of the following theorem:

Theorem 5.3 (Algorithm termination) Under hypotheses H0 to H4, the application of Rule (1) or Rule (2) on a big or a bad configuration φ always leaves the interpoint distance greater than μ0 = 9 , and if φ is a bad configuration then there exists x ∈ P (φ, p) such that ψ p (x) is a good point. Since M is a compact manifold this implies that the refinement algorithm terminates and returns a point sample P which is an μ0 -sparse -sample of the manifold M.

We will prove that at every step the algorithm maintains the following two invariants:

Sparsity: Whenever a refinement rule inserts a new point x = ψ p (y), the distance between x and the existing point set P is greater than μ0 .

Good points: For a bad configuration φ refined by Rule (2), there exists a set of positive volume

G ⊆ P (φ, p) such that if x ∈ G, then ψ p (x) is a good point.
The Termination Theorem 5.3 is a direct consequence of these two algorithmic invariants. We first prove the sparsity invariant in Section 5.1.1, using an induction argument that relies on the fact that the algorithm only inserts good points. The existence of good points is then established in Section 5.1.2, using the sparsity invariant and a volumetric argument. Termination must follow since M is compact and therefore can only support a finite number of sample points satisfying a minimum interpoint distance.

The sparsity invariant

The proof of the sparsity invariant employs the following observation, which serves to bound the distance between a point inserted by Rule (2) and the existing point set:

Inria

Lemma 5.4 Assume Hypotheses H0 to H4. Let φ = σ m ∈ star(p) or φ = p m+1 * σ m ∈ cosph δ0 (p) be a bad configuration being refined by Rule (2). Then for all x ∈ P (φ, p) we have

d R N (c p (σ m ), ψ p (x)) < (α + 4.5˜ 0 )R p (σ m ) and d R N (ψ p (x), P) > (1 -α -4.5˜ 0 )R p (σ m ) > R p (σ m ) 3 .

Proof

Using the facts that α < 1 2 , and ˜ < ˜ 0 , and R p (σ m ) < , we have that for all x ∈ P (φ, p)

p -x < (1 + α)R p (σ m ) < 3 2 < 3 ˜ 0 2 < 1 4 ,
and so we may apply Lemma B.2 to get

x -ψ p (x) ≤ 2 p -x 2 rch(M) ≤ 2(1 + α) 2 R p (σ m ) 2 rch(M) ≤ 4.5 ˜ 0 R p (σ m ), and
c p (σ m ) -ψ p (x) ≤ c p (σ m ) -x + x -ψ p (x) ≤ (α + 4.5 ˜ 0 ) R p (σ m ). (15)
Let S p = ∂B R N (c p (σ m ); R p (σ m )). From Eq. ( 15) we have for

x ∈ P (φ, p) d R N (ψ p (x); P) ≥ d R N (ψ p (x), S p ) > (1 -α -4.5 ˜ 0 ) R p (σ m ) > R p (σ m ) 3 ,
where the final inequality follows from H0 and the definition of ˜ 0 .

We introduce some additional terminology to facilitate the demonstration of the sparsity invariant. An abstract simplex in the initial sample set σ ⊂ P 0 is called an original simplex , otherwise σ ⊂ P is called a created simplex .

Let φ be an unfit configuration that was refined by inserting a point x. We say that x created σ if x ∈ σ and x is the last inserted vertex of the simplex σ, i.e., σ \ {x} already existed just before the refinement of the unfit configuration φ. The unfit configuration φ is called the parent of σ and will be denoted p(σ).

Let σ denote the simplex being refined by the refinement algorithm. We will denote by e(σ) the distance between the point newly inserted to refine σ and the current sample set.

The sparsity invariant is demonstrated by induction. We use a case analysis according to the type of unfit configuration being refined; it is necessary to consider sub-cases. The induction hypothesis is employed only in the sub-case Case 2(b)(ii) and the implicit similar Case 3(b)(ii). The base for the induction hypothesis, i.e., the insertion of the first point, cannot involve Case 2(b) or Case 3(b).

Case 1. Let φ = σ m ∈ star(p) be a big configuration being refined by Rule (1). Since P 0 (⊆ P) is an -net, we have from the fact that ˜ ≤ ˜ 0 < 1 16 and Lemma 4.2, R p (σ m ) ≤ 4 . Rule (1) will refine σ by inserting ψ p (c p (σ m )). Using the fact that ˜ < ˜ 0 < 1 16 , R p (σ m ) ≤ 4 and R p (σ m ) ≥ (since σ m is being refined by Rule (1)), and Lemma B.2, the distance between ψ p (c p (σ m )) and any vertex inserted before ψ p (c p (σ m )) is not less than

R p (σ) -c p (σ) -ψ p (c p (σ)) ≥ R p (σ) - 2R p (σ) 2 rch(M) > (1 -8 ˜ 0 ) > 2 ,
which establishes the sparsity invariant for this case.

Case 2. Consider now the case where φ = p m+1 * σ m ∈ cosph δ0 (p), with σ m ∈ star(p), is being refined by Rule (2). In this case, recalling Lemma 2.5, we have

• R p (σ m ) < , and

• there exists a face of φ that is a Γ 0 -flake.

Let σ 1 ⊆ φ denote a face of φ that is a Γ 0 -flake. We have to now consider two cases:

(a) σ 1 is an original simplex (b) σ 1 is a created simplex
Case 2(a). If σ 1 is an original simplex then σ 1 ⊆ P 0 , and since P 0 is an -net, L(σ 1 ) ≥ . Since a flake must have at least three vertices, σ 1 and σ m must share at least two vertices, and therefore R(σ m ) ≥ /2.

Let x = ψ p (x ) be point inserted to refine φ where x ∈ P (φ, p). Using Lemma 5.4 and the fact that R(σ m ) ≥ /2, we therefore have

d R N (x, P) > (1 -α -4.5 ˜ 0 )R p (σ m ) ≥ (1 -α -4.5 ˜ 0 )R(σ m ) ≥ (1 -α -4.5 ˜ 0 ) 2 > μ0 .
where the final inequality follows from Inequality (14). Hence the sparsity invariant is maintained on the refinement of φ if σ 1 is an original simplex.

Case 2(b) We will now consider the case when σ 1 is a created simplex. We denote by p(σ 1 ) the parent simplex whose refinement gave birth to σ 1 .

We will bound the distance between x = ψ p (x ), where x ∈ P (φ, p), and the point set P. Let x * denote the point whose insertion killed p(σ 1 ). By definition x * is a vertex of σ 1 , and hence also of φ since σ 1 ≤ φ. We distinguish the following two cases:

Case 2(b)(i) Suppose p(σ 1 ) was a big configuration refined by the application of Rule (1). According to Case 1, the lengths of the edges incident to x * in σ 1 are greater than /2.

Inria Therefore d R N (x, P) ≥ (1 -α -4.5 ˜ 0 )R p (σ m ) by Lemma 5.4 ≥ (1 -δ 2 0 )(1 -α -4.5 ˜ 0 )∆(φ) 2 by Lemma 4.4 ≥ (1 -δ 2 0 )(1 -α -4.5 ˜ 0 )∆(σ 1 ) 2 since σ 1 ≤ φ > (1 -δ 2 0 )(1 -α -4.5 ˜ 0 ) 4 by Case 1 > μ0 Inequality 14,
and the sparsity invariant is maintained.

Case 2(b)(ii) Suppose p(σ 1 ) was a bad configuration refined by Rule (2). Thus p(σ 1 ) was either an m-simplex σ m 2 ∈ star(q) or an (m + 1)-simplex q m+1 * σ m 2 ∈ cosph δ0 (q) with σ m 2 ∈ star(q). Consider the elementary weight function ω σ1 = ω φ | σ1 , where ω φ is the weight function (13) identifying φ as a member of cosph δ0 (p). From Lemma 4.1(1), and Lemma 4.4 we have that R p (σ m ) ≥ R(σ 1 , ω σ1 ). We also have that R(σ 1 , ω σ1 ) ≥ β R q (σ m 2 ). Indeed, otherwise σ 1 \ {x * } would be a hitting set for x * , contradicting the hypothesis that p(σ 1 ) was refined according to Rule (2) by the insertion of a good point x * . Thus we have

d R N (x, P) > (1 -α -4.5˜ 0 )R p (σ m ) ≥ (1 -α -4.5˜ 0 )R(σ 1 , ω σ1 ) ≥ (1 -α -4.5˜ 0 ) β R q (σ m 2 ) ≥ (1 -δ 2 0 )(1 -α -4.5˜ 0 )β∆(σ m 2 ) 2 from Lemma 4.4 > ∆(σ m 2 ) from Hypotheses H1 on β > μ0 ,
where the last inequality follows from the induction hypothesis. Again the sparsity invariant is maintained after refinement of φ.

Case 3 The proof for the case of a bad configuration φ = σ m ∈ star(p) to be refined by Rule (2) is similar to Case 2, and the lower bound on the interpoint distances is the same.

This completes the demonstration of the sparsity invariant.

The good points invariant

We will now show that the good point invariant is maintained if φ is a bad configuration being refined by Rule (2). Without loss of generality, we will assume that φ is either equal to σ m ∈ star(p) or to q * σ m ∈ cosph δ0 (p), with σ m ∈ star(p).

Recall the picking region P (φ, p) introduced in Definition 4.6. We will show that there exists y ∈ P (φ, p) such that x = ψ p (y) is a good point. Let Y ⊆ P (φ, p) be the set of points that ψ p maps to a point with a hitting set:

Y = {y ∈ P (φ, p) | ψ p (y) is not a good point}.
We will show that the volume of P (φ, p) exceeds the volume of Y . To this end, we will first bound the number of simplices that could hit some point in ψ p (Y ). Then we will bound the volume that each potential hitting set can contribute to Y .

In order to bound the number of hitting sets, we will use the sparsity invariant together with the following lemma [BG10, Lemma 4.7] to bound the number of points that can be a vertex of a hitting set: Lemma 5.5 (Bound on sparse points) For a point p ∈ M and R > 0, let V be a maximal set of points in B R N | M (p, R) such that the smallest interpoint distance is not less than 2r. There exists ξ that depends on m and rch(M), and A that depends on m, such that if R + r ≤ ξ, then

|V | ≤ 1 + Aξ 1 -Aξ R r + 1 m .
We obtain the following bound on the number of hitting sets:

Lemma 5.6 Let S(φ) denote the set of simplices contained in B + ∩ P that can hit a point in ψ p (Y ). Then

|S(φ)| ≤ E m+1 2 m , (16)
where

E def = 2 1 + Aξ 1 -Aξ (18(α + 2β + 6.5˜ 0 ) + 1) m . (17)
Proof Suppose σ ⊆ P is a hitting set of a point x = ψ p (y), where y ∈ P (φ, p), with |σ| = k and k ≤ m + 1. Let σ = x * σ, and let ω σ denote the corresponding hitting map (see Definition 4.7). Therefore, we have R(σ, ω σ ) < βR p (σ m ), and it follows from Lemma 4.1(2) that

∆(σ) ≤ 2β 1-δ 2 0 R p (σ m
). Thus from Lemma 5.4 and the Triangle inequality we have σ ⊂ B -def = B R N (c p (σ m ), r -), where

r -= 4.5˜ 0 R p (σ m ) + αR p (σ m ) + 2β 1 -δ 2 0 R p (σ m ).
Let c = ψ p (c p (σ m )). Then using Lemma B.2 and the fact that R p (σ m ) < < ˜ 0 rch(M) we have

c p (σ m ) -c ≤ 2R p (σ m ) 2 rch(M) < 2˜ R p (σ m ) ≤ 2˜ 0 R p (σ m ).
Using δ 2 0 < 2 4 ˜ 0 from Lemma 5.1, and β = β 1-2 4 ˜ 0 , and R p (σ m ) < , we find

c p (σ m ) -c + r -≤ αR p (σ m ) + 6.5˜ 0 R p (σ m ) + 2β 1 -δ 2 0 R p (σ m ) ≤ (α + 2β + 6.5 ˜ 0 ) def = R. Thus B -⊆ B + def = B R N (c, R)
, and y ∈ Y if and only if there exists σ ⊂ B + ∩ P such that σ hits ψ p (y).

Using Lemma 5.5 we will bound the number of sample points in B + ∩ P. Set r = μ0 2 = 18 and observe that

R + r = 1 18 + α + 2β + 6.5 ˜ 0 ≤ (2β + 1) ≤ ξ,
Inria by Hypothesis H4. The sparsity invariant and Lemma 5.5 then yields

|B + ∩ P| ≤ 1 + Aξ 1 -Aξ × (α + 2β + 6.5˜ 0 ) 1/18 + 1 m = E 2 .
Since the number of k-simplices is less than E 2 k+1 , and the maximum dimension of a hitting set is m, we have

|S(φ)| ≤ E m+1 2 m .
We now turn to the problem of bounding the volume of Y . We will consider the contribution of each σ ∈ S(φ). The following definition characterises the set of points in M that can be hit by σ: Definition 5.7 (Forbidden region) For a k-simplex σ with vertices in M with k ≤ m and parameter t < , the forbidden region, F (σ, t), is the set of points x ∈ M such that σ 1 = x * σ satisfies the following conditions:

• L(σ 1 ) > t 9 • σ 1 is a Γ 0 -flake • there exists an elementary weight function ω σ1 s.t. R(σ 1 , ω σ1 ) < βt
We will use the following lemma, which is proved in Appendix C. It bounds the volume of the set of points that can be hit by a given simplex: Lemma 5.8 (Volume of forbidden region) Let σ be a k-simplex with vertices on M and

k ≤ m. If 1. Γ 0 ≤ 1 B+1 , 2. ˜ ≤ min{ ξ 4β rch(M) , Γ m+1 0 8β } and 3. δ 2 0 ≤ min{Γ m+1 0 , 1 4 }, then vol(F (σ, t)) ≤ D Γ 0 R(σ) m ,
where D depends on m and β.

Lemma 5.8, together with Lemma 5.6, yields a bound on the set of points Y in the picking region that do not map to a good point:

Lemma 5.9 The volume of the set Y ⊂ P (φ, p) of points that do not map to a good point is bounded as follows:

vol(Y ) ≤ E m+1 β m D Γ 0 R p (σ m ) m . Proof Let t 0 = R p (σ m ) < .
For a given σ ∈ S(φ), let Y σ ⊆ Y be the set of points y for which σ hits x = ψ p (y). Then from Hypotheses H0 to H4 and Lemma 5.8, we have

vol(Y σ ) ≤ vol(π p (F (σ, t 0 ))) ≤ vol(F (σ, t 0 )) since π p is a projection map on T p M ≤ D Γ 0 R(σ) m . (18)
Let σ 1 = x * σ, and let ω σ1 be the corresponding hitting map. From the definition of hitting sets and hitting maps, we have R p (σ m ) < , and R(σ 1 , ω σ1 ) < βR p (σ m ) and σ is Γ k 0 -thick. Define ω σ = ω σ1 | σ . Then, using Lemma 4.1 (3) and the fact that R(σ, ω σ ) ≤ R(σ 1 , ω σ1 ) < βR p (σ m ), we have

R(σ) ≤ R(σ, ω σ ) 1 - δ 2 0 Υ(σ) -1 ≤ R(σ, ω σ ) 1 - δ 2 0 Γ m 0 -1 since Υ(σ) ≥ Γ k 0 ≥ Γ m 0 ≤ 2R(σ, ω σ ) since δ 2 0 Γ m 0 ≤ Γ 0 < 1 2 from Hyp. H2, H3 < 2βR p (σ m ). (19)
The inequalities ( 18) and ( 19) together yield

vol(Y σ ) ≤ 2 m β m D Γ 0 R p (σ m ) m ,
and so using Lemma 5.6 we have

vol(Y ) = vol   σ∈S(φ) Y σ   ≤ σ∈S(φ) vol(Y σ ) ≤ E m+1 β m D Γ 0 R p (σ m ) m .
By the definition of the picking region, we have that

vol(P (φ, p)) = V m α m R p (σ m ) m .
By Hypothesis H2, E m+1 β m D Γ 0 R p (σ m ) m is less than vol(P (φ, p)), the volume of the picking region of φ. Thus with Lemma 5.9, this proves the existence of points y in the picking region P (φ, p) of φ such that ψ p (y) is a good point. The proof of Theorem 5.3 is complete.

Output quality

We will now show that if Hypothesis H5 is satisfied, in addition to Hypotheses H0 to H4, then the output to the refinement algorithm will meet the demands imposed by Theorem 3.5, thus yielding Theorem 5.2. The main task is to ensure that every m-simplex in Del T M (P) has, for each vertex, a δ 2 0 μ2 0 2power-protected Delaunay ball centred on the tangent space of that vertex. This is achieved in two steps. First we establish conditions to ensure that cosph δ0 (p) = ∅ for every p ∈ P. As noted by Lemma 4.5, this ensures that every simplex in star(p) has a δ 2 0 μ2 0 2 -power-protected Delaunay ball centred on T p M. Next we show conditions such that if σ m ∈ Del T M (P), then σ m ∈ star(p) for every vertex p ∈ σ m . In each step the required conditions impose an additional constraint on the sampling radius, and this leads to Hypothesis H5.

As a starting point, we observe the following direct consequence of the Termination Theorem 5.3: Inria Corollary 5.10 Under Hypotheses H0 to H4, for all p ∈ P, the output of the algorithm satisfies the following:

1. σ m ∈ star(p) =⇒ R p (σ m ) < and σ m is a Γ 0 -good simplex, and 2. all σ m+1 ∈ cosph δ0 (p) are Γ 0 -good.

We will show that for an appropriate sampling radius, there cannot be a Γ 0 -good simplex in cosph δ0 (p). We exploit the following bound on the thickness of a small (m + 1)-simplex: Lemma 5.11 (Small (m + 1)-simplices are not thick) Let σ m+1 be an (m + 1)-simplex with vertices in M and ∆(σ m+1 ) < rch(M). For distinct vertices p, q ∈ σ m+1 define θ = ∠(aff(σ q ), T p M). Then

Υ(σ m+1 ) ≤ ∆(σ m+1 ) 2 rch(M) + sin θ .

Proof

We will bound the altitude D(q, σ m+1 ). Let be the line through p and q. Using Lemma B.1 and the fact that ∠(aff(σ q ), T p M) = θ, we get

D(q, σ m+1 ) = d R N (q, aff(σ q )) = sin ∠( , aff(σ q )) × d R N (p, q) ≤ (sin ∠( , T p M) + sin ∠(aff(σ q ), T p M)) × d R N (p, q) ≤ d R N (p, q) 2 rch(M) + sin θ × d R N (p, q) ≤ ∆(σ m+1 ) 2 rch(M) + sin θ × ∆(σ m+1 ).
Therefore we have

Υ(σ m+1 ) ≤ ∆(σ m+1 ) 2 rch(M) + sin θ .
Also, Whitney's Lemma 2.1 implies that a Γ 0 -good simplex in star(p) makes a small angle with the tangent space at p:

Lemma 5.12 If σ m ∈ star(p) is Γ 0 -good with R p (σ m ) < , then sin θ < 2 Γ m 0 rch(M)
,

where θ = ∠(aff(σ m ), T p M). Proof Let ζ = max x∈σ m d R N (x, T p M) where x is a vertex of σ m . From Lemma B.1, we have ζ = max x∈σ m d R N (x, T p M) ≤ max x∈σ m d R N (p, x) 2 2rch(M) ≤ ∆(σ m ) 2 2 rch(M) .
Lemma 5.15 Let P be a μ0 -sparse -sample of M with μ0 ≤ 1 independent of . We further assume δ 0 ≤ 1 and

(1) for all p ∈ P, every σ m ∈ star(p) is a Γ 0 -good simplex with R p (σ m ) < , and

(2) for all p ∈ P, cosph δ0 (p) = ∅.

If ˜ ≤ δ 2 0 μ2 0 Γ m 0 36 ,
then star(p) = star(p; Del T M (P)) for all p in P.

Proof

For p ∈ P, let σ m ∈ star(p) and q ( = p) be a vertex of σ m . We will show that σ m is also in star(q).

Let θ = max ∠(aff(σ m ), T x M) where the max is taken over the vertices x of σ m . Since

˜ ≤ δ 2 0 μ2 0 Γ m 0 36 < Γ m 0 4 ,
Lemma 5.12 yields

sin θ ≤ 2˜ Γ m 0 def = c 1 ˜ < 1 2 .
It follows that cos θ > √ 3/2 and so tan θ ≤ 2c 1 ˜ .

Recall that N (σ m ) denotes the affine space orthogonal to aff(σ m ) and passing through C(σ m ). Let c be the unique point in N (σ m ) ∩ T q M, and let R = d R N (c, p).

Using the fact that ∠(aff(σ m ), T q M) ≤ θ, we have

d R N (C(σ m ), c) ≤ R(σ m ) tan θ ≤ 2c 1 ˜ R(σ m ),
and likewise

d R N (C(σ m ), c p (σ m )) ≤ 2c 1 ˜ R(σ m ).
It follows that R ≤ (1 + 2c 1 ˜ )R(σ m ), and d R N (c p (σ m ), c) ≤ 4c 1 ˜ R(σ m ). From the above observations, and using the fact that R(σ m ) ≤ R p (σ m ) < , we get

B R N (c, R) ⊆ B R N (c p (σ m ), (1 + 6c 1 ˜ )R(σ m )) ⊆ B R N (c p (σ m ), R p (σ m ) + 6c 1 ˜ ).
Since cosph δ0 (p) = ∅, and P is μ0 -sparse, we have that σ m is δ 2 0 μ2 0 2 -power protected on T p M (Lemma 4.5). This means that

B R N (c p (σ m ), R p (σ m ) + ∆) ∩ (P \ σ m ) = ∅, where ∆ = R p (σ m ) 2 + δ 2 0 μ2 0 2 -R p (σ m ) = δ 2 0 μ2 0 2 R p (σ m ) 2 + δ 2 0 μ2 0 2 + R p (σ m ) > δ 2 0 μ2 0 1 + δ 2 0 μ2 0 + 1 > δ 2 0 μ2 0 3 def = c 2 .

RR n°8273

Since 6c 1 ˜ ≤ c 2 , by our hypothesis on ˜ , we have

B R N (c, R) ⊂ B R N (c p (σ m ), R p (σ m ) + ∆),
and thus the m-simplex σ m belongs to star(q).

The consequence of Lemma 5.15, together with Lemma 4.5 is that every m-simplex in Del T M (P) has, for each vertex, a δ 2 0 μ2 0 2 -power-protected Delaunay ball centred on the tangent space of that vertex:

Corollary 5.16 Let P be a μ0 -sparse -sample of M with μ0 being independent of . Under the hypotheses in Lemma 5.15, for all p ∈ P, all the m-simplices σ m in star(p; Del T M (P)) are δ 2 0 μ2 0 2 -power protected on T p M. I.e, for all σ m ∈ star(p; Del T M (P)) there exists a c p (σ m ) ∈ N (σ m ) ∩ T p M such that for all q ∈ P \ σ m

d R N (q, c p (σ m )) 2 > d R N (p, c p (σ m )) 2 + δ 2 0 μ2 0 2 .
We are now in a position to show that Hypothesis H5, when added to Hypotheses H0 to H4, results in the output of the algorithm meeting the demands of Theorem 3.5.

Recalling that μ0 = 1 9 , Hypotheses H3 yields the following consequence of H5:

˜ ≤ δ 2 0 Γ 2m 0 1.1 × 10 9 ≤ min Γ 2m+1 0 4 , δ 2 0 μ2 0 Γ m 0 36 , δ 2 0 μ3 0 Γ 2m 0 1.5 × 10 6 .
In other words, the sampling radius bounds demanded by Corollary 5.14, Lemma 5.15, and Theorem 3.5 are all simultaneously satisfied. Corollary 5.10 together with Corollary 5.14 ensure that the hypotheses of Lemma 5.15 are satisfied, and so it follows that the m-simplices of Del T M (P) are power-protected as described by Corollary 5.16. Thus all the requirements of Theorem 3.5 are satisfied, and we obtain Theorem 5.2.

Conclusions

We have described an algorithm which meshes a manifold according to extrinsic sampling conditions which guarantee that the intrinsic Delaunay complex coincides with the restricted Delaunay complex, and that it is homeomorphic to the manifold. The algorithm constructs the tangential Delaunay complex, which is also shown to be equal to the intrinsic Delaunay complex, and in this way we are able to exploit existing structural results [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF] to obtain the homeomorphism guarantee. This approach relies on an embedding of M in R N . In future work we aim to develop algorithms and structural results which enable the construction of an intrinsic Delaunay triangulation in the absence of an embedding in Euclidean space.

A An obstruction to intrinsic Delaunay triangulations

When meshing Riemannian manifolds of dimension 3 and higher using Delaunay techniques, flake simplices pose problems which cannot be escaped simply by increasing the sampling density. In particular, developing an example on a 3-manifold presented by Cheng et al. [START_REF] Cheng | Manifold reconstruction from point samples[END_REF], Boissonnat et al. [START_REF] Boissonnat | Manifold reconstruction in arbitrary dimensions using witness complexes[END_REF]Lemma 3.1] show that the restricted Delaunay triangulation need not be homeomorphic to the original manifold, even with dense well separated sampling.

In this appendix we develop this example from the perspective of the intrinsic metric of the manifold. It can be argued that this is an easier way to visualize the problem, since we confine our viewpoint to a three dimensional space and perturb the metric, without referring to deformations into a fourth ambient dimension. This viewpoint also provides an explicit counterexample to the results announced by Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF]: In general the nerve of the intrinsic Voronoi diagram is not homeomorphic to the manifold. The density of the sample points alone cannot guarantee the existence of a Delaunay triangulation.

We explicitly show how density assumptions based upon the strong convexity radius cannot escape the problem. The configuration considered here may be recognised as essentially the same as that which was described qualitatively in Section 2.4.3, but here we consider the Voronoi diagram rather than Delaunay balls. We work exclusively on a three dimensional domain, and we are not concerned with "boundary conditions"; we are looking at a coordinate patch on a densely sampled compact 3-manifold.

A.1 Sampling density alone is insufficient

We will now construct a more explicit example to demonstrate that the problem of neardegenerate configurations cannot be escaped with the kind of sampling criteria proposed by Leibon and Letscher [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF].

Leibon and Letscher [LL00, p. 343] explicitly assume that the points are generic which they state as Definition A.1 The set P ⊂ M, is generic if M is an m-manifold and m + 2 points never lie on the boundary of a round ball.

Here a round ball refers to a geodesic ball. This definition of genericity is natural, and corresponds to Delaunay's original definition [START_REF] Delaunay | Sur la sphère vide[END_REF], except Delaunay only imposed the constraint on empty balls. A question that Delaunay addressed explicitly, but which was not addressed by Leibon and Letscher, is whether or not such an assumption is a reasonable one to make. Delaunay showed that any (finite or periodic) point set in Euclidean space can be made generic through an arbitrarily small affine perturbation. That a similar construction of a perturbation can be made for points on a compact Riemannian manifold has not been explicitly demonstrated. However, in light of the construction we now present, it seems that the question is moot when m > 2, because an arbitrarily small perturbation from degeneracy will not be sufficient to ensure a triangulation.

Leibon and Letscher proposed adaptive density requirements based upon the strong convexity radius. These requirements are somewhat complicated, but they will be satisfied if a simple constant sampling density requirement is satisfied. Exploiting a theorem [Cha06, Thm. IX.6.1], that relates the strong convexity radius to the injectivity radius, inj(M), and a positive bound on the sectional curvatures, they arrive at the following: Claim A.2 ([LL00, Lemma 3.3]) Suppose K 0 is a positive upper bound on the sectional curvatures of M, and

η(M) = min inj(M) 10 , π 10 √ K 0 . ( 20
)
If P is an η(M)-sample set for M with respect to d M , then |Del M (P)| ∼ = M.

In fact, we will show that no sampling conditions based on density alone will be sufficient to guarantee a homeomorphic Delaunay complex in general, even when a sparsity assumption is also demanded. An ˜ -net is an ˜ -sparse, ˜ -sample set. We will show:

Theorem A.3 With η(M) as defined in Equation (20), for any > 0, there exists a compact Riemannian manifold M, and and a finite set P ⊂ M, such that P is an ( η(M))-net for M, with respect to the metric d M , but Del M (P) is not homeomorphic to M.

A.1.1 A counter-example

We will construct the counter-example by considering a perturbation of a Euclidean metric. This is a local operation, and the global properties of the manifold are only relevant in so far as they affect η(M) of Equation (20). We may assume, for example, that the manifold is a 3-dimensional torus M ∼ = S 1 × S 1 × S 1 , initially with a flat metric. Thus assume there is some 0 such that any compact Riemannian manifold may be triangulated by the intrinsic Delaunay complex when P is an 0 η(M)-net. For convenience, we choose a system of units so that 0 η(M) = 1. We will first construct a point configuration and metric perturbation that leads to a problem, and then we will show that the sampling assumptions are indeed met. We introduce a number of parameters which we will manipulate to produce the counter-example. We are exploiting the fact that the genericity assumption allows configurations that are arbitrarily close to being degenerate. The assumed 0 has been fixed.

We will work within a coordinate chart on M, where the metric is Euclidean. We will perturb this metric by constructing a metric tensor g, and we will denote by M the manifold with with this new metric.

Consider points u, v, w, p in the xz-plane arranged with u and v at ±a on the z axis, and w and p at ±(a + ξ) on the x axis, with a = 3 4 , and 0 < ξ < r 0 γ, where r 0 and γ will be specified below. The Voronoi diagram of these points in the xz-plane is shown in Figure 2. The main point here is that the Voronoi boundary between V M (u) and V M (v) may be arbitrarily small with respect to the distance between the sites, i.e., ξ will be very very small.

The three dimensional Voronoi diagram is the extension of this in the horizontal y-direction, so that every crosssection looks the same. Note that since the points are not co-circular, they do not represent a degeneracy by Delaunay's criteria [START_REF] Delaunay | Sur la sphère vide[END_REF], but this is irrelevant; we will also argue that the points will not represent a degenerate configuration with respect to the new metric.

We now introduce a small localized metric perturbation so as to change the Voronoi diagram near the origin. For example, we can demand that the matrix of the metric tensor in our coordinate system has the form

g(p) =   1 -f (|p|) 0 0 0 1 0 0 0 1   ,

Inria

where |p| is the parametric distance from p to the origin. The radial function f is non-negative, and it and its first two derivatives are bounded, e.g.,

f (r), |f (r)| , |f (r)| ≤ β. (21)
We also demand that there exists a positive γ ≤ β such that f (r) ≥ γ when r ≤ r 0 , and that f (r) = 0 if r ≥ 2r 0 . The parameter r 0 , defines the radius of the ball bounding the perturbed region. Now we have d(w, p) < d(u, v) when ξ < r 0 γ. Since γ may be arbitrarily small compared to β, standard arguments supply a function f meeting these conditions. For example, the C ∞ construction described by Munkres [Mun68,p. 6] may be multiplied by a scalar sufficiently small to meet our needs.

The vertical y = 0 cross-section of the perturbed Voronoi diagram will look something like Figure 3: V M(p) and V M(w) now meet in the xz-plane, and V M(u) and V M(v) do not. However, since geodesics which do not intersect the ball B R 3 (0, 2r 0 ) will remain straight lines in the parameter space, the Voronoi diagram is unchanged outside of a neighbourhood of the origin. Thus looking from above at the slice of the Voronoi diagram in the xy-plane, we will see something like Figure 4(a). Two Voronoi vertices have been introduced, the red and blue points in Figure 4. These are the centres of distinct empty geodesic circumballs for {p, u, v, w}. Since they cannot lie in the region unaffected by the perturbation, a quick calculation shows that the parametric distance of these Voronoi vertices from the origin is bounded by 4r 0 , when r 0 ≤ 1 4 , and it follows from another small calculation that the parametric distance from these Voronoi vertices to any of the four sample points is bounded by a(1+

3ξ+16r 2 0 a 2
). The distances between these Voronoi vertices and the sample points in the new metric will also be subjected to the same bound, since no distances increase. Also, The sparsity condition will not be affected by the perturbation. Thus, since we can make r 0 as small as we please, and ξ is chosen such that ξ < r 0 γ, it follows that the radius of these balls may be made arbitrarily close to a = 3 4 = 3 4 0 η(M). We will argue next that we can make η(M)η( M) as small as desired by reducing the size of β in Equation (21). Then other sample points may be placed on the manifold so that the density criteria are met, and no degenerate configuration (violation of Definition A.1) need be introduced. This means that the Delaunay complex, defined as the nerve of the Voronoi diagram, will not be a triangulation of the manifold M. As observed by Boissonnat et al. [START_REF] Boissonnat | Manifold reconstruction in arbitrary dimensions using witness complexes[END_REF], the triangle faces {p, w, u} and {p, w, v} will be adjacent to only a single tetrahedron, namely {p, u, v, w}. Thus Del M(P ) is not a manifold complex as defined in Section 2. This is clearly a problem if the original manifold has no boundary.

Although it is in some sense close to being degenerate, we emphasise that this configuration represents a problem that cannot be escaped by an arbitrarily small perturbation of the sample points. An argument based on the triangle inequality shows that in order to effect a change in the topology of the Voronoi diagram, a displacement of the points by a distance of Ω(r 0 γξ) is required.

More specifically, we observe that the configuration {p, u, v, w} may be placed in an otherwise well behaved point set P such that within a small ball centred at the origin in our coordinate chart, all points will have {p, u, v, w} as the four closest points in P, and this would remain the case even if the point positions were perturbed a small amount. We may further assume that the other Delaunay simplices are well shaped, so that stability results [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] can be used to argue that they cannot be destroyed with an arbitrarily small perturbation. Then we argue that in order to obtain a triangulation by a perturbation P → P , we must ensure that the Voronoi cell V M({p , w }) must vanish: the edge {p , w } will never be incident to any tetrahedron other than {p , u , v , w }. Then an argument based on the triangle inequality shows that for a ρ-perturbation with ρ < r0γ-ξ 6 , there will be a point in V M({p , w }) within a distance of 2ρ of the origin.

A.1.2 The sizing function under perturbation

We need to establish that the metric manipulation that we performed in order to construct the counter-example, does not have a dramatic effect on the sizing function η( M). This follows from the fact that we have bounded gg together with its first and second deriviatives. Since the sectional curvature may be described as a continuous function of g and it's first and second derivatives [dC92, pp. 56 & 93], the effect of our perturbation on the sectional curvatures can be made arbitrarily small by reducing β in Equation ( 21).

Since we started with a flat metric anyway, the bound K 0 can be made arbitrarily small, and so the second term in Equation (20) will not be the smallest. We need to bound the change in the injectivity radius as well.

This follows from results in the literature [START_REF] Ehrlich | Continuity properties of the injectivity radius function[END_REF][START_REF] Sakai | On continuity of injectivity radius function[END_REF], which state that for a compact manifold, inj(M) depends continuously on the metric and its first and second derivatives. Specifically, Inria Lemma A.4 (Ehrlich) Let M be the space of C 3 Riemannian metric structures g on a compact manifold M, and endow M with the C 2 topology. The function g → inj g (M) is continuous in this topology. This means that for any desired bound on η(M)η( M) , there will be a β that will satisfy the bound.

The construction of the counter-example is complete.

A.2 Discussion

We have shown that for constructing a Delaunay triangulation for an arbitrary Riemannian manifold, a sampling density requirement is not sufficient in general. The solution we propose in the body of this paper, is to constrain the kind of sample sets that we consider. Another approach would be to constrain the kind of metrics that are assumed. However, even with a purely Euclidean metric, allowing configurations to be arbitrarily close to degeneracy means that arbitrarily poorly shaped simplices are to be expected. When the metric is no longer Euclidean, the "shape" of a simplex no longer has an obvious meaning, but the problems associated with point configurations near degeneracy will certainly be present. Our analysis relied on the ability to make the support of the perturbation small. This is unlikely to be a necessary feature of the construction, but it facilitates our simplistic analysis.

Clarkson [START_REF] Clarkson | Building triangulations using epsilon-nets[END_REF] remarked that an implication of Leibon and Letscher's claim [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF] is that for four points close enough together, there is a unique circumsphere with small radius. Our counter-example shows that circumcentres need not be unique under these conditions. In fact the existence of unique circumcentres does not follow from the triangulation result: In our work we do not claim that the m-simplices have a unique circumcentre in the intrinsic metric. However, the argument sketched out by Leibon and Letscher claimed that the intrinsic Voronoi diagram is a cell complex (i.e., it satisfies the closed ball property [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]), and this does imply unique circumcentres for the top dimensional simplices.

It is worth emphasising that the problems discussed here only arise when the dimension is greater than 2. The same sampling criteria for two dimensional manifolds has been fully validated [START_REF] Leibon | Random Delaunay triangulations, the Thurston-Andreev theorem, and metric uniformization[END_REF][START_REF] Dyer | Surface sampling and the intrinsic Voronoi diagram[END_REF], however these works both assume genericity in the sample set, without demonstrating that it is a reasonable assumption.

B Background results for manifolds

The tangent space at p ∈ M is denoted T p M, and we identify it with an m-flat in the ambient space. The normal space, N p M, is the orthogonal complement of T p M in T p R N , and we likewise treat it as the affine subspace of dimension mk orthogonal to T p M ⊂ R N .

A ball B = B R N (c, r) is a medial ball at p if B ∩ M = ∅, it is tangent to M at p, and it is maximal in the sense that any ball which contains B either coincides with B or intersects M. The local reach at p is the infimum of the radii of the medial balls at p, and the reach of M, denoted rch(M), is the infimum of the local reach over all points of M. In order to approximate the geometry and topology with a simplical complex, manifolds with small reach require a higher sampling density than those with a larger reach. As is typical, an upper bound on our sampling radius will be proportional to rch(M). Since M ⊂ R N is a smooth, compact embedded submanifold, it has positive reach.

An estimate of how the tangent space locally deviates from the manifold is given by an observation of Federer [Fed59, Theorem 4.8(7)] (see also Giesen and Wagner [GW04, Lemma 6]):

Using the fact that pp k+1 ≤ δ 0 L(σ 1 ), we get pp ≤ D(p k+1 , σ 1 ) + δ 0 L(σ 1 )

We will now consider the triangle σ 2 = [p q r]. Note that C(σ 1 , ω σ1 ), R(σ 1 , ω σ1 ) are the circumcenter and radius of σ 2 respectively. Also, C(σ, ω σ ) is the midpoint of the line segment [q r] with 2R(σ, ω σ ) = qr and D(p, σ 2 ) = pp . From the definition of q, we have pr ≥ pq . Using the fact qr = 2R(σ, ω σ ), we have pr ≥ qr 2 = R(σ, ω σ ).

This implies from Lemma C.1

pq = 2R(σ 2 )D(p, σ 2 ) pr ≤ 2R(σ 1 , ω σ1 )D(p, σ 2 ) R(σ, ω σ ) as R(σ 2 ) ≤ R(σ 1 , ω 1 ) and pr ≥ R(σ, ω σ ) ≤ 48βD(p, σ 2 ) = 48β pp as Eq. (24) (28)

From Eq. ( 26), ( 27) and (28) p k+1q ≤ p k+1p + pq ≤ δ 0 L(σ 1 ) + 48β(D(p k+1 , σ 1 ) + δ 0 L(σ 1 ))

def = η 1 (29)
Using the fact that σ is Γ k 0 -thick (since σ 1 is a Γ 0 -flake), and the bound δ 2 0 L(σ 1 ) 2 on the differences of the squared distances between C(σ, ω σ ) and the vertices of σ, we obtain a bound [BDG12, Lemma 4.1] on the distance from C(σ, ω σ ) to C(σ): C(σ) -C(σ, ω σ ) ≤ δ 2 0 L(σ 1 ) 2 2Υ(σ)∆(σ)

≤ δ 2 0 R(σ) Υ(σ) as L(σ 1 ) ≤ ∆(σ) ≤ 2R(σ) ≤ δ 2 0 R(σ) Γ k 0 since σ is Γ k 0 -thick, Υ(σ) ≥ Γ k 0 ≤ δ 2 0 R(σ) Γ m 0 as Γ 0 ≤ 1 def = η 2 (30)
Since k ≥ 1, there exists p i ∈ σ such that

p i ∈ B R N (C(σ), R(σ)) ∩ B R N (C(σ, ω σ ), R(σ, ω σ )) ∩ aff(σ).

Inria

Figure 1 :

 1 Figure 1: In three dimensions, three closed geodesic balls can all touch three points, u, v, p, on their boundary and yet no one of them is contained in the union of the other two.

  is the centre of a Delaunay ball for σ. Thus every Voronoi face corresponds to a Delaunay simplex. The Voronoi cells give a decomposition of U , denoted Vor d (P), called the Voronoi diagram. Our definition of the Delaunay complex of P corresponds to the nerve of the Voronoi diagram.

  Proposition 3.9 (Equating Delaunay complexes) Suppose P ⊂ M is an -sample set with respect to d R N | M , and that for every p ∈ P, in the local Euclidean metric on W = B R N | M (p, 5 ), every m-simplex in star(p; Del(P W )) is secure, where P W = P ∩ W , and δ = ν 0 . If ≤ Υ 0 µ 0 ν 0 rch(M) 20700 then star(p; Del(P W )) = star(p; Del R N | M (P W )) = star(p; Del M (P W )). (8) Thus Del R N | M (P) = Del M (P), and they are manifold complexes.

  then σ = π p (σ) has a δ-protected Delaunay ball B R m (c, r) with respect to the local Euclidean metric d R m for p on any admissible neighbourhood W that contains B R N | M (p, 3 ), and δ = ν 0 , with

  , then Insert(Pick_valid(σ m+1 , p)); Output Del T M (P) = ∪ p∈P star(p); Algorithm 2 Pick_valid(σ, p) // Assume that σ is either equal to σ m ∈ star(p) // or σ m+1 = p m+1 * σ m ∈ cosph δ0 (p) with σ m ∈ star(p) Step 1. Pick randomly y ∈ P (σ m , p) (or P (σ m+1 , p)); // Recall that ψ p projects points from T p M onto M along N p M Step 2. x ← ψ p (y); Step 3. Avoid hitting sets: // |σ| denotes the cardinality of σ if ∃ σ ⊂ P, with |σ| ≤ m + 1, which is a hitting set of x, then discard x, and go back to Step 1; Step 4. Return x; Algorithm 3 Insert(p)

Figure 2 :

 2 Figure 2: A vertical slice: the xzplane of the initial Voronoi diagram, seen from the negative y axis.

  Figure 4(b) shows the yz-plane.

Figure 3 :

 3 Figure 3: The y = 0 slice of the perturbed Voronoi diagram.

Figure 4 :

 4 Figure 4: Looking at cross-sections; the positive y-direction is to the right. The four points, p, u, v, w, admit two small circumballs with distinct centres (the red and blue points).

  Also, C(σ)p i = R(σ) and C(σ, ω σ )p i = R(σ, ω σ ).Using the facts that R(σ) = C(σ)p i and R(σ, ω σ ) = C(σ, ω σ )p i , and the Triangle inequality, we getR(σ) -C(σ) -C(σ, ω σ ) ≤ C(σ, ω σ )p i ≤ R(σ) + C(σ) -C(σ, ω σ ) R(σ)η 2 ≤ R(σ, ω σ ) ≤ R(σ) + η 2(31)

RR n°8273

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

Acknowledgements

This work was partially supported by the CG Learning project. The project CG Learning acknowledges the financial support of the Future and Emerging Technologies (FET) programme within the Seventh Framework Programme for Research of the European Commission, under FET-Open grant number: 255827.

Inria

Boissonnat, Dyer, & Ghosh

Using Lemma 2.1 and the facts that R(σ m ) ≤ R p (σ m ) < and Υ(σ m ) ≥ Γ m 0 (since σ m is a Γ 0 -good simplex), we have

Using Lemmas 5.11 and 5.12 we get that no (m + 1)-dimensional simplices in cosph δ0 (p) can be Γ 0 -good when is sufficiently small:

since R p (σ m ) < and δ 2 0 ≤ 1 2 . Then from Lemmas 5.11 and 5.12 we get

from the hypothesis on ˜ .

We emphasise the consequence of Lemma 5.13:

and all the simplices in cosph δ0 (p) are Γ 0 -good, then

Now we proceed to the second step of the analysis. Assuming that cosph δ0 (p) = ∅ for all p in P, the following lemma says that if σ ∈ star(p), then also σ ∈ star(q) for every vertex q ∈ σ, provided the appropriate constraints are met.

2rch(M) , and thus sin α ≤ r 2rch(M) , where α is the angle between [x, y] and T x M.

A complementary result bounds the distance to the manifold from a point on a tangent space [BG10, Lemma 4.3]:

rch(M) . The previous two lemmas lead to a convenient bound on the angle between nearby tangent spaces. We prove here a variation on previous results [NSW08, Prop. 6.2] [BG11, Lemma 5.5]:

4

, and let α be the angle between T x M and T y M. Then, sin α < 6r rch(M) .

Proof

Let v ∈ T y M ⊂ R N with vy = r. We will bound the angle between vy and T x M. We have

where

rch(M) . For the third term in Equation ( 22), we find

and so we may apply Lemma B.1 to obtain d R N (v, T x M) ≤ 25r 2 8rch(M) . Putting these observations back into Equation ( 22) we find

.

The following observation is a direct consequence of results established by Niyogi et al. [NSW08, Lemma 5.4]:

] that the Jacobian of f is nonsingular on W , so that W is a covering space for U = f (W ) ⊂ T p M. The Morse-theory argument of Boissonnat and Chazals [BC01, Proposition 12 ] can be applied to demonstrate that W is a topological ball. It follows that U is connected, since any path in W projects to a path in U . Thus W must be a single-sheeted cover of U , since f -1 (0) = {p}. Indeed, if q ∈ W with q = p and f (q) = 0, then [p, q] would be perpendicular to T p M, contradicting Lemma B.1. Thus f : W → U is a diffeomorphism.

Inria

Constructing intrinsic DTs

45

Niyogi et al [NSW08, Prop 6.3] demonstrate a bound on the geodesic distance between nearby points, with respect to the ambient distance. We will use a modified statement of this result:

Proof

The announced result states

under the same hypothesis on x and y. Rearranging, we have

where the second inequality is obtained by squaring away the radical.

C Forbidden volume calculation

In this appendix we demonstrate: Lemma 5.8(Volume of forbidden region) Let σ be a k-simplex with vertices on M and

where D depends on m and β.

We will use the following lemmas in the proof of Lemma 5.8:

Lemma C.1 (Triangle altitude bound) For any non-degenerate triangle σ = [p, q, r], we have

Proof Let α = ∠prq and observe that sin α = pq 2R(σ) .

Since D(p, σ) = pr sin α, the result follows.

RR n°8273

Boissonnat, Dyer, & Ghosh

and the following conditions are satisfied:

(1) L(σ 1 ) > t 9 , (2) R(σ 1 , ω σ1 ) < βt,

(3) σ 1 is a Γ 0 -flake, and

where S = B R N (C(σ), R(σ)) ∩ aff(σ) and

Proof

Let ω σ = ω σ1 | σ . Note that ω σ : σ → [0, ∞) is an elementary weight function, and C(σ, ω σ ) is the orthogonal projection of C(σ 1 , ω σ1 ) onto aff(σ).

From Lemma 4.1 (2) and the fact that δ 2 0 ≤ 1 4 , we have

and

Therefore, from Lemma 2.7, we have

Let p be the point closest to

Let q be the point closest to p on ∂B R N (C; R) ∩ aff(σ), p be the projection of p onto aff(σ), and and let r denotes the intersection of the line aff([q C(σ, ω σ )]) with ∂B R N (C; R). Note that C(σ 1 , ω σ1 ), C(σ, ω σ ), p k+1 , p, p , q and r lie on the same 2-dimensional affine space.

Inria

Boissonnat, Dyer, & Ghosh

The last equation follows from Eq. (30).

Let S and S denote B R N (C(σ), R(σ))∩aff(σ) and B R N (C(σ, ω σ ), R(σ, ω σ )) respectively. From Eq. ( 30) and (31), we have d R N (∂S , ∂S) ≤ η 1 + 2η 2 . This implies that there exists q ∈ ∂S such that qq ≤ 2η 2 .

(32) Therefore from Eq. ( 29) and (32), we get

Using the facts that δ 2 0 ≤ Γ m+1 0 ≤ Γ 2 0 (from hyp. (4) of the lemma and

≤ 2 7 3 2 β 2 Γ 0 (from Eq. ( 25)), and Eq. ( 29) and (30), we get

where

We will use the following lemma from [START_REF] Boissonnat | Triangulating smooth submanifolds with light scaffolding[END_REF] to bound the volume of F (σ).

Lemma C.3 Let p be a point on M. There exists ξ that depends on rch(M) and m, and A that depends only on m such that, for all r = t ≤ ξ, we have

where V m is the volume of the m-dimensional unit Euclidean ball.

Proof of Lemma 5.8 For the rest of the proof we define

Consider the following elementary weight function:

Let p be a vertex of σ. Let c be the point closest to C(σ) on T p M and c * be the point closest to c on M (see Fig. 5).

Inria

Constructing intrinsic

From Lemma 2.1, and the facts that Υ(σ) ≥ Γ m 0 and ∆(σ) ≤ 2R(σ) ≤ 4βt, we have

and from Lemma 4.2

Let x ∈ F (σ, t) and x * be the point closest to

Using the fact that C(σ)x * = R(σ), we get Therefore

The last inequality follows from the fact that (1 + x) m -(1x) m ≤ 2 m x for all x ∈ [0, 1]. From hyp. 2 and the fact that Γ 0 < 1, we have

The lemma now follows from Eq. ( 37) and (38).