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

Ramsay Dyer
email: ramsay.dyer@inria.fr

Arijit Ghosh

Delaunay triangulation of manifolds

Keywords: Delaunay complex, triangulation, manifold, protection, perturbation MSC2010: 57R05 (Primary), 52B70, 54B15 (Secondary)

We present an algorithm for producing Delaunay triangulations of manifolds. The algorithm can accommodate abstract manifolds that are not presented as submanifolds of Euclidean space. Given a set of sample points and an atlas on a compact manifold, a manifold Delaunay complex is produced for a perturbed point set provided the transition functions are bi-Lipschitz with a constant close to 1, and the original sample points meet a local density requirement; no smoothness assumptions are required. If the transition functions are smooth, the output is a triangulation of the manifold.

The output complex is naturally endowed with a piecewise flat metric which, when the original manifold is Riemannian, is a close approximation of the original Riemannian metric. In this case the output complex is also a Delaunay triangulation of its vertices with respect to this piecewise flat metric.

Introduction

We present an algorithm for computing Delaunay triangulations of Riemannian manifolds. Not only is this the first algorithm guaranteed to produce a Delaunay triangulation of an arbitrary compact Riemannian manifold, it also provides the first theoretical demonstration of the existence of such triangulations on manifolds of dimension greater than 2 with nonconstant curvature.

The Delaunay complex is a natural structure to consider when seeking to triangulate a space equipped with a metric. It plays a central role in the development of algorithms for meshing Euclidean domains. In applications where an anisotropic mesh is desired, a standard approach is to consider a Riemannian metric defined over the domain and to construct an approximate Delaunay triangulation with respect to this metric [START_REF] Labelle | Anisotropic Voronoi diagrams and guaranteed-quality anisotropic mesh generation[END_REF]BWY15,CG12]. In this context we can consider the domain to be a Riemannian manifold that admits a global coordinate parameterisation. The algorithm we present here encompasses this setting, modulo boundary considerations.

In the case of surfaces, it has been shown that a Delaunay triangulation exists if the set of sample points is sufficiently dense [START_REF] Leibon | Random Delaunay triangulations, the Thurston-Andreev theorem, and metric uniformization[END_REF]DZM08]. What is perhaps surprising is that, contrary to previous claims [START_REF] Leibon | Delaunay triangulations and Voronoi diagrams for Riemannian manifolds[END_REF], this is no longer true in higher dimensional manifolds; sample density alone is not sufficient to ensure that the Delaunay complex is a triangulation [BDG13a, App. A].

In Euclidean space, R m , the Delaunay complex on a set of points P is a triangulation provided the points are generic, i.e., no ball empty of points contains more than m + 1 points of P on its boundary [START_REF] Delaunay | Sur la sphère vide[END_REF]. A point set that is not generic is said to be degenerate, and such configurations can be avoided with an arbitrarily small perturbation. However, when the metric is no longer homogeneous, an arbitrarily small perturbation is not sufficient to guarantee a triangulation. In previous work [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] we have shown that genericity can be parameterised, with the parameter, δ, indicating how far the point set is from degeneracy. A δ-generic point set in Euclidean space yields a Delaunay triangulation that is quantifiably stable with respect to perturbations of the metric, or of the point positions. We later produced an algorithm [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] that, given an initial point set P ⊂ R m , generates a perturbed point set P that is δ-generic. The algorithm we present here adapts this Euclidean perturbation algorithm to the context of compact manifolds equipped with a metric that can be locally approximated by a Euclidean metric. In particular, this includes Riemannian metrics, as well as the extrinsic metric on submanifolds, which defines the so-called restricted Delaunay complex, variations of which have been exploited in algorithms for reconstructing submanifolds of Euclidean space from a finite set of sample points [CDR05, [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF].

The simplicial complex produced by our algorithm is naturally equipped with a piecewise flat metric that is a quantifiably good approximation to the metric on the original manifold. The stability properties of the constructed Delaunay triangulation yield additional benefits. In particular, the produced complex is a Delaunay triangulation of its vertices with respect to its own intrinsic piecewise-flat metric: a self-Delaunay complex. Such complexes are of interest in discrete differential geometry because they provide a natural setting for discrete exterior calculus [START_REF] Bobenko | A discrete Laplace-Beltrami operator for simplicial surfaces[END_REF][START_REF] Dyer | Self-Delaunay meshes for surfaces[END_REF]HKV12].

In Section 2 we review the main ideas involved in the perturbation algorithm [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] for producing δ-generic point sets in Euclidean space. The extension of the algorithm to general manifolds is described in Section 3, where we also state our main results. Sections 3.3 and 4 describe the analysis that leads to these results. All constants are explicitly expressed but they are generally, and somewhat unconventionally, rounded to powers of 2.

Contributions In this paper we provide the first proof of existence of Delaunay triangulations of arbitrary compact Riemannian manifolds. The proposed algorithm is the first triangulation algorithm that can accommodate abstract manifolds that are not presented as submanifolds of Euclidean space. (Although Nash's embedding theorem ensures that all Riemannian manifolds may be realised as submanifolds of Euclidean space, constructing such an embedding is considerably more complicated than constructing local coordinate charts.) The output complex is a good geometric approximation to the original manifold, and also possesses the self-Delaunay property. These results are summarised in Theorem 3.

The algorithm accommodates more general inputs than Riemannian manifolds: strong bi-Lipschitz constraints are required on the transition functions, but they need not be smooth. In this case we cannot guarantee a triangulation, but, as stated in Theorems 1 and 2, the output complex is a manifold Delaunay complex (it is the nerve of the Voronoi diagram of the perturbed points P ⊂ M, and it is a manifold).

The framework encompasses and unifies previous algorithms for constructing anisotropic meshes, and for meshing submanifolds of Euclidean space, and the algorithm itself is conceptually simple (if not the analysis).

The perturbation strategy

We outline here the main ideas behind the Euclidean perturbation algorithm [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF], upon which the current algorithm is based. Given a set P ⊂ R m , that algorithm produces a perturbed point set P that is δ-generic. This means that the Delaunay triangulation of P will not change if the metric is distorted by a small amount [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF].

Thickness and protection We consider a finite set P ⊂ R m . A simplex σ ⊂ P is a finite collection of points: σ = {p 0 , . . . , p j }, where j is the dimension of σ (one less than the number of points in σ). We work with abstract simplices, and in particular x ∈ σ means x is a vertex of σ. Although we prefer abstract simplices, we freely talk about standard geometric properties, such as the longest edge length, ∆(σ), and the length of the shortest edge L(σ).

For p ∈ σ, σ p is the facet opposite p, and D(p, σ) is the altitude of p in σ, i.e., D(p, σ) = d(p, aff(σ p )), where d(p, q) is the standard Euclidean distance between p, q ∈ R m , and for A ⊂ R m , d(p, A) = inf x∈A d(p, x). The thickness of σ is a measure of the quality of σ, and is denoted Υ(σ). If σ is a 0-simplex, then Υ(σ) = 1. Otherwise Υ(σ) is the smallest altitude of σ divided by j∆(σ), where j is the dimension of σ. The factor of j in the denominator was introduced in [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] to simplify certain bounds, and we continue to employ it here for convenience. If Υ(σ) = 0, then σ is degenerate. We say that σ is Υ 0 -thick, if Υ(σ) ≥ Υ 0 . If σ is Υ 0 -thick, then so are all of its faces.

A circumscribing ball for a simplex σ is any m-dimensional ball that contains the vertices of σ on its boundary. A degenerate simplex may not admit any circumscribing ball. If σ admits a circumscribing ball, then it has a circumcentre, C(σ), which is the centre of the unique smallest circumscribing ball for σ. The radius of this ball is the circumradius of σ, denoted R(σ).

For any set A ⊂ R m , A denotes its closure, ∂A its boundary, conv(A) its convex hull, and aff(A) its affine hull. A ball B(c, r) ⊂ R m is open (and B(c, r) is its closure). The Delaunay complex, Del(P) is the (abstract) simplicial complex defined by the criterion that a simplex belongs to Del(P) if it has a circumscribing ball whose intersection with P is empty. For p ∈ P, the star of p is the subcomplex star(p; Del(P)) consisting of all simplices that contain p, as well as the faces of those simplices. An m-simplex σ

m ∈ Del(P) is δ-protected if B(C(σ m ), R(σ m ) + δ) ∩ P = σ m . The point set P ⊂ R m is δ-generic if all the m-simplices in Del(P) are δ-protected.

Forbidden configurations

The essential observation that leads to the perturbation algorithm is that if P ⊂ R m is such that there exists a Delaunay m-simplex that is not δ-protected, then there is a forbidden configuration: a (possibly degenerate) simplex τ ⊂ P characterised by the properties we describe in Lemma 2 below. We emphasise that a forbidden configuration need not be a Delaunay simplex. The perturbation algorithm guarantees that the Delaunay m-simplices will be δ-protected by ensuring that each point is perturbed to a position that is not too close to the circumsphere of any of the nearby simplices in the current (perturbed) point set. A volumetric argument shows that this can be achieved.

Let D ⊂ R m be a bounded set, and P ⊂ R m a finite set of points. We say P is -dense for D if d(x, P) < for all x ∈ D. We refer to as the sampling radius. For our purposes, D will usually be a domain (open and simply connected), or the closure of a domain. The set P is µ 0 -separated if d(p, q) ≥ µ 0 for all p, q ∈ P, and P is a (µ 0 , )-net for D if it is µ 0 -separated, and -dense for D. If no set D is explicitly specified, P is a (µ 0 , )-net if it is a (µ 0 , )-net for

D (P) = {x ∈ conv(P) | d(x, ∂conv(P)) ≥ }.
(1)

Note that the bounds indicated by the sampling parameters (µ 0 , ) are not assumed to be tight. In particular, if P is a (µ 0 , )-net, then it is also a (μ 0 , ˜ )-net for any μ0 ≤ µ 0 and ˜ ≥ . In this work we always assume µ 0 ≤ 1, even though a larger value is possible. This does not limit the point sets that can be considered, although it does mean that some of the bounds could be tightened for very well-separated point sets.

Given an initial (µ 0 , )-net P, the goal is to produce a perturbed point set P that contains no forbidden configurations. A ρ-perturbation of a (µ 0 , )-net P ⊂ R m is a bijective application ζ : P → P ⊂ R m such that d(ζ(p), p) ≤ ρ for all p ∈ P. Unless otherwise specified, a perturbation will always refer to a ρ-perturbation, with ρ = ρ 0 for some ρ 0 ≤ µ 0 4 . We also refer to P itself as a perturbation of P. We generally use p to denote the point ζ(p) ∈ P , and similarly, for any point q ∈ P we understand q to be its preimage in P. We observe [BDG14, Lemma 2.2] that P is necessarily a (µ 0 , )-net, with ≤ 5 4 and µ 0 ≥ 2 5 µ 0 . Henceforth, we shall always assume these inequalities:

0 < ≤ ≤ 5 4 1 ≥ µ 0 ≥ µ 0 ≥ 2 5 µ 0 .
(2) This is not restrictive, even though it is possible that P could satisfy a sampling radius smaller than and a separation parameter larger than µ 0 . In this case P is still a (µ 0 , )-net with (µ 0 , ) satisfying (2). Given a positive parameter Γ 0 ≤ 1, we say that σ is Γ 0 -good if for all j-dimensional faces σ j ⊆ σ, we have Υ(σ j ) ≥ Γ j 0 , where Γ j 0 is the j th power of Γ 0 . A Γ 0 -flake is a simplex that is not Γ 0 -good, but whose facets all are. A flake may be degenerate. The altitudes of a flake are all subjected to an upper bound proportional to Γ 0 .

If a simplex is not Γ 0 -good, then it necessarily contains a face that is a flake. This follows easily from the observation that Υ(σ) = 1 if σ is a 1-simplex. If σ m ∈ Del(P) is not δ-protected, then there is a q ∈ P \ σ m that is within a distance δ of the circumsphere of σ m . Since {q} ∪ σ m is (m + 1)-dimensional, it is degenerate, and therefore has a face τ that is a Γ 0 -flake; it may be that τ = {q} ∪ σ m (i.e., not a proper face). Such a τ is a forbidden configuration.

If a simplex σ has a circumcentre, we define the diametric sphere as the boundary of the smallest circumscribing ball: S m-1 (σ) = ∂B(C(σ), R(σ)), and the circumsphere: S(σ) = S m-1 (σ) ∩ aff(σ). If σ ⊂ τ , then S(σ) ⊆ S(τ ), and if dim σ = m and σ is nondegenerate, then S(σ) = S m-1 (σ).

The bound on the altitudes, together with the stability property of circumscribing balls of thick simplices, allows us to demonstrate (subject to

S(σ) σ p ≤ α 0 R(σ)
Figure 1: The α 0 -hoop property illustrated for a tetrahedron τ = {p} ∪ σ. The vertex p lies within a distance α 0 R(σ) from the circumcircle of its opposite facet. This property applies to any vertex q ∈ τ . constraints on δ 0 and Γ 0 ) that forbidden configurations have the α 0 -hoop property. A k-simplex τ has the α 0 -hoop property if for every (k -1)-facet σ ⊂ τ we have

d(p, S(σ)) ≤ α 0 R(σ) < ∞,
where p is the vertex of τ not in σ (see Figure 1), and α 0 is a constant that depends on Γ 0 (see Lemma 2 P1 below). The idea of the perturbation algorithm is to ensure that no simplicies with the α 0 -hoop property exist, thus ensuring that there are no forbidden configurations. We are concerned with forbidden configurations in the perturbed point set P . In addition to the two parameters that describe a (µ 0 , )-net, forbidden configurations depend on the flake parameter Γ 0 , as well as the parameter δ 0 , which governs the protection via the requirement δ = δ 0 µ 0 .

Definition 1 (Forbidden configuration) Let P ⊂ R m be a (µ 0 , )-net. A (k + 1)-simplex τ ⊆ P , is a (δ 0 , Γ 0 )-forbidden configuration in P if it is a Γ 0 -flake, with k ≤ m, and there exists a p ∈ τ such that τ p has a circumscribing ball B = B(C, R) with R < , and |d(p, C) -R| ≤ δ, where δ = δ 0 µ 0 . When the parameters (δ 0 , Γ 0 ) can be inferred from context, we simply speak of a forbidden configuration.

Notice that if τ is a Γ 0 -flake that has a circumradius less than , then it is automatically a forbidden configuration. Thus the removal of forbidden configurations ensures that the Delaunay triangulation of P has only Γ 0 -good simplices of all dimensions, and further, that the m-simplices are δ-protected, where δ = δ 0 µ 0 . Definition 1 itself is awkward, but for most purposes we can simply refer to the following summary [BDG14, Theorem 3.10] of properties of forbidden configurations in P in terms of the parameters of the original point set P. In particular, forbidden configurations have the α 0 -hoop property for an α 0 that depends on Γ 0 .

Lemma 2 (Properties of forbidden configurations) Suppose that P ⊂ R m is a (µ 0 , )-net and that P is a ρ 0 -perturbation of P, with ρ 0 ≤ µ 0 4 . If

δ 0 ≤ Γ m+1 0 and Γ 0 ≤ 2µ 2 0 75 , (3)
then every (δ 0 , Γ 0 )-forbidden configuration τ ⊂ P satisfies all of the following properties:

P1 Simplex τ has the α 0 -hoop property, with α 0 = 2 13 Γ 0

µ 3 0 . P2 For all p ∈ τ , R(τ p ) < 2 . P3 ∆(τ ) < 5 2 (1 + 1 2 δ 0 µ 0 ) . P4 Every facet of τ is Γ 0 -good.
The algorithm focuses on Property P1 of forbidden configurations. The bound on Γ 0 imposed in (3) is much larger than the bound required by the algorithm; the latter bound is proportional to µ m 2 0 (see [BDG14, Thm. 5.6], and also Hypothesis 2 and Lemma 5 below). A critical aspect of the α 0 -hoop property is its symmetric nature; if we can ensure that τ has one vertex that is not too close to its opposite facet, then τ cannot be a forbidden configuration.

Using Property P3, we can find, for each p ∈ P, a complex S p consisting of all simplices σ ∈ P such that after perturbations {p} ∪ σ could be a forbidden configuration.

The algorithm proceeds by perturbing each point p ∈ P in turn, such that each point is only visited once. The perturbation p → p is found by randomly trying perturbations p → x until it is found that x is a good perturbation. A good perturbation is one in which d(x, S(σ)) > 2α 0 for all σ ∈ S p (P ), where S p (P ) is the complex in the current perturbed point set whose simplices correspond to those in S p . By Property P1, {x} ∪ σ cannot be a forbidden configuration.

Assuming a sufficiently small Γ 0 , a volumetric argument based on the finite number of simplices in S p , the small size of α 0 , and the volume of the ball B(p, ρ 0 ) of possible perturbations of p, reveals a high probability that p → x will be a good perturbation, and thus ensures that the algorithm will terminate.

Upon termination there will be no forbidden configurations in P , because every perturbation p → p ensures that there are no forbidden configurations incident to p in the current point set, and no new forbidden configurations are introduced.

Overview and main results

The extension of the perturbation algorithm to the curved setting is accomplished by performing the perturbations, and the analysis, in local Euclidean coordinate patches. The main idea is that forbidden configurations exhibit some stability with respect to small changes in the Euclidean metric. By this we mean that if τ is a (δ 0 , Γ 0 )-forbidden configuration in one coordinate patch, then it will appear as a ( δ0 , Γ0 )-forbidden configuration in a different coordinate patch, where δ0 and Γ0 are close to δ 0 and Γ 0 respectively, assuming that the transition function has low metric distortion. In particular, if we ensure there are no forbidden configurations in some region of one Euclidean coordinate patch, then, assuming a slightly smaller hoop parameter α 0 , there will be no forbidden configurations in the corresponding region of any nearby coordinate patch. This means that the perturbed point set will be δ-generic in any local Euclidean coordinate patch, and the resulting stability of the local Euclidean Delaunay triangulations ensures that they will agree on neighbouring patches.

We assume we have a finite set of points P in a compact manifold M. It is convenient to employ an index set N of unique (integer) labels for P, thus we employ a bijection ι : N → P ⊂ M. We assume that P is sufficiently dense that we may define an atlas {(W i , ϕ i )} i∈N for M such that the coordinate charts ϕ i : W i → U i ⊂ R m have low metric distortion, as defined in Section 3.2. The set W i is required to contain a sufficiently large ball centred at the point indexed by i. We refer to U i as a coordinate patch.

We work exclusively in the Euclidean coordinate patches U i , exploiting the transition functions ϕ ji = ϕ j • ϕ -1 i to translate between them. We define P i = ϕ i (W i ∩ P), but given these sets, the algorithm itself makes no explicit reference to either P or to the coordinate charts ϕ i , except to keep track of the labels of the points. We employ the discrete map φ i = ϕ i • ι to index the elements of the set P i .

The idea is to perturb p i = φ i (i) ∈ U i in such a way, p i → p i , that not only are there no forbidden configurations incident to p i in P i = ϕ i (W i ∩ P ), but there are no forbidden configurations incident to ϕ ji (p i ) ∈ P j ⊂ U j either, where j is the index of any sample point near p i .

Before detailing the requirements of the input data in Section 3.2, we briefly discuss the implicit and explicit properties of the underlying manifold M in Section 3.1. A summary of the analysis and main results is presented in Section 3.3.

Manifolds represented by transition functions

The essential input data for the algorithm are the transition functions, and the sample points in the coordinate patches; we do not explicitly use the coordinate charts or the metric on the manifold. However, given that the transition functions can be defined by an atlas on a manifold, this manifold is essentially unique: If M has an atlas {( Wi , φi )} i∈N such that φi ( Wi ) = U i and ϕ ji = φj • φ-1 i for all i, j ∈ N , then M and M are homeomorphic. Indeed, we define the map f :

M → M by f (x) = φ-1 i • ϕ i (x) if x ∈ W i . The map is well defined, because ϕ j = ϕ ji • ϕ i on W i ∩ W j , and φ-1 j = φ-1 i • ϕ -1 ji on U ji = ϕ j (W i ∩ W j ).
It can be verified directly from the definition that f is a homeomorphism, since it is bijective and locally a homeomorphism. Although the algorithm does not make explicit reference to a metric on the manifold M, the metric distortion bounds required on the transition functions imply a metric constraint. Implicitly we are using a metric on the manifold for which the coordinate charts have low metric distortion. If a metric on the manifold is not explicitly given then, at least in the case where the transition functions are smooth, we can be sure that such a metric exists: Given the coordinate charts an appropriate Riemannian metric on the manifold can be obtained from the coordinate patches by the standard construction employing a partition of unity subordinate to the atlas (e.g., [Boo86, Thm. V.4.5]).

Thus, although the manifold may be presented abstractly in terms of coordinate patches and transition functions between them, this information essentially characterises the manifold. The algorithm we present is not a reconstruction algorithm; it is an algorithm to triangulate a known manifold.

The setting and input data

We take as input a finite index set N = {1, . . . , n}, which we might think of as an abstract set of points (without geometry), together with the geometric data we will now introduce. The details of the arguments that lead to our choices in the size of the domains are given in Section 4.3.

Coordinate patches.

For each i ∈ N we have a neighbourhood set N i ⊆ N , a sampling radius i > 0, and an injective application

φ i : N i → U i ⊆ R m , such that P i = φ i (N i ) is a (µ 0 , i )-net for B(p i , 8 i ) ⊆ U i , where p i = φ i (i).
The separation parameter µ 0 is globally defined to be the same on all coordinate patches, but the sampling radius i may be different in different patches, subject to a mild constraint described below. We call the standard metric on U i ⊆ R m the local Euclidean metric for i, and we will denote it by d i to distinguish between the different local Euclidean metrics. Similarly, B i (c, r) denotes a ball with respect to the metric d i .

Transition functions. For each

p j ∈ B i (p i , 6 i ) ⊂ U i we require a neigh- bourhood U ij ⊆ U i such that B i (p i , 6 i ) ∩ B i (p j , 9 i ) ⊆ U ij and U ij ∩ φ i (N i ) = φ i (N i ∩ N j ).
The set U ij is the domain of the transition function ϕ ji , which is a homeomorphism

ϕ ji : U ij → U ji ⊆ U j , such that ϕ ji = ϕ -1 ij and ϕ ji • φ i = φ j on N i ∩ N j .
These transition functions are required to have low metric distortion:

|d i (x, y) -d j (ϕ ji (x), ϕ ji (y))| ≤ ξ 0 d i (x, y) for all x, y ∈ U ij , (4)
where ξ 0 > 0 is a small positive parameter that quantifies the metric distortion. We say that ϕ ji is a ξ 0 -distortion map.

In order to ease the notational burden, φ i (j) ∈ U ik , and φ k (j) ∈ U ki , are denoted by the same symbol, p j . Ambiguities are avoided by distinguishing between the Euclidean metrics d i and d k . Although d i is the canonical metric on U ik , we may consider the pullback of d k from the homeomorphic domain U ki . Thus for x, y ∈ U ik the expression d k (x, y) is understood to mean d k (ϕ ki (x), ϕ ki (y)), but we also occasionally employ the latter, redundant, notation.

Using symmetry, we observe that Equation (4) implies that

|d i (x, y) -d j (x, y)| ≤ ξ 0 min{d i (x, y), d j (x, y)}.
Our analysis will require that ξ 0 be very small. For standard coordinate charts, ξ 0 can be shown to be O( ), where is a sampling radius on the manifold. For example, this is the case when considering a smooth submanifold of R N , and using the orthogonal projection onto the tangent space as a coordinate chart [BDG13a, Lemma 3.7]. Thus ξ 0 may be made as small as desired by increasing the sampling density.

Adaptive sampling. We will further require a constraint on the difference between neighbouring sampling radii:

| i -j | ≤ 0 min{ i , j },
whenever d i (p i , p j ) ≤ 6 i .
This allows us to work with a constant sampling radius in each coordinate frame, while accommodating a globally adaptive sampling radius. For example, suppose : M → R is a positive, ν-Lipschitz function, with respect to the metric d M on the manifold. Then may be used as an adaptive sampling radius on M, i.e., P ⊂ M is -dense if d M (x, P) < (x) for all x ∈ M. A popular example of such a function is (x) = νf (x), where f is the (1-Lipschitz) local feature size [START_REF] Amenta | Surface reconstruction by Voronoi filtering[END_REF].

Using the ν-Lipschitz continuity of , we can define, for any p i ∈ M, a constant ˜ i , such that P is ˜ i -dense in some neighbourhood of p i . In fact, given c > 0, with c < ν -1 , we find that P is ˜ i -dense within the ball B M (p i , c˜ i ), where

˜ i = (p i ) 1 -cν .
For any

p j ∈ B M (p i , c˜ i ), we obtain |˜ i -˜ j | ≤ 0 ˜ i , where 0 = cν 1 -cν , (5)
and if ν ≤ 1 2c , then 0 ≤ 1. Similarly, if P is μ0 -separated, i.e., if d M (p, q) ≥ μ0 max{ (p), (q)} for all p, q ∈ P, then it will be μ0 ˜

i -separated on B M (p i , c˜ i ), provided μ0 ≤ (1 -2cν)μ 0 , which is positive if ν < 1/(2c).
In our framework here, the local constant sampling radii are applied to the local Euclidean metric, rather than the metric on the manifold, but the same idea applies. Although Equation (5) indicates that 0 is expected to become small as the sampling radius decreases, our analysis does not demand this. As explained in Section 3.3, we only require that 0 be mildly bounded.

We summarise the assumptions on the input to the extended algorithm as Hypotheses 1, where the parameters 0 and ξ 0 are left free to be constrained by subsequent hypotheses.

Hypotheses 1 (Input assumptions) We have a finite index set N representing the sample points. For each i ∈ N there is associated a subset of neighbours N i ⊆ N . The geometry is imposed by 1. Coordinate patches. For each i ∈ N , there is a coordinate patch U i ⊆ R m , and an injective application φ i :

N i → U i such that P i = φ i (N i ) is a (µ 0 , i )-net for B i (p i , 8 i ) ⊆ U i .
We introduce a parameter 0 ≥ 0, and demand that, if d i (p i , p j ) ≤ 6 i , then

| i -j | ≤ 0 min{ i , j }. 2. Transition functions. Each p j ∈ B i (p i , 6 i ) ⊂ U i lies in the domain U ij ⊆ U i of the transition function ϕ ji : U ij ∼ =
-→ U ji . The domains must be suffiently large:

B i (p i , 6 i ) ∩ B i (p j , 9 i ) ⊆ U ij ,
and the transition functions must satisfy the compatibility conditions ϕ ji = ϕ -1 ij and ϕ ji • φ i = φ j on N i ∩ N j . Furthermore, the metric distortion of the transition functions is bounded by a parameter ξ 0 :

|d i (x, y) -d j (ϕ ji (x), ϕ ji (y))| ≤ ξ 0 d i (x, y) for all x, y ∈ U ij .
The extended algorithm The algorithm we present here is the same in spirit as the algorithm for the Euclidean setting [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] described in Section 2, and we refer to it as the extended algorithm. It takes an input satisfying Hypotheses 1. For each i ∈ N , a ρ 0 i -perturbation is repeatedly applied to the point p i ∈ P i ⊂ U i until a good perturbation p i is found. The Definition 3 of a good perturbation involves something closely resembling the hoop property P1 with a parameter α0 > α 0 . When a good perturbation p i is selected, then the affected point sets P j are updated, as well as P i itself. By demonstrating stability of the hoop property with respect to small changes in the Euclidean metric, we are able to show that when the extended algorithm terminates, there will be no forbidden configurations in the region of interest of any local Euclidean coodinate patch. Then, assuming appropriate constraints on ξ 0 and 0 , a manifold simplicial complex whose vertex set is N is constructed by defining the star of i to correspond to star(p i ; Del(P i )). The stability of these stars ensures that this complex is indeed a manifold (Theorem 1), and we call it Del(P ), as justified by Theorem 2. We refer to Del(P ) as the output of the extended algorithm, thus we assume that the extended algorithm includes a final step of computing all the stars after the perturbation algorithm has completed.

Outline of the analysis

We have defined the point sets P i = φ i (N i ) in the coordinate patch for p i . We will let P i denote the corresponding perturbed point set at any stage in the algorithm: P i changes during the course of the algorithm, and we do not rename it according to the iteration as was done in the original description of the algorithm for flat manifolds [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF]. The perturbation of a point p i → p i is performed in the coordinate patch U i , and then all the discrete maps must be updated so that if i ∈ N j , then φ j (i) = ϕ ji (p i ). However, we will refer to the point as p i regardless of which coordinate frame we are considering. The discrete maps φ i will change as the algorithm progresses, but φ i will always refer to the initial map.

In order to ensure that we maintain a (µ 0 , j )-net in each local Euclidean coordinate patch U j , we need to constrain the point perturbation so that the cumulative perturbation is a ρ0 j -perturbation with ρ0

≤ µ 0 /4 (see [BDG14, Lemma 2.2]). If the perturbation p i → ζ(p i ) is such that d i (ζ(p i ), p i ) ≤ ρ = ρ 0 i , then d j (ζ(p i ), p i ) ≤ (1 + ξ 0 )ρ ≤ (1 + ξ 0 )(1 + 0 )ρ 0 j , and we have ρ0 = (1 + ξ 0 )(1 + 0 )ρ 0 .
Thus we demand that ρ0 ≤ µ 0 4 .

In order to facilitate the analysis, we want an explicit constant to bound the ratio between ρ 0 and ρ0 . We ensure that

(1 + 0 )(1 + ξ 0 ) ≤ 2 (7)
by imposing the mild constraint that

0 ≤ 1 -ξ 0 1 + ξ 0 . (8)
We will keep the definition of forbidden configuration as in the flat case. In other words a forbidden configuration is that which satisfies the four properties described in Lemma 2, where refers to the local sampling radius i .

We do not attempt to remove the forbidden configurations from all of P i . Rather, we define Q i = P i ∩ B i (p i , 6 i ) as our region of interest. The reasoning behind this choice appears in Section 4.3, where we also show that [BDG14, Lemma 3.6] implies:

Lemma 17 (Protected stars) If there are no forbidden configurations in Q i , then all the m-simplices in star(p i ; Del(Q i )) are Γ 0 -good and δ-protected, with δ = δ 0 µ 0 i . This allows us to exploit the Delaunay metric stability result [BDG13b, Theorem 4.17], which we show (Section 4.3) may be stated in our current context as:

Lemma 18 (Stable stars) If ξ 0 ≤ Γ 2m+1 0 µ 2 0 2 12 ,
and there are no forbidden configurations in Q i , then for all p j ∈ star(p i ; Del(P i )), we have star(p i ; Del(P i )) ∼ = star(p i ; Del(P j )).

The main technical result we develop in the current analysis is the bound on the distortion of the hoop property (Section 4.4) due to the transition functions:

Lemma 20 (Hoop distortion) If ξ 0 ≤ Γ 2m+1 0 4 2 ,
then for any forbidden configuration τ = {p j } ∪ σ ⊂ Q i , there is a simplex σ = ϕ ji (σ) ⊂ P j such that d j (p j , S m-1 (σ)) ≤ 2 α0 j , where

α0 = 2 16 m 3 2 Γ 0 µ 3 0 .
The idea of the extended algorithm is to perturb each point such that the simplices σ described by Lemma 20 do not exist. This ensures that there will be no forbidden configurations in any of the point sets Q i .

The proof of Lemma 20 relies heavily on the thickness bound (Property P4) for the facets of a forbidden configuration. In Section 4.1 we show bounds on the changes of the intrinsic properties, such as thickness and circumradius, of a Euclidean simplex subjected to the influence of a transition function. This leads, as shown in Section 4.2, to bounds on circumcentre displacement under small changes of a Euclidean metric. These bounds could not be recovered directly from earlier work [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF], because they involve simplices that are not full dimensional. With these results in place, the proof of Lemma 20 is assembled in Section 4.4.

By considering the diameter of a forbidden configuration subjected to metric distortion, we can determine the size of the neighbourhood of p i that must be considered when checking whether a perturbation p i → p i creates conflicts.

Suppose τ ⊂ Q j ⊂ U j is a forbidden configuration with p i ∈ τ . By Lemma 2, Property P3, we have ∆(τ )

< 5 2 (1 + 1 2 δ 0 µ 0 ) j . It follows then that if τ = ϕ ij (τ ) ⊂ P i , then ∆(τ ) < (1 + ξ 0 )(1 + 0 ) 5 2 1 + 1 2 δ 0 µ 0 i ≤ 5 1 + 1 2 δ 0 µ 0 i ,
and we find, as in [BDG14, Lemma 4.4], that if δ 0 ≤ 2 5 , then

(φ j ) -1 (τ ) ⊂ φ -1 i (B i (p i , r) ∩ P i ), where r = 5 + 3µ 0 2 i .
Indeed, this is ensured by the fact that P i is a (µ 0 , i )-net for B i (p i , 8 i ), and 8 ir > i .

Let S i denote all the m-simplices in N i whose vertices are contained in

φ -1 i (B i (p i , r) ∩ P i ) \ {i} where r = 5 + 3µ 0 2 i .
Then the simple packing argument demonstrated in [BDG14, Lemma 5.1] yields

#S i < 14 µ 0 m 2 +m . (9)
We strengthen the definition of a good perturbation:

Definition 3 (Good perturbation) For the extended algorithm, we say that p i → x is a good perturbation of p i ∈ U i if there are no simplices σ ∈ φ i (S i ) of dimension ≤ m such that d i (x, S m-1 (σ)) ≤ 2α 0 i , where α0 is defined in Lemma 20.

Lemma 20 bounds the distance from p j to the diametric sphere S m-1 (σ). Using [BDG14, Lemma 3.14], this will yield a bound on the distance to the circumsphere: d j (p j , S(σ)) ≤ 2α 0 j . If such an α0 were used in Definition 3 instead of α0 , then it would be sufficient to only consider the m-simplices, because if σ is a non-degenerate j-simplex, with j < m, then it is the face of some non-degenerate m-simplex τ , and S(σ) ⊂ S m-1 (τ ).

Using Definition 3 for a good perturbation, the extended algorithm yields the analogue of [BDG14, Lemma 4.3]: Lemma 4 After the extended algorithm terminates there will be no forbidden configurations in Q i , for every i ∈ N .

Proof

We argue by induction that after the i th iteration, for any j ≤ i, and any k ∈ N , there are no forbidden configurations in Q k that have p j as a vertex. For i = 1, the assertion follows directly from Definition 3, and Lemma 20. Assume the assertion is true for i -1. Suppose τ is a forbidden configuration in Q k , after the i th iteration. Then since p i is a good perturbation, according to Definition 3, τ cannot contain p i . Also, τ cannot contain any p j with j < i, for that would contradict the induction hypothesis. Thus the hypothesis holds for all i ≥ 1.

In order to quantify the conditions under which the algorithm is guaranteed to terminate, we need to show that there is a positive probability that any given perturbation will be a good perturbation, i.e., that the volume occupied by the good perturbations of any p i is a significant fraction of vol B i (p i , ρ 0 i ) (which is the volume of possible perturbations). We use the same volumetric analysis that is demonstrated in the proof of [BDG14, Lemma 5.4], with the only modifications being a change in two of the constants involved in the calculation. In particular, the number of simplices involved is now given by Equation ( 9), and we use the bound on α0 given by Lemma 20, which is 2 3 m 3 2 times the bound on α 0 used in the original calculation. This calculation, coupled with the criterion for Lemma 20, yields a constraint on ξ 0 with respect to the perturbation parameter ρ 0 . Together with Equations ( 6) and (8), this gives us all the constraints on the parameters that will ensure the existence of good perturbations, and therefore the termination of the algorithm:

Hypotheses 2 (Parameter constraints) Define ρ0 = (1 + 0 )(1 + ξ 0 )ρ 0 , We require 0 ≤ 1 -ξ 0 1 + ξ 0 , and ρ0 ≤ µ 0 4 , and
ξ 0 ≤ 1 2 4 ρ 0 C 4m+2 ,
where

C = m 3 2 2 µ 0 4m 2 +5m+21
.

Thus, using Lemma 17, the main result [BDG14, Theorem 4.1] of the original perturbation algorithm can be adapted to the context of the extended algorithm as:

Lemma 5 (Algorithm guarantee) If Hypotheses 1 and 2 are satisfied, then the extended algorithm terminates, and for every i ∈ N , the set Q i is a (µ 0 , i )-net such that there are no forbidden configurations with

Γ 0 = ρ 0 C , and δ = Γ m+1 0 µ 0 i ,
where µ 0 = µ 0 -2ρ 0 1+ρ 0 , and i = (1 + ρ0 ) i .

This allows us to apply Lemma 18, and we can define the abstract complex Del(P ) by the criterion that φ i (star(i; Del(P ))) = star(p i ; Del(P i )) for all i ∈ N . This is a manifold piecewise linear1 simplicial complex. The bound on ξ 0 imposed by Lemma 18 is met by the one imposed by Hypotheses 2 and we arrive at our first main result:

Theorem 1 (Manifold mesh) Given an input satisfying Hypotheses 1 and 2, the extended algorithm produces a manifold abstract simplicial complex Del(P ) defined by star(i; Del(P )) ∼ = star(p i ; Del(P i )).

The algorithm itself makes no explicit reference to the underlying manifold M or point set P ⊂ M, but we need to consider these in order to justify the name Del(P ) for the output of the extended algorithm (see Theorem 2 below). Also, without further assumptions, we can provide no guarantee that Del(P ) is homeomorphic to M. When the transition functions are smooth (so that, if necessary, a Riemannian metric can be constructed on M, as mentioned on page 9), then a natural homeomorphism |Del(P )| → M is given by the barycentric coordinate map [DVW15], provided a sufficient initial sampling density is met (see Theorem 3 below).

Given P ⊂ M, we define the set P ⊂ M to be the perturbed point set produced by the algorithm, i.e., P → P is given by p → p = ϕ -1 i (p i ), where i ∈ N is the label associated with p ∈ P. If the metric on M is such that the coordinate maps ϕ i themselves have low metric distortion, then the constructed complex Del(P ) is in fact the Delaunay complex of P ⊂ M. This follows from the fact that in the local Euclidean coordinate frames we have ensured that the points have stable Delaunay triangulations. Thus, using Γ 0 = ρ 0 C given by Lemma 5, the stability result [BDG13b, Thm 4.17] leads, by the same reasoning that yields Lemma 18, to the following:

Theorem 2 (Delaunay complex) Suppose that {(W i , ϕ i )} i∈N is an atlas for the compact m-manifold M, and the finite set P ⊂ M is such that Hypotheses 1 and 2 are satisfied. Suppose also that M is equipped with a metric d M , such that

|d i (ϕ i (x), ϕ i (y)) -d M (x, y)| ≤ ηd i (ϕ i (x), ϕ i (y)), whenever x and y belong to ϕ -1 i (B(p i , 6 i )). If η ≤ µ 2 0 2 12 ρ 0 C 2m+1 ,
then Del(P ) is the Delaunay complex of P ⊂ M with respect to d M .

The required bound on η is weaker than the bound on ξ 0 demanded by Hypotheses 2. In the standard scenario, the metric distortion of the transition functions is bounded by bounding the metric distortion of the coordinate charts, and in this case the bound on η required by Theorem 2 is automatically met when Hypotheses 2 are satisfied.

The Riemannian setting

If M is a Riemannian manifold, then, assuming the initial point set P is a (µ 0 , )-net with sufficiently small, Del(P ) is a Delaunay triangulation of M and is equipped with a piecewise flat metric that is a good approximation of d M . This follows from recent results [DVW15] that provide a homeomorphism in this setting. We emphasise that this triangulation result only applies to the perturbed point set P . In general, we cannot expect the Delaunay complex of the initial (µ 0 , )-net P to be homeomorphic to M, regardless of how small may be [BDG13a, App. A].

We use the exponential map to define the coordinate charts. Proposition 17 and Lemma 11 of [DVW15] directly imply that if

ϕ ji = exp -1 ι(j) • exp ι(i) ,
then on B i (p i , r) we have

|d j (ϕ ji (x), ϕ ji (y)) -d i (x, y)| ≤ 6Λr 2 d i (x, y),
where Λ is a bound on the absolute value of the sectional curvatures of M.

Here we will assume a constant sampling radius, i.e., 0 = 0 and j = i for all j ∈ N . For our purposes, we need r = 6 i , and thus ξ 0 = 6 3 Λ 2 i , and in order to satisfy Hypotheses 2 we require 6

3 Λ 2 i ≤ 1 2 4 (Γ 2m+1 0 ) 2 , which is satisfied if i ≤ Γ 2m+1 0 2 6 √ Λ . (10)
We exploit [DVW15, Theorem 3], which guarantees that the output complex is homeomorphic and with a metric quantifiably close to d M . This demands that the star of each p ∈ P ⊂ M be contained in a geodesic ball

B M (p, h) with h = min inj(M) 4 , Γ m 0 6 √ Λ .
Since exp p preserves the radius of a ball centred at p, we have that h = 2 i , and we see that the constraint h ≤

Γ m 0 6 √
Λ is automatically satisfied when i satisfies (10).

We wish to express the required sampling conditions in terms of the intrinsic metric d M . If is the sampling radius with respect to d M , we require an upper bound on such that the needed bound on i is attained when accounting for metric distortion. The Rauch Theorem ([DVW15, Lemma 9]) bounds the metric distortion of the exponential map, and it implies that within a ball of radius r

d i (ϕ i (x), ϕ i (y)) ≤ 1 + Λr 2 3 d M (x, y). (11)
In order to ensure that P i meets the density requirement of item 1 of Hypotheses 1, we demand that B i (p i , 9 i ) ⊆ U i . Then, using (10) to bound r = 9 i , we use (11) to find the bound on required to ensure (10). In fact, the correction is so small that it is already accommodated by the adjustment made when we rounded the constant in (10) to a power of 2. In other words, the right hand side of (10) is already sufficient as a bound on . We also need to ensure that the conditions of Hypotheses 2 are met. In particular, if P is a (µ 0 , )-net with respect to d M , then the effective separation parameter with respect to d i will be slightly smaller, due to the metric distortion of the coordinate charts. In order to compute this correction, we again use the Rauch Theorem [DVW15, Lemma 9], and we find, for p, q ∈ P

d i (ϕ i (p), ϕ i (q)) ≥ 1 - Λr 2 2 µ 0 ≥ 1 -Λr 2 2 1 + Λr 2 3 µ 0 i ,
where r = 9 i , as above. Using (10), and the constraint on Γ 0 imposed by Lemma 5, we find that the correction is indeed extremely small:

d i (ϕ i (p), ϕ i (q)) ≥ (1 -2 -200 )µ 0 i .
We make crude adjustments to the constraint on ρ 0 and the constant defining Γ 0 to accommodate this. We can now formulate [DVW15, Theorem 3] in our context:

Theorem 3 (Riemannian Delaunay triangulation) Suppose M is a Riemannian manifold, and P ⊂ M is a (µ 0 , )-net with respect to the metric d M , with

≤ min inj(M) 4 , 1 2 6 √ Λ ρ 0 C 2m+1 ,
where Λ is a bound on the absolute value of the sectional curvatures and inj(M) is the injectivity radius, and

C = m 3 2 2 µ 0 5m 2 +5m+21
.

If the coordinate charts are defined by

ϕ i = exp -1 ι(i) : B M (ι(i), 10 ) → U i ,
and ρ 0 ≤ µ 0 5 , then the output Del(P ) of the extended algorithm is a Delaunay triangulation: there is a natural homeomorphism H :

|Del(P )| → M that satisfies |d M (H(x), H(y)) -d PL (x, y)| ≤   2 8 Λ C ρ 0 2m 2   d PL (x, y),
where d PL is the natural piecewise flat metric on Del(P ) defined by geodesic distances between vertices in M. In addition, Del(P ) is self-Delaunay: it is a Delaunay triangulation of its vertices with respect to its intrinsic metric d PL .

The general smooth case. The homeomorphism result in the Riemannian setting can be exploited whenever the transition functions are smooth (or at least C 3 ), even if there is no explicit Riemannian metric associated with the input. The reason is that we can construct a Riemannian metric on the manifold by the standard trick using a partition of unity subordinate to the atlas [Boo86, Thm. V.4.5]. Then a short exercise shows that the metric distortion of the coordinate charts is bounded by ξ 0 . It follows then that the constructed Riemannian metric satisfies Theorem 2 if Hypotheses 1 and 2 are satisfied. Thus we can guarantee the existence of a Delaunay triangulation with respect to any smooth metric that can be locally approximated by a Euclidean metric to any desired accuracy (i.e., with arbitrarily small metric distortion): At a sufficiently high sampling density Del(P ) will satisfy Theorem 2 with respect to the given metric, as well as both of Theorems 2 and 3 for the constructed Riemannian metric. The primary example of such a smooth, non-Riemannian metric is the metric of the ambient space R N restricted to a submanifold M ⊂ R N . The associated Delaunay complex is often called the restricted Delaunay complex. Sampling conditions that ensure that the Delaunay complexes associated with the restricted ambient metric and the induced Riemannian metric coincide are worked out in detail from the extrinsic point of view in [START_REF] Boissonnat | Constructing intrinsic Delaunay triangulations of submanifolds[END_REF].

Details of the analysis

In this section we provide details to support the argument made in Section 3.3.

Simplex distortion

Our transition functions introduce a metric distortion when we move from one coordinate chart to another. The geometric properties of a simplex will be slightly different if we consider it with respect to the Euclidean metric d i than they would be if we are using a different Euclidean metric d j . We wish to bound the magnitude of the change of such properties as the thickness and the circumradius of a simplex that is subjected to such a distortion. This is an exercise in linear algebra.

We wish to compare two Euclidean simplices with corresponding vertices, but whose corresponding edge lengths differ by a relatively small amount. The embedding of the simplex in Euclidean space (i.e., the coordinates of the vertices) is not relevant to us. Previous results often only consider the case where the vertices of a given simplex are perturbed a small amount to obtain a new simplex. Lemma 10 demonstrates the existence of an isometry that allows us to also consider the general situation in terms of vertex displacements.

We will exploit observations on the linear algebra of simplices developed in previous work [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]. A k-simplex σ = {p 0 , . . . , p k } in R m can be represented by an m × k matrix P , whose i th column is p ip 0 . We let s i (A) denote the i th singular value of a matrix A, and observe that P = s 1 (P ) ≤ √ k∆(σ).

We are particularly interested in bounds on the smallest singular value of P , which is the inverse of the largest singular value of the pseudo-inverse P † = (P T P ) -1 P T . If the columns of P are viewed as a basis for aff(σ), then the rows of P † may be viewed as the dual basis. The magnitude of a dual vector is equal to the inverse of the corresponding altitudes in σ, and this leads directly to the desired bound on the smallest singular value of P , which is expressed in the following Lemma [BDG13b, Lemma 2.4]:

Lemma 6 (Thickness and singular value) Let σ = [p 0 , . . . , p k ] be a nondegenerate k-simplex in R m , with k > 0, and let P be the m × k matrix whose i th column is p ip 0 . Then the i th row of P † is given by w T i , where w i is orthogonal to aff(σ p i ), and

w i = D(p i , σ) -1 .
We have the following bound on the smallest singular value of P :

s k (P ) ≥ √ kΥ(σ)∆(σ).
We will also have use for a lower bound on the thickness of σ, given the smallest singular value for the representative matrix P . We observe that P was constructed by arbitrarily choosing one vertex, p 0 , to serve as the origin. If there is a vertex p i , different from p 0 , such that D(p i , σ) is minimal amongst all the altitudes of σ, then according to Lemma 6, w i = (kΥ(σ)∆(σ)) -1 , and it follows then that s 1 (P † ) ≥ (kΥ(σ)∆(σ)) -1 , and therefore

s k (P ) ≤ kΥ(σ)∆(σ), (12)
in this case.

We are going to be interested here in purely intrinsic properties of simplices in R m ; properties that are not dependent on the choice of embedding in R m . In this context it is convenient to make use of the Gram matrix P T P , because if Q T Q = P T P , then there is an orthogonal transformation O such that P = OQ. This assertion becomes evident when considering the singular value decompositions of P and Q. Indeed, the entries of the Gram matrix can be expressed in terms of squared edge lengths, as observed in the proof of the following: Lemma 7 Suppose that σ = {p 0 , . . . , p k } and σ = {p 0 , . . . , pk } are two k-simplices in R m such that

| p i -p j -pi -pj | ≤ ξ 0 ∆(σ),
for all 0 ≤ i < j ≤ k. Let P be the matrix whose i th column is p ip 0 , and define P similarly. Consider the Gram matrices, and let E be the matrix that records their difference:

P T P = P T P + E. If ξ 0 ≤ 2 3 , then the entries of E are bounded by |E ij | ≤ 4ξ 0 ∆(σ) 2 , and in particular E ≤ 4kξ 0 ∆(σ) 2 . ( 13
)
Proof Let v i = p ip 0 , and ṽi = pi -p0 . Expanding scalar products of the form (v jv i ) T (v jv i ), we obtain a bound on the magnitude of the coefficients of E:

|ṽ T i ṽj -v T i v j | ≤ 1 2 | ṽi 2 -v i 2 | + | ṽj 2 -v j 2 | + | ṽj -ṽi 2 -v j -v i 2 | ≤ 3 2 (2 + ξ 0 )ξ 0 ∆(σ) 2 ≤ 4ξ 0 ∆(σ) 2 .
This leads us to a bound on s 1 (E) = E . Indeed, the magnitude of the column vectors of E is bounded by √ k times a bound on the magnitude of their coefficients, and the magnitude of s 1 (E) is bounded by √ k times a bound on the magnitude of the column vectors. We obtain Equation (13).

Lemma 7 enables us to bound the thickness of a distorted simplex: Lemma 8 (Thickness under distortion) Suppose that σ = {p 0 , . . . , p k } and σ = {p 0 , . . . , pk } are two k-simplices in R m such that

| p i -p j -pi -pj | ≤ ξ 0 ∆(σ)
for all 0 ≤ i < j ≤ k. Let P be the matrix whose i th column is p ip 0 , and define P similarly.

If

ξ 0 ≤ ηΥ(σ) 2 then s k ( P ) ≥ (1 -η 2 )s k (P ), and
Υ(σ)∆(σ) ≥ 1 √ k (1 -η 2 )Υ(σ)∆(σ), and
Υ(σ) ≥ 4 5 √ k (1 -η 2 )Υ(σ).

Proof

The equation P T P = P T P + E implies that

|s k ( P ) 2 -s k (P ) 2 | ≤ s 1 (E),
and so

|s k ( P ) -s k (P )| ≤ s 1 (E) s k ( P ) + s k (P ) ≤ s 1 (E) s k (P ) . Thus s k ( P ) ≥ s k (P ) - s 1 (E) s k (P ) = s k (P ) 1 - s 1 (E) s k (P ) 2 .
From Lemma 7 and the bound on ξ 0 we have

s 1 (E) ≤ η 2 kΥ(σ) 2 ∆(σ) 2 ,
and so s 1 (E) s k (P ) 2 ≤ η 2 by Lemma 6, and we obtain s k ( P ) ≥ (1η 2 )s k (P ). For the thickness bound we assume, without loss of generality, that there is some vertex different from p0 that realises the minimal altitude in σ (our choice of ordering of the vertices is unimportant, other than to establish the correspondence between σ and σ). Thus Equation (12) and Lemma 6, give the inequalities kΥ(σ)∆(σ) ≥ s k ( P ), and

s k (P ) ≥ √ kΥ(σ)∆(σ),
and we obtain kΥ

(σ)∆(σ) ≥ (1 -η 2 ) √ kΥ(σ)∆(σ).
The final result follows since ∆(σ) ∆(σ) ≥ 1 1+ξ 0 ≥ 4 5 .

In order to obtain a bound on the circumradius of σ with respect to that of σ, it is convenient to find an isometry that maps the vertices of σ close to the vertices of σ. Choosing p0 and p 0 to coincide at the origin, the displacement error for the remaining vertices is minimised by taking the orthogonal polar factor of the linear transformation A = P P -1 that maps σ to σ. In other words, if the singular value decomposition of A is A = U A Σ A V T A , then A = ΦS, where S = V A Σ A V T A , and Φ = U A V T A is the desired linear isometry. We have the following result, which is a special case of a theorem demonstrated by Jiménez and Petrova [JP13]:

Lemma 9 (Close alignment of bases) Suppose that P and P are nondegenerate k × k matrices such that

P T P = P T P + E. ( 14
)
Then there exists a linear isometry Φ : R k → R k such that

P -ΦP ≤ s 1 (P )s 1 (E) s k (P ) 2 .

Proof

Multiplying by P -T := (P T ) -1 on the left, and by P -1 on the right, we rewrite Equation ( 14) as

A T A = I + F, (15)
where A = P P -1 , and F = P -T EP -1 . Using the singular value decomposition A = U A Σ A V T A , we let Φ = U A V T A , and we find P -ΦP = (A -Φ)P = U A (Σ A -I)V T A P.

From Equation (15) we deduce that s 1 (A) 2 ≤ 1 + s 1 (F ), and also that

s k (A) 2 ≥ 1 -s 1 (F ). It follows that max i |s i (A) -1| ≤ s 1 (F ) 1 + s i (A) ≤ s 1 (F ),
and thus

Σ A -I ≤ s 1 (F ) ≤ s 1 (P -1 ) 2 s 1 (E) = s k (P ) -2 s 1 (E).
The result now follows from Equation ( 16).

Recalling that an upper bound on the norm of a matrix also serves as an upper bound on the norm of its column vectors, we obtain the following immediate consequence of Lemma 9, using Lemma 7 and Lemma 6:

Lemma 10 (Close alignment of simplices) Suppose that σ = {p 0 , . . . , p k } and σ = {p 0 , . . . , pk } are two k-simplices in R m such that

| p i -p j -pi -pj | ≤ ξ 0 ∆(σ),
for all 0 ≤ i < j ≤ k. Let P be the matrix whose i th column is p ip 0 , and define P similarly. If ξ 0 ≤ 2 3 , then there exists an isometry Φ : R m → R m such that

P -ΦP ≤ 4 √ kξ 0 ∆(σ) Υ(σ) 2 ,
and if σ = Φσ = {p 0 , . . . , pk }, then p0 = p0 , and

pi -pi ≤ 4 √ kξ 0 ∆(σ) Υ(σ) 2 for all 1 ≤ i ≤ k.
Using Lemma 10 together with [BDG13b, Lemma 4.3] we obtain a bound on the difference in the circumradii of two simplices whose edge lengths are almost the same:

Lemma 11 (Circumradii under distortion) Suppose that σ = {p 0 , . . . , p k } and σ = {p 0 , . . . , pk } are two k-simplices in R m such that

| p i -p j -pi -pj | ≤ ξ 0 ∆(σ), for all 0 ≤ i < j ≤ k. If ξ 0 ≤ Υ(σ) 4 2 , then |R(σ) -R(σ)| ≤ 16k 3 2 R(σ)ξ 0 Υ(σ) 3 .

Proof

We define σ = Φσ, where Φ : σ → aff(σ) is the isometry described in Lemma 10. Since p0 = p0 , and

R(σ) = R(σ), we have |R(σ) -R(σ)| ≤ C(σ) -C(σ)
. By Lemma 10, and the hypothesized bound on ξ 0 , the distances between C(σ) and the vertices of σ are all bounded by

R(σ) + 4 √ kξ 0 ∆(σ) Υ(σ) 2 ≤ 1 + √ k 2 R(σ) ≤ 3 √ k 2 R(σ),
and these distances differ by no more than

8 √ kξ 0 ∆(σ) Υ(σ) 2 . It follows then from [BDG13b, Lemma 4.3] that C(σ) -C(σ) ≤ 3 √ k 2 R(σ) Υ(σ)∆(σ) 8 √ kξ 0 ∆(σ) Υ(σ) 2 ≤ 12kR(σ)ξ 0 3 4 √ k Υ(σ) 3 by Lemma 8, with η = 1 2 ≤ 16k 3 2 R(σ)ξ 0 Υ(σ) 3 .

Circumcentres and distortion maps

It is convenient to introduce the affine space N (σ), which is the space of centres of circumscribing balls for a simplex σ ∈ R m . If σ is a non-degenerate k-simplex, then N (σ) is an affine space of dimension mk perpendicular to aff(σ) and containing C(σ).

The transition functions introduce a small metric distortion, which motivated our interest in the properties of perturbed simplices. In order to extend the perturbation algorithm [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] to the setting of curved manifolds, we are interested in quantifying how the test for the hoop property behaves under a perturbation of the interpoint distances. Specifically, if a point p is at a distance α 0 R from the diametric sphere of a simplex σ in one coordinate frame, what can we say about the distance of p from S m-1 (σ) when measured by the metric of another coordinate frame? To this end, we are interested in the behaviour of the circumcentre under the influence of a mapping that is not distance preserving. As a first step in this direction, we observe another consequence of [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]Lemma 4. Proof By the perturbation bounds on φ, the distances between φ(c) and the vertices of σ differ by no more than 2ξ 0 r, and these distances are all bounded by (1 + ξ 0 )r. In this context [BDG13b, Lemma 4.3] says that there exists a c ∈ N (σ) such that d(φ(c), c) ≤

(1 + ξ 0 )r2ξ 0 r Υ(σ)∆(σ) .

We apply Lemma 8, using η = 1 2 , to obtain Υ(σ)∆(σ) ≥ 3 and similarly for the upper bound.

We will find it convenient to have a bound on the circumradius of a simplex, relative to its thickness and longest edge length:

Υ 3 + 6k 1 2 R Υ∆ + 3 • 32k 2 ξ 0 R Υ 4 ∆ + (2 + ξ 0 ) ξ 0 R 2 < 16k 3 Υ 4 + 32k 3 2 Υ 3 + 3k 1 2
Υ 2 + 3k 2 Υ 3 + 3 ξ 0 R 2 using Lemma 13

≤ 57k 3 Υ 4 ξ 0 R 2 .

The size of the domains

The domains U ij on which the transition functions are defined need to be large enough to accommodate two distinct requirements. First, the domain of the transition function ϕ ji must contain a large enough neighbourhood

Proof

Since P i is a (µ 0 , i )-net for B i (p i , 8 i ), it follows that Q i is a (µ 0 , i )net. Thus if there are no forbidden configurations in Q i , then by [BDG14, Lemma 3.6], all the m-simplices in Del | (Q i ) will be Γ 0 -good and δ-protected, with δ = δ 0 µ 0 i , where Del | (Q i ) is the subcomplex of Del(Q i ) consisting of simplices that have an empty circumscribing ball centred in D i (Q i ) (as defined in Equation ( 1)).

The sampling criteria ensure that every point on ∂conv(Q i ) must be at a distance of less than 2 i from ∂B i (p i , 6 i ). According to Lemma 2 P3, if τ is a forbidden configuration in Q i , then ∆(τ ) < 15 4 i , and it follows that if p j ∈ τ , then τ ⊂ B i (p j , 4 i ). We will require that each potential forbidden configuration in Q i lies within the domain of any transition function associated with one of its vertices. Thus we demand that

B i (p j , 4 i ) ∩ B i (p i , 6 i ) ⊂ U ij if p j ∈ B i (p i , 6 i ). (20)
For simplicity we accommodate Equations ( 19) and (20) by demanding that

B i (p j , 9 i ) ∩ B i (p i , 6 i ) ⊂ U ij if p j ∈ B i (p i , 6 i ). (21)
In summary, Lemmas 15, 16, and 17 combine to yield the desired equivalence of stars (18), under the assumption that Q i has no forbidden configurations. We take U = B i (p i , 6 i ) in Lemma 15, and Equation (4) yields ξ ≤ ξ 0 12 i . Using δ = Γ m+1 0 µ 0 i (see Lemma 5) and (2), we require ξ 0 ≤ Γ 2m+1 0 µ 2 0 /(5 • 12 • 36), so we obtain:

Lemma 18 (Stable stars) If

ξ 0 ≤ Γ 2m+1 0 µ 2 0 2 12 ,
and there are no forbidden configurations in Q i , then for all p j ∈ star(p i ; Del(P i )), we have star(p i ; Del(P i )) ∼ = star(p i ; Del(P j )).

We have established minimal requirements on the size of the domains U ij , but these requirements may implicitly demand more. Although ϕ ji : U ij → U ji is close to an isometry, i may be almost twice as large as j . Thus the

,

  3]:Lemma 12 (Circumscribing balls under distortion) Suppose φ : R m ⊃ U → V ⊂ R m is a homeomorphism such that, for some positive ξ 0 ,|d(x, y)d(φ(x), φ(y))| ≤ ξ 0 d(x, y)for all x, y ∈ U.Suppose also that σ ⊂ U is a k-simplex, and that B(c, r) is a circumscribing ball for σ with c ∈ U . Let σ = φ(σ). If then there is a circumscribing ball B(c, r) for σ such thatd(φ(c), c) ≤ 3 √ kr 2 ξ 0 Υ(σ)∆(σ) ,and|r -r| ≤ 5 √ kr 2 ξ 0 Υ(σ)∆(σ).

4 √k

 4 Υ(σ)∆(σ). We findd(φ(c), c) ≤ 8 √ k(1 + ξ 0 )r 2 ξ 0 3Υ(σ)∆(σ) .The announced bound on d(φ(c), c) is obtained by observing that ξ 0 ≤ 1 16 . Choosing a vertex p = φ(p) ∈ σ, the bound on the difference in the radii follows:r = d(p, c) ≥ d(p, φ(c))d(φ(c), c) ≥ rξ 0 r -3 √ kr 2 ξ 0 Υ(σ)∆(σ) ≥ r -5 √ kr 2 ξ 0 Υ(σ)∆(σ),

Lemma 13 Figure 2 :

 132 Figure 2: Diagram for the proof of Lemma 14.

  Thus d i (p i , ∂conv(Q i )) > 6 i -2 i ≥ 14 5 i . Also, d i (p i , p i ) ≤ i 4, and we find thatd i (p i , ∂conv(Q i )) ≥ 51 20 i . Thus, since d i (p i , C(σ)) < i if σ is in star(p i ; Del(Q i )), we have star(p i ; Del(Q i )) ⊆ Del | (Q i ), and hence the result.

A manifold simplicial complex that admits an atlas of piecewise linear coordinate charts from the stars of the vertices is called piecewise linear . There exists manifold simplical complexes that are not piecewise linear [Thu97, Example 3.2.11].

with η 2 ≤ 1,

Acknowledgements

This research has been partially supported by the 7th Framework Programme for Research of the European Commission, under FET-Open grant number 255827 (CGL Computational Geometry Learning). Partial support has also been provided by the Advanced Grant of the European Research Council GUDHI (Geometric Understanding in Higher Dimensions).

Arijit Ghosh is supported by Ramanujan Fellowship number SB/S2/RJN-064/2015. Part of this work was done when he was a researcher at the Max-Planck-Institute for Informatics, Germany supported by the IndoGerman Max Planck Center for Computer Science (IMPECS). Part of this work was also done while he was a visiting scientist at the Advanced Computing and Microelectronics Unit, Indian Statistical Institute, Kolkata, India.

We gratefully profited from discussions with Mathijs Wintraeken. We also thank the reviewers for the careful reading and thoughtful comments that significantly improved the final manuscript.

Proof

Let σ = {p 0 , . . . , p k }, We work in R k ∼ = aff(σ) ⊂ R m , and let P be the k × k matrix whose i th column is p ip 0 . Then, by equating C(σ)p 0 2 with C(σ)p i 2 and expanding, we find a system of equations that may be written in matrix form as

where the i th component of the vector b is 1 2 ( p i 2p 0 2 ). Choosing p 0 as the origin, we have C(σ) = R(σ), and b ≤ 1 2 √ k∆(σ) 2 . From (17) and Lemma 6 we also have

and the result follows.

Using the bound on d(φ(C(σ)), N (σ)) given by Lemma 12, together with the circumradius bound of Lemma 11, we obtain a bound on d(φ(C(σ)), C(σ)) by means of the Pythaogrean theorem: Lemma 14 (Circumcentres under distortion) Suppose φ : R m ⊃ U → V ⊂ R m is a homeomorphism such that |d(x, y)d(φ(x), φ(y))| ≤ ξ 0 d(x, y) for all x, y ∈ U.

Suppose also that σ ⊂ U is a k-simplex with C(σ) ∈ U , and let σ = φ(σ). If

Proof

Let c be the closest point in N (σ) to φ(C(σ)), and let w be the distance from c to φ(C(σ)). Setting z as the distance between c and C(σ), we have that d(φ(C(σ)), C(σ)) 2 = z 2 + w 2 ; see Figure 2. Let ĉ be the orthogonal projection of φ(C(σ)) into aff(σ). Then, letting R = R(σ), and R = R(σ), and choosing p = φ(p) ∈ σ, we have

of p i that we can apply the metric stability result of [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] to ensure that star(p i ; Del(P i )) will be the same as star(p i ; Del(P j )) whenever p j ∈ star(p i ; Del(P i )). The second requirement is that any potential forbidden configuration in the region of interest must lie entirely within the domain of the transition function associated with each of its vertices. We need some more terminology and notation in order to recall the stability result that we will use. The (µ 0 , i )-net Q i ⊂ U ij will be introduced shortly. A point q ∈ Q i is an interior point if it is not on the boundary of conv(Q i ). For an abstract simplicial complex K whose vertices are Q i , and S ⊂ Q i , we denote by star(S; K) the subcomplex of K consisting of simplices that include a point in S, together with the faces of these simplices. If d is a metric on U ij , then Del d (Q i ) is the Delaunay complex computed with respect to the metric d. Since d i is the standard Euclidean metric on

In this notation, the stability result [BDG13b, Theorem 4.17] reads as follows.

Lemma 15 (Delaunay stability under metric perturbation) Suppose

For our purposes, d j is the pullback by ϕ ji of the Euclidean metric on U j . Thus we have the identification star(S; Del

We will use S = {p i }. Some argument is required to ensure that Lemma 15 provides a route to the desired equivalence star(p i ; Del(P i )) ∼ = star(p i ; Del(P j )), when p j ∈ star(p i ; Del(P i )). ( 18)

We first establish our "region of interest". We demand, for all i ∈ N , that P i be a (µ 0 , i )-net for B i (p i , 8 i ), and we define Q i = P i ∩ B i (p i , 6 i ). Since P i changes as the algorithm progresses, points may come and go from Q i , but we will ensure that when the algorithm terminates, Q i will contain no forbidden configurations.

Lemma 16 For all i ∈ N we have

and if p ∈ star(p i ; Del(P i )), then

Proof

The density assumption guarantees that if

, establishing the first equality, and since 3 i > 2 i , we also obtain the bound on the distance from p to ∂B i (p i , 6 i ).

The second equality follows from two observations. First we show that if σ m ∈ star(p i ; Del(ϕ ji (Q i ))), then R(σ) < j . Since ϕ ji (Q i ) ⊂ P j , and P j is j -dense for B = B j (p j , 8 j ), it is sufficient to show that d j (p i , ∂B) ≥ 2 j . Since p j ∈ star(p i ; Del(P i )), we have

Thus since d j (p j , p j ) ≤ 1 4 j , we have d j (p j , ∂B) ≥ 8 j -21 4 j = 11 4 j ≥ 11 5 j . This establishes that the Delaunay ball for σ m must remain empty when points outside of B are considered.

The second observation required to establish the second equality is that if q ∈ P j is such that d j (p i , q ) < 2 j , then q ∈ ϕ ji (B i (p i , 6 i )). Indeed, we have

If p j ∈ star(p i ; Del(P i )), then d i (p i , p j ) < 1 2 i + 2 i ≤ 3 i . Thus Lemma 16 establishes the first requirement on U ij , namely

The second requirement arises from the fact that we wish to ensure that there are no forbidden configurations in Q i . This will be sufficient for us to apply Lemma 15.

Lemma 17 (Protected stars) If there are no forbidden configurations in Q i , then all the m-simplices in star(p i ; Del(Q i )) are Γ 0 -good and δ-protected, with δ = δ 0 µ 0 i .

requirement on U ij may imply that U ji = ϕ ji (U ij ) is significantly larger than Equation ( 21) demands. Clearly we must have

We have also explicitly demanded that P i be a (µ 0 , i )-net for B i (p i , 8 i ). We will assume that B i (p i , 8 i ) ⊂ U i .

Hoop distortion

We will rely primarily on Properties P1 and P4 of forbidden configurations (Lemma 2), and the stability of the circumcentres exhibited by Lemma 14.

We have the following observation about the properties of forbidden configurations under the influence of the transition functions:

has a radius satisfying

Proof

The bound for Γk 0 follows immediately from Lemma 8, and the fact that σ is Γ k 0 -thick (Lemma 2 P4). Likewise, the radius bound is a direct consequence of Lemma 11 and Lemma 2 P2.

The bound on α0 is obtained from Property P1 with the aid of Lemmas 11 and 14. We have d i (p i , S m-1 (σ)) = |d i (p i , C(σ)) -R(σ)|, and we are able to get a tighter upper bound on R(σ)d i (p i , C(σ)), than we can for

We have abused the notation slightly because τ = ϕ ij (τ ) need not actually satisfy the α0 -hoop property definition d i (p, S(τ p )) ≤ α0 R(τ p ), because R(τ ) may be smaller than 2 i . However we are not concerned with the α0 -hoop property for τ ; instead we desire a condition that will permit the extended algorithm to emulate the original Euclidean perturbation algorithm [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF], and guarantee that forbidden configurations such as τ cannot exist in any of the sets Q j .

The bounds in Lemma 19 can be further simplified. We have announced them in this intermediate state in order to elucidate the roles played by ξ 0 and 0 . In particular, there is no need to significantly constrain 0 . The original perturbation algorithm for points in Euclidean space [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] extends to the case of a non-constant sampling radius simply by replacing α 0 by α0 ≤ (1 + 0 )α 0 ≤ 2α 0 , as can be seen by setting ξ 0 = 0 in the expression for α0 . This shows the utility of 0 and the explicit local sampling radii: such a variable sampling radius could not be constructed just by varying the local coordinate charts due to our requirement that ξ 0 be extremely small.

In the general case of interest here, we see from the expression for α0 presented in Lemma 19, that ξ 0 must be considerably more constrained with respect to Γ 0 if we are to obtain an expression for α0 that goes to zero as Γ 0 goes to zero. For the purposes of the algorithm, we do not require the bounds on the radius or the thickness.