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

Siargey Kachanovich
email: siargey.kachanovich@inria.fr

Mathijs Wintraecken

Triangulating Submanifolds: An Elementary and Quantified Version of Whitney's Method

Keywords: Triangulations, Manifolds, Coxeter triangulations 1 Introduction

We quantise Whitney's construction to prove the existence of a triangulation for any C 2 manifold, so that we get an algorithm with explicit bounds. We also give a new elementary proof, which is completely geometric.

Theorem 1.1 Every compact n-dimensional C 2 manifold M embedded in R d admits a triangulation.

We note that C 2 -manifolds have positive reach, see [START_REF] Federer | Curvature measures[END_REF]. The reach rch M was introduced by Federer [START_REF] Federer | Curvature measures[END_REF], as the minimal distance between a set M (in this paper always a manifold) and its medial axis.

By more quantitative, we mean that instead of being satisfied with the existence of constants that are used in the construction, we want to provide explicit bounds in terms of the reach of the manifold, which we shall assume to be positive. The medial axis consists of points in ambient space that do not have a unique closest point on M . Federer [START_REF] Federer | Curvature measures[END_REF]Rem. 4.20] also mentions that manifolds are of positive reach if and only if they are C 1,1 . It is not too difficult to generalise the precise quantities to the setting where the manifold is C 1,1 (instead of C 2 ) at a small cost, see Appendix C.

Note that Theorem 1.1 implies that any C 1 manifold admits a triangulation. This is because any C 1 manifold can be smoothed (see for example [START_REF] Hirsch | Differential Topology[END_REF]) and Whitney's own embedding theorem [START_REF] Whitney | Geometric Integration Theory[END_REF]Sect. IV.1] gives a smooth embedding in R d .

Triangulations in computational geometry and topology are most often based on Voronoi diagrams and their dual Delaunay triangulations of the input point set, see for example [START_REF] Berger | State of the art in surface reconstruction from point clouds[END_REF][START_REF] Boissonnat | Meshing of surfaces[END_REF][START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF][START_REF] Cheng | Delaunay Mesh Generation[END_REF][START_REF] Dey | Curve and Surface Reconstruction[END_REF] for general references in low dimensions and more recent work on manifolds embedded in higher dimensional spaces [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF][START_REF] Cheng | Manifold reconstruction from point samples[END_REF]. Whitney's construction is of a quite different nature. He uses an ambient triangulation and constructs the triangulation of the manifold M based on the intersections of M with this triangulation. In this paper, we have chosen this ambient triangulation T to be (a perturbation of) a Coxeter triangulation T of type Ãd . A Coxeter triangulation of type Ãd is Delaunay protected, a concept we will recall in detail in Sect. [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF]. Delaunay protection gives that the triangulation is stable under perturbations. This property simplifies the proof, which in fact was one of the motivations for our choice. Moreover, Coxeter triangulations can be stored very compactly, in contrast with previous work [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF][START_REF] Cheng | Manifold reconstruction from point samples[END_REF] on Delaunay triangulations.

The approach of the proof of correctness of the method, that we present in this paper, focuses on proving that after perturbing the ambient triangulation the intersection of each d-simplex in the triangulation T with M is a slightly deformed n-dimensional convex polytope, more precisely the intersection is piecewise smoothly homeomorphic to a polytope. Proving this is the core of the homeomorphism proof in Sect. 7. The triangulation K of M consists of a barycentric subdivision of a straightened version of these polytopes. This may remind the reader of the general result on CW-complexes, see [START_REF] Lundell | The Topology of CW Complexes[END_REF], which was exploited by Edelsbrunner and Shah [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF] for their triangulation result.

In this paper we construct 'normals' and a tubular neighbourhood for K that is compatible with the ambient triangulation T and prove that the projection along these 'normals' is a homeomorphism. This interpretation of Whitney's triangulation method is different from Whitney's original proof where the homeomorphism is given by the closest point projection and uses techniques which we also exploited in [START_REF] Boissonnat | Local criteria for triangulation of manifolds[END_REF]. The homeomorphism we give in this paper is in fact piecewise smooth. We stress that this result is stronger than if we had based our work on the closed ball property of Edelsbrunner and Shah, with given criteria for a homeomorphism, but not for a 123 piecewise linear/smooth homeomorphism nor an explicit map. We also believe that the tubular neighbourhood we construct is of independent interest. Because we have a bound on the size of the tubular neighbourhood of K and M lies in this neighbourhood, we automatically bound the Hausdorff distance between the two. A bound on the difference between the normals of K and M is also provided. Thanks to our choice of ambient triangulation and our homeomorphism proof, this entire paper is elementary in the sense that no topological results are needed, all arguments are geometrical.

In addition to the more quantitative/algorithmic approach, the purely geometrical homeomorphism proof, the link with the closed ball property, the tubular neighbourhood for the triangulation K , and a bound on the Hausdorff distance, we also give different proofs for a fair number of Whitney's intermediate results.

In spite of this paper not being a review, the authors hope that it will serve to spread awareness of the classical work by Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF] in the computational geometry and applied math communities. The main reason for this is that a large number of authors has reintroduced (weaker) versions of Whitney's concepts and results, without having been aware of the original.

The marching cube algorithm and some of its variants [START_REF] Allgower | An algorithm for piecewise-linear approximation of an implicitly defined manifold[END_REF][START_REF] Doi | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[END_REF][START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] provide ways to approximate a manifold that is the zeroset of a function. We will call such a manifold an isomanifold. These algorithms use a subdivision of the ambient space into simplices or cubes and constructing a piecewise linear approximation of the isomanifold inside each simplex or cube. This coincides with Whitney's approach where he subdivides the ambient space into cubes, which he then subdivides into simplices and then approximates the manifold inside each simplex. The main difference is that Whitney needs a perturbation of the ambient triangulation to guarantee topological correctness, while (with the exception of [START_REF] Boissonnat | Isotopic implicit surface meshing[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] in two and three dimensions) no topological correctness (homeomorphism) is proved for the marching cube algorithms. Whitney is also more general because he treats general manifolds and not just isomanifolds. Moreover, Allgower and Georg [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF]Thm. 15.4.1] assume that the isomanifold avoids simplices in the ambient triangulation whose dimension is strictly less than the codimension of the isomanifold to prove that the piecewise linear approximation of the manifold is itself a manifold. This idea also originates from Whitney, and will be discussed in detail below.

Whitney's idea of using a subdivision of ambient space as a scaffold to build a triangulation has also been adopted outside of the marching cube community, see for example [START_REF] Cheng | Dynamic skin triangulation[END_REF]. In [START_REF] Cheng | Dynamic skin triangulation[END_REF] the scaffolding is based on the Voronoi diagram of a point sample. This is unlike the ambient triangulation used by Whitney. The focus on three dimensional ambient space and a specific type of surface, instead of general manifolds of arbitrary dimension and codimension, further distinguishes it from Whitney's work. As mentioned above, the idea to use barycentric subdivision to construct a triangulation has also been often used, e.g. in [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF][START_REF] Lundell | The Topology of CW Complexes[END_REF].

The part of the algorithm described in this paper that constructs the triangulation (see part 2 of the algorithm in Sect. 2.1) and the data structure to store the ambient triangulation have been implemented, see [START_REF] Boissonnat | Sampling and meshing submanifolds in high dimension[END_REF] and [START_REF] Kachanovich | Meshing Submanifolds Using Coxeter Triangulations[END_REF]. The implementation of the perturbation scheme (see part 1 of the algorithm in Sect. 2.1) is not yet complete at the moment of writing.

The Algorithm (Based on Whitney's Construction)

The algorithm takes as input an n-dimensional C 2 manifold M ⊂ R d with reach rch M , and outputs the triangulation K of M . The algorithm based on Whitney's construction consists of two parts: We will refer to the first part as the perturbation algorithm. The perturbation algorithm perturbs the vertices of the ambient triangulation which ensures that the intersection of the ambient simplices with the manifold is nice (the intersection is piecewise smoothly homeomorphic to a polytope as we will prove in Sect. 7). The second part is where the triangulation is constructed and is based on barycentric subdivision of polytopes.

Part 1 (the perturbation algorithm): This part of the algorithm outputs a perturbed version of a Coxeter triangulation of R d of type Ãd (see Sect. 4 for the precise definition) and consists of two steps. In these two steps we have to carefully choose a significant 123 number of parameters, which we will not discuss in detail in the global description of the algorithm. An overview of the most important parameters and notation can be found in Appendix A.

-Choose a Coxeter triangulation T of type Ãd of R d that is sufficiently fine. Here by fine we mean as determined by the longest edge length L. The longest edge length L is linear in the reach and depends in a rather intricate manner on the thickness (minimal altitude over longest edge length) of the top dimensional simplices in Ãd and the dimension and codimension of the manifold. The precise expression will be given in [START_REF] Boissonnat | Meshing of surfaces[END_REF]. -Perturb the vertices of T slightly into a T (with the same combinatorial structure), such that all simplices in T of dimension at most dn -1 are sufficiently far away from the manifold. Here slightly is in terms of quality, protection (see Sect. 4 for the definitions) of the ambient triangulation as well as the longest edge length, separation, and dimension; see [START_REF] Boissonnat | Sampling and meshing submanifolds in high dimension[END_REF] for the precise bounds. Sufficiently far means some small fraction of the longest edge length and thus even smaller fraction of the reach of M . The precise bound can be found in [START_REF] Boissonnat | Delaunay stability via perturbations[END_REF]. This is done as follows:

One maintains a list Ti of vertices and simplices, starting with an empty list and adding perturbed vertices while keeping the combinatorial structure of T intact. This means that if τ = {v j 1 , . . . , v j k } is a simplex in T and ṽ j 1 , . . . , ṽ j k ∈ Ti , where ṽi denotes the perturbed vertex v i , then τ = {ṽ j 1 , . . . , ṽ j k } is a simplex in Ti . We shall think of Ti simultaneously as a list, a simplicial complex, and a triangulation of a subset of R d . We shall think of i as the index of the vertex that was added last. To this list Ti , one first adds all vertices v i of T such that d(v i , M ) ≥ 3L/2, as well as the simplices with these vertices (see Case 1 of Sect. 5.2). For a vertex v i such that d(v i , M ) < 3L/2 (Case 2), one goes through the following procedure. We first pick a point p ∈ M that is not too far from v i .

We then consider all τ j ⊂ Ti-1 of dimension at most dn -2, such that the join1 v i * τ j lies in Ti . For all such τ j we consider span(τ j , T p M ) and we pick our perturbed v i , that is ṽi , so that it lies sufficiently far from the union of these spans, but also not too far from v i (as we mentioned at the beginning). Here sufficiently far means a very small fraction of the longest edge length, see [START_REF] Cairns | On the triangulation of regular loci[END_REF]. The existence of such a point can by proved by volume estimates and is shown in Lemma 5.6.

The fact that such a perturbation ensures that the (d -n -1)-skeleton is sufficiently far away from the manifold is non-trivial and is proved in Lemma 5.7.

We note that for a curve in two dimensions, as depicted in Fig. 1, or more generally a manifold of codimension 1, the set of all τ j ⊂ Ti-1 of dimension at most dn -2 is the empty set and span(τ j , T p M ) is T p M . The perturbation therefore ensures that ṽi lies far from T p M .

Note that we only require limited knowledge of the manifold. Given a vertex v i we need to be able to find a point on M that is close to v i or know if v i is far from M and we need access to T M in a finite sufficiently dense set of points (so that for every point v i that is close to M we have a linear approximation of M ). We assume we have two oracles for the two operations. There are no fundamental difficulties in including small uncertainties in our knowledge of the close points or the tangent spaces, but the analysis would be more complicated. If we can sample M densely finding close points is algorithmically not difficult. Methods to estimate the tangent space have been described in [START_REF] Aamari | Nonasymptotic rates for manifold, tangent space and curvature estimation[END_REF]. The same paper also describes estimates on the curvature. The estimate of the reach is discussed in [START_REF] Dey | Normal and feature approximations from noisy point clouds[END_REF] in three dimensions and [START_REF] Aamari | Estimating the reach of a manifold[END_REF] in high dimensions.

Complexity of part 1: The complexity of the perturbation (per vertex) of the algorithm is dominated by the number of simplices τ j that we have to consider. This number is bounded by the number of simplices of dimension at most dn -2 in the star of a vertex in a Coxeter triangulation plus 1, see (4) below. The number of simplices in turn is bounded by Lemma 4.11. This compares favourably with the complexity of the perturbation method in [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] for Delaunay triangulations, which is of order O(2 d 2 ). A full analysis of the complexity of the algorithm, including basic operations on Coxeter triangulations, will be reported upon in a separate paper.

(d -n) d d d-n , see
Part 2 (the triangulation construction): The construction of the triangulation of M is now straightforward barycentric subdivision; for each τ k ∈ T , of dimension k, that contains a part of M , we pick a point v(τ k ) in τ k , see [START_REF] Coxeter | Discrete groups generated by reflections[END_REF]. For any sequence

τ d-n ⊂ τ d-n+1 ⊂ • • • ⊂ τ d ,
such that all simplices in the sequence intersect M we add a simplex {v(τ d-n ), . . . , v(τ d )} to a simplicial complex K . If we have done this for all simplices that contain M , K is a triangulation of M . For this second part we need an oracle that is able to tell us if the intersection between M and τ d-n ∈ T is non-empty and if so, gives us the point of intersection. As we will see in Sect. 6.1, it would in fact suffice to be able to find intersections between tangent planes and simplices.

A Nice Byproduct

The triangulation algorithm does not only provide a triangulation of the manifold itself, with simplices whose quality is lower bounded. It in fact immediately gives that the barycentric subdivision of the ambient triangulation contains a triangulation of the manifold. To ensure that the triangulation of the manifold is geometrically close to the manifold, we need to shift (some of the) vertices to the position that is computed by the algorithms above, see Fig. 2. Because the simplices of the triangulation of the manifold have good quality, we find a triangulation of the ambient space whose simplices have good quality. This byproduct may be of interest for finding numerical solutions to partial differential equations, in particular for space time methods [START_REF] Behr | Simplex space-time meshes in finite element simulations[END_REF][START_REF] Von Danwitz | Simplex space-time meshes in compressible flow simulations[END_REF][START_REF] Shewchuk | What is a good linear finite element? Interpolation, conditioning, anisotropy, and quality measures[END_REF]. This also serves as a first step in generalising the work on the triangulation of general stratifolds in three dimensions [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF][START_REF] Dey | Localized Delaunay refinement for volumes[END_REF][START_REF] Dey | Localized Delaunay refinement for piecewise-smooth complexes[END_REF][START_REF] Oudot | Meshing volumes bounded by smooth surfaces[END_REF][START_REF] Rineau | Meshing Volumes Bounded by Piecewise Smooth Surfaces[END_REF], which may be of interest given the effort that went into the detection of strata in arbitrary dimension, see for example [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF][START_REF] Bendich | Stratification learning through homology inference[END_REF][START_REF] Brown | Sheaf-theoretic stratification learning[END_REF].

123 Fig. 2 The same triangulation as depicted in Fig. 1 with the addition of simplices of the barycentric subdivision of the ambient triangulation added for the simplices that intersect the manifold

Outline and Overview of the Proof

This paper is dedicated to the correctness proof of the algorithm presented in Sect. 2.1. After some background sections dedicated to manifolds of positive reach and Coxeter triangulations and their stability under perturbations, we continue with the perturbation algorithm.

In Sect. 3 we recall some results on the geometry of manifolds of positive reach. Coxeter triangulations, Delaunay protection, and the combinatorial stability of a triangulation under perturbations is the topic of Sect. 4.

In Sect. 5, we both give the details of the perturbation of the vertices and some geometric consequences for the triangulation. In Sect. 6, the triangulation K of M is defined and an important quality bound for the simplices is given. Section 7 is dedicated to proving that K is a triangulation of M . The proof is quite different from the approach Whitney described, which uses the closest point projection onto M . Here we construct a tubular neighbourhood and 'normals' around the triangulation K , which is adapted to the ambient triangulation T . We then prove that the projection using these 'normals' gives a piecewise smooth homeomorphism from

τ d ∩ M to τ d ∩ K , where τ d ∈ T is d-dimensional.
Because the construction is compatible on the faces of d-dimensional simplices, the global result immediately follows. A more detailed overview of the homeomorphism proof is given in Sect. 7.

Manifolds, Tangent Spaces, Distances, and Angles

In this section, we discuss some general results that will be of use. The manifold M ⊂ R d is a compact C 2 manifold with reach rch M . We adhere as much as possible to the same notation as used in [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF]. The tangent bundle will be denoted by T M , while the tangent space at a point p is written as T p M . Similarly, N M is the normal bundle and N p M the normal space. Distances on the Fig. 3 The cylinder C (T p M , r 1 , r 2 ), with the manifold and tangent space manifold will be indicated by d M ( • , • ), while we write d( • , • ) for distances in the ambient Euclidean space, and | • | for the length of vectors. A ball centred at x with radius r is denoted by B(x, r ). For a point x in the ambient space such that d(x, M ) < rch M , the closest point projection onto M is denoted by π M (x). The orthogonal projection onto the tangent T p M is denoted by π T p M (x).

We will use a result from [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF], which improves upon previous works such as Niyogi et al. [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]: Lemma 3.1 [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF]Lem. 6 and Corr. 3]) Suppose that M is C 2 and let p, q ∈ M , then

∠(T p M , T q M ) ≤ d M ( p, q) rch M and sin ∠(T p M , T q M ) 2 ≤ | p -q| 2 rch M .
In Lemma 3.2 we prove that the projection onto the tangent space is a diffeomorphism in a neighbourhood of size the reach of the manifold. This improves upon previous results by Niyogi et al. [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] in terms of the size of the neighbourhood, and is a more quantitative version of results by Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF].

We first recall some notation. Similarly to [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF], we let C (T p M , r 1 , r 2 ) denote the 'filled cylinder' given by all points that project orthogonally onto a ball of radius r 1 in T p M and whose distance to this ball is at most r 2 . We write C (T p M , r 1 , r 2 ) for the open cylinder. We refer to Fig. 3 for an illustration. We now have:

Lemma 3.2 Suppose that M is C 2 and p ∈ M , then for all r < rch M , the projec- tion π T p M onto the tangent space T p M restricted to M ∩ C (T p M , r , rch M ) is a diffeomorphism onto the open ball B T p M (r ) of radius r in T p M , centred at p.
Proof Apart from Lemma 3.1, we will be using the following results from [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF]: For a minimising geodesic γ on M with length parametrised by arc length, with γ (0) = p and γ ( ) = q, we have

∠( γ (0), γ (t)) ≤ t rch M . ( 1
)
If we also write v p = γ (0), we see that

γ ( ), v p = 0 d dt γ (t), v p dt = 0 γ (t), t 0 dt ≥ 0 cos t rch M dt (using (1)) = rch M • sin rch M ≥ rch M • sin ∠(T p M , T q M ) (using Lemma 3.1)
as long as < (π/2) rch M . Because v p ∈ T p M and γ ( ) = q, we have

| p -π T p M (q)| ≥ rch M • sin ∠(T p M , T q M ).
This means in particular that for all q such that | p -π T p M (q)| < rch M and |qπ T p M (q)| ≤ rch M the angle between T p M and T q M is less than 90 degrees. This in turn implies that the Jacobian of projection map in non-degenerate. Note that the condition on mentioned above is satisfied by a combination of Theorem 1 and Lemma 11 of [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF].

It is clear by considering the sphere that this result is tight, in the sense that r cannot be chosen equal to rch M for general manifolds. See Appendix C for some remarks on these results in the C 1,1 setting. Definition 3. [START_REF]Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] We shall write π p as an abbreviation for the restriction (of the domain) of π T p M to M ∩ C (T p M , rch M , rch M ) and π -1 p for its inverse.

We now also immediately have a quantitative version of [52, Lem. IV.8a]:

Corollary 3.4 Suppose that M is C 2 and p ∈ M , then for all r < rch M ,

d( p, M \ C (T p M , r , rch M )) = d p, M \ π -1 p (B T p M (r )) ≥ r . Proof Lemma 3.2 implies that π -1 p (B T p M (r ))) = M ∩ C (T p M , r , rch M )
. By definition of the filled cylinder we have that d( p, R d \ C (T p M , r , rch M )) = r . The result now follows.

We shall also need the following bound on the (local) distance between a tangent space and the manifold. Lemma 3.5 (distance to manifold [START_REF] Boissonnat | The reach, metric distortion, geodesic convexity and the variation of tangent spaces[END_REF]Lem. 11]) Let M be a manifold of positive reach. Suppose that w ∈ T p M and |w -p| < rch M . Let π -1 p (w) be as in Definition 3.3. Then

|π -1 p (w) -w| ≤ ⎛ ⎝ 1 -1 - |w -p| rch M 2 ⎞ ⎠ rch M .
This is attained for the sphere of radius rch M .

Coxeter Triangulations, Delaunay Protection and Stability

Coxeter triangulations [START_REF] Coxeter | Discrete groups generated by reflections[END_REF] of Euclidean space play a significant role in our work. They combine many of the advantages of cubes with the advantages of triangulations. They are also attractive from the geometrical perspective, because they provide simplices with very good quality and some particular Coxeter triangulations are Delaunay protected and thus very stable Delaunay triangulations. We will now very briefly introduce both the concepts of Coxeter triangulations and Delaunay protection, but refer to [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] for more details on Coxeter triangulations and to [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF][START_REF] Boissonnat | Delaunay stability via perturbations[END_REF] for Delaunay protection. This definition imposes very strong constraints on the geometry of the simplices, implying that there are only a small number of such triangulations in each dimension. Most of these triangulations are part of four families for which there is one member for (almost) every dimension d. We will focus on one such family, Ãd , which is Delaunay protected. We refer to Fig. 4 for an illustration of the Ã2 and Ã3 triangulations.

The simplest and shortest definition of a Coxeter triangulation of type Ãd is to give it as a triangulation of a d-dimensional linear subspace of R d+1 by rotation.

Definition 4.2 Let

P = (x i ) ∈ R d+1 | i x i = 0 and consider the d-simplex with vertices u k in P. u 0 = 0 {d+1} , u k = - d + 1 -k d + 1 {k} , k d + 1 {d+1-k} , k ∈ [d],
where x {k} denotes k consecutive coordinates x. The Coxeter triangulation of type Ãd in P is found by consecutively reflecting the simplex in its faces.

Protection

Definition 4. [START_REF]Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] The protection of a d-simplex σ in a Delaunay triangulation on a point set P is the minimal distance of points in P \ σ to the circumscribed ball of σ :

δ(σ ) = inf p∈P\σ d( p, B(σ )),
where B(σ ) is the circumscribed ball of σ . The protection δ of a Delaunay triangulation T is the infimum over the d-simplices of the triangulation: δ = inf σ ∈T δ(σ ). A Delaunay triangulation with a positive protection is called protected.

The proof that Ãd triangulations are protected can be found in [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]Sect. 6]. We shall denote the triangulation of this type by T .

Ã2 Ã3

Fig. [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF] The vertex sets of the Coxeter triangulations in dimensions two and three are the triangular lattice and the body-centred cubic lattice, respectively

Stability

In the triangulation proof below we need that a perturbation T of our initial ambient triangulation (T of type Ãd ) is still a triangulation of R d . We shall refer to this as (combinatorial) stability. Because Whitney did not use a protected Delaunay triangulation, he needs a non-trivial topological argument to establish this, see [52, App. Sect. II.16]. The argument for stability of triangulations for à type Coxeter triangulations is much simpler, because it is a Delaunay triangulation and is δ-protected, see [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]. Before we can recall this result we need to introduce some notation and a definition:

-The minimal altitude or height, denoted by min alt, is the minimum over all vertices of the altitude, that is the distance from a vertex to the affine hull of the opposite face. t(τ ) denotes the thickness of a simplex τ , that is the ratio of the minimal altitude to the maximal edge length. We write t(T ) for infimum of the thickness over all simplices in T . -We can think of the vertices of T as an ( , μ)-net. Here μ is the separation (for Coxeter triangulations, the shortest edge length in T ), and the sampling density (which is the circumradius of the simplices in the Coxeter triangulation). We write μ 0 for the normalised separation, that is μ 0 = μ/ . -For any complex K , L(K ) denotes the longest edge length in K . We use the abbreviations

L = L(T ) and L = L( T ). -A perturbation of the vertices {v i } to { ṽi } is called an -perturbation if |v i -ṽi | ≤ , for all i.
From [13, Thm. 4.14] we immediately get:

Corollary 4. [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF] The triangulation T is (combinatorially) stable under a cL-perturbation as long as

cL ≤ t(T )μ 0 18d δ. ( 2
)
We claim the following concerning the behaviour of c.

Lemma 4.5 c ≤ t(T )μ 0 18d • δ L ≤ √ 2 √ d 2 + 2d + 24 - √ d 2 + 2d 9d 3/2 (d + 1) √ d + 2 ∼ √ 32 3d 4 ,
where ∼ denotes equality up to the leading order in the asymptotic development.

Proof Choudhary et al. [25,App. B] provide explicit values of all the quantities mentioned in Corollary 4.4 for a Coxeter triangulation of type Ã, with the exception of μ, which can be easily derived from a more general result. If we fix the scale (which in [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] we did by a convenient choice of coordinates for the vertices), we have

L(σ ) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ √ d + 1 2 if d is odd, 1 2
d (d + 2) d + 1 if d is even, t(σ ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 2 d if d is odd, 2(d + 1) d (d + 2) if d is even, (3)
= d (d + 2) 12(d + 1) , δ ( σ ) = √ d 2 + 2d + 24 - √ d 2 + 2d √ 12(d + 1
) .

The value of μ easily follows from the general expression for edge lengths (see [25, App. B, Ãd , item 5]) and is equal to

μ = √ d/(d + 1). From (3), we get that μ 0 = μ/ = √ 12/(d + 2). The bound in (2) is therefore c ≤ t(σ )μ 0 18d • δ L = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ √ 2 √ d 2 + 2d + 24 - √ d 2 + 2d 9d 3/2 (d + 1) √ d + 2 if d is odd, 2(d + 1) √ d 2 + 2d + 24 - √ d 2 + 2d 9d 2 (d + 2) 3/2 if d is even. ≤ √ 2 √ d 2 + 2d + 24 - √ d 2 + 2d 9d 3/2 (d + 1) √ d + 2 ∼ √ 32 3d 4 ,
where we used that

√ 1 + x ∼ 1 + x/2 if x is close to zero.

Thickness and angles

The quality of simplices and the control over the alignment of the simplices with the manifold is an essential part of the triangulation proof, for which we need two basic results. Similar statements can be found in [52, Sects. IV.14 and IV.15]. Let us remind the following.

Lemma 4.6 (thickness under distortion [START_REF] Dyer | Riemannian simplices and triangulations[END_REF]Lem. 7]) Suppose that σ = {v 0 , . . . , v k } and σ = {ṽ 0 , . . . , ṽk } are two k-simplices in

R d such that ||v i -v j | -| ṽi -ṽ j || ≤ c 0 L(σ ) for all 0 ≤ i < j ≤ k. If c 0 ≤ t(σ ) 2 /4, then t( σ ) ≥ 4 5 √ k 1 - 4c 0 t(σ ) 2 t(σ ).
123

We can now state a variation of Whitney's angle bound result, see [START_REF] Whitney | Geometric Integration Theory[END_REF]Sect. IV.15].

Lemma 4.7 (Whitney's angle bound) Suppose σ is a j-simplex of R d , j < d, whose vertices all lie within a distance d max from a k-dimensional affine space A

0 ⊂ R d with k ≥ j. Then sin ∠(aff σ, A 0 ) ≤ ( j + 1)d max min alt σ .
Proof We first notice that the barycentre c b of a simplex σ j is at least a distance (min alt σ j )/( j + 1) removed from the faces of the simplex. This means that the ball in aff σ j centred at c with radius (min alt σ j )/( j + 1), denoted by B aff σ j (c, (min alt σ j )/( j + 1)), is contained in σ j . We now consider any diameter, that is a line segment connecting a pair of antipodal points of ∂ B aff σ j (c, (min alt σ j )/ ( j + 1)). This diameter is contained in a d max neighbourhood of A 0 and thus sin ∠( , A 0 ) ≤ ( j + 1)d max min alt σ .

The result now follows, because is arbitrarily chosen.

Simplices in a star in a triangulation of type Ãd

The precise number of simplices in the star of a vertex plays an important role in the volume estimates in Sect. 5. We will now give an explicit bound on this number.

In general the (d -k)-faces of a Voronoi cell correspond to the k-faces in the Delaunay dual. The triangulation T is Delaunay and the dual of a vertex is a permutahedron, see [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]. We recall that the permutahedron is defined as follows: Definition 4.8 (permutahedron) A d-permutahedron is a d-dimensional polytope, which is the convex hull P of all points in R d+1 , the coordinates of which are permutations of {1, . . . , d + 1}.

We also remind the following definition, see [START_REF]Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF], and corollary, see [START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF]: Definition 4.9 Let S(d, k) be the Stirling number of the second kind, which is the number of ways to partition a set of d elements into k non-empty subsets, that is By duality, the lemma immediately gives us the number N k of k-faces that contain a given vertex in T , N k = k! S(d + 1, k). We also write

S(d, k) = 1 k! k j=0 (-1) j k j (k -j) d .
N ≤k = 2 + k j=1 j! S(d + 1, j), ( 4
)
which is an upper bound on the total number of faces of dimension less or equal to k that contain a given vertex. We have added 2 because we want to have a safety margin if we have to consider the empty set (as will be apparent in ( 18)), and have a strict inequality. We now claim the following:

Lemma 4.11 We have N ≤k k d d k .
Proof [START_REF] Rennie | On Stirling numbers of the second kind[END_REF]Thm. 3] gives us that for d ≥ 2 and 1 ≤ j ≤ d -1,

j 2 + j + 2 2 j d-j-1 -1 ≤ S(d, j) ≤ j d-j 2 d j .
Furthermore, Stirling's theorem and the binomial theorem give that j! ∼ j j and k j=0 d j d k , respectively. We now see that

N ≤k = 2 + k j=1 j! S(d + 1, j) k j=1 j! d + 1 j j d+1-j k d k j=1 d j k d d k .
It is clear that if k is much smaller than d that then k d dominates.

Perturbing the Ambient Triangulation

This section is dedicated to the perturbation of the Coxeter triangulation such that the manifold is sufficiently far from the simplices of dimension at most dn -1 in T .

-In Sect. 5.1, we prove that it is possible to perturb the points as described in the second step of part 1 of the algorithm. This involves a significant amount of volume estimates, which are completely quantised. We also indicate how fine the ambient triangulation T has to be compared to rch M ; the longest edge length is linear in terms of the reach (the dependence on the dimension and codimension is rather complicated). -In Sect. 5.2, we define the perturbation and prove that this in fact gives a triangulation for which the low dimensional simplices lie sufficiently far from the manifold.

The proofs of the results in Sect. 5.2 rely on Appendix B. We shall indicate the corresponding sections in Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF], when appropriate.

The Complex T

Before we can dive into the algorithmic construction of the perturbed complex T , we need to fix some constants and give some explicit bounds on them. This subsection corresponds to [52, Sect. IV.18].

Balls and exclusion volumes

Let B d (r ) be any ball in R d of radius r . We now define ρ1 > 0 as follows: For any two parallel (d -1)-hyperplanes whose distance apart is less than 2 ρ1 r , the intersection of the slab between the two hyperplanes with the ball B d (r ) is denoted by S . Now, ρ1 is the largest number such that the volume (vol) of any S satisfies

vol S ≤ vol B d (r ) 2N ≤d-n-1 ,
with N ≤d-n as in [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF]. A precise bound on ρ1 can be given, see Remark 5.2 below. We will use an easier bound ρ 1 , at the cost of weakening the result:

Lemma 5.1 We have ρ1 ≥ ρ 1 = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 2 2k-2 (k!) 2 π(2k)! N ≤d-n-1 if d = 2k, (2k)! 2 2k+2 k!(k -1)! N ≤d-n-1 if d = 2k -1.
(

) 5
Note that

ρ 1 ∼ 1 √ d N ≤d-n-1 .
Proof We can bound the volume of the slab S by the volume of cylinder with base B d-1 (r ) and height 2ρ 1 r , that is

2ρ 1 r d π (d-1)/2 Γ ((d + 1)/2) . This means that vol S vol B d (r ) < 2ρ 1 r • vol B d-1 (r ) vol B d (r ) = 2ρ 1 π (d-1)/2 /Γ ((d -1)/2 + 1) π d/2 /Γ (d/2 + 1) = 2ρ 1 Γ (d/2 + 1) √ π Γ ((d -1)/2 + 1) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ π(2k)! 2 2k-1 (k!) 2 ρ 1 if d = 2k, 2 2k+1 k!(k -1)! (2k)! ρ 1 if d = 2k -1.
using the standard formulae for the volume of the ball, see for example [34, p. 622].

Note that the inequality is strict because ρ 1 > 0. We see that therefore ρ 1 may be chosen to be as in [START_REF] Allgower | An algorithm for piecewise-linear approximation of an implicitly defined manifold[END_REF]. From Wendel's bound on the ratio of Gamma functions [START_REF] Wendel | Note on the gamma function[END_REF], we immediately see that for a fixed constant a, Γ (x + a)/Γ (x) ∼ x a . This means that

2ρ 1 Γ (d/2 + 1) √ π Γ ((d -1)/2 + 1) ∼ 2ρ √ π d 2 + 1 2 1/2 ∼ √ d.
We now see that

ρ 1 ∼ 1 √ d N ≤d-n-1 .
Remark 5.2 Because of symmetry, the largest volume S can attain is when both delimiting hyperplanes are equidistant to the centre of B d (r ). The volume of S is given by the integral

r d ρ 1 -ρ1 vol B d-1 1 -h 2 dh = π (d-1)/2 Γ ((d + 1)/2) r d ρ 1 -ρ 1 1 -h 2 d-1 dh,
where B d-1 (r ) denotes the ball of dimension d -1 with radius r and Γ denotes the Euler gamma function. This integral can be expressed using special functions such as the hypergeometric function or beta functions. This gives an explicit value for ρ1 .

The coarseness of T As mentioned, we perturb the vertices of a Coxeter triangulation. The maximal distance that we allow between an unperturbed vertex v i and the associated perturbed vertex ṽi is cL. We define c as

c = min t(T )μ 0 δ 18d L , t(T ) 2 24 . ( 6
)
The reasons for this particular choice will be discussed after [START_REF] Boissonnat | Sampling and meshing submanifolds in high dimension[END_REF] below. We stress that ( 6) is independent of L because δ scales linearly with L. Notice that because t(T ) ≤ 1, by definition of the thickness of a simplex, we have

c ≤ 1 24 . ( 7
)
We are now ready to introduce the demands on the triangulation of ambient space. We start by bounding the scale of the Coxeter triangulation T by bounding the longest edge length. We do this by giving some constants. We define α 1 and α k by a recursion relation as follows:

α 1 = 4 3 ρ 1 c, 2 3 α k-1 cρ 1 = α k , ( 8
)
that is α k = 2 k+1 ρ k 1 ck /3 k . These definitions play an essential role in the volume estimates for the perturbation of the vertices, that are necessary to guarantee quality. Note that α k is extremely small. In particular, we shall have that

α k ≤ 1 18 k , ( 9
)
because c ≤ 1/24, as we have seen in [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF]. ρ 1 is also very small, as a direct consequence of Lemma 5.1. Furthermore we notice that α k < α k-1 . To make sure the formulae do 123 not become too big, we introduce the notation

ζ = 8 15 √ d d d -n • (1 + 2 c) 1 - 8 c t(T ) 2 t(T ). ( 10
)
Note that ζ depends on both the ambient and intrinsic dimension, and the perturbation parameter c. Because c ≤ t(T ) 2 /24 and t(T ) ≤ 1, we see that ζ ≤ 1. We set the coarseness of the ambient triangulation by demanding that L satisfies

⎛ ⎝ 1 -1 - 6 L(T ) rch M 2 ⎞ ⎠ rch M = α 4+2n d-n 6(n + 1) 2 ζ 2n L, ( 11
)
or equivalently

L rch M = α 4+2n d-n ζ 2n /(3(n + 1) 2 ) (α 4+2n d-n ζ 2n /(6(n + 1) 2 )) 2 + 6 2 . ( 12
)
Note that

L rch M < α 4+2n d-n 54(n + 1) 2 ζ n < α 2 d-n 54 , α 4+2n d-n 6(n + 1) 2 ζ 2n < α 2 d-n 3 ≤ α d-n 3 , ( 13
)
where we used that ζ ≤ 1, which will often be used below to simplify expressions.

Remark 5. [START_REF]Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF] We have to choose the right hand side in [START_REF] Boissonnat | Meshing of surfaces[END_REF] very small, because the bounds on the quality of the simplices that will make up the triangulations are very weak. The details of these estimates can be found in Lemma 6.7.

(dn -1)-skeleton safe triangulations We shall denote the simplices by τ and σ . We will use lower indices to distinguish simplices, while upper indices will stress the dimension, for example τ k j is a simplex of dimension k.

Definition 5.4 ((dn -1)-skeleton safe triangulations) We say that a perturbed triangulation

T of T in R d is (d -n -1)-skeleton safe with respect to the n- dimensional manifold M if d(τ k , M ) > α k L, ( 14
)
for all faces τ k in T , with k ≤ dn -1, and

L < 13 12 L, ( 15
)
t( T ) ≥ 4 5 √ d 1 - 8 c t(T ) 2 t(T ). ( 16
)

Perturbing the Vertices

We now discuss the details of the perturbation scheme that we described in the algorithm section. The perturbation scheme follows Whitney [52, Sect. IV.18] and is inductive.

Construction of T Let v 1 , v 2 , .
. . be the vertices of T . We are going to inductively choose new vertices ṽ1 , ṽ2 , . . . for T , with

|v i -ṽi | ≤ cL = min t(T )μ 0 δ 18d , t(T ) 2 L 24 , ( 17
)
using the notation of Sect. 4. With this bound we have that ( 15) is satisfied, because the two vertices of an edge are perturbed by at most cL and thus the triangle inequality yields L ≤ (1 + 2 c)L. We also claim the following:

Lemma 5.5 T has the same combinatorial structure as T . Moreover, ( 16) is satisfied.

Proof Because we assume that the perturbation is sufficiently small compared to the protection, as given in the first condition of ( 17), ( 2) is satisfied and T will have exactly the same combinatorial structure as T .

By the third condition of [START_REF] Boissonnat | Sampling and meshing submanifolds in high dimension[END_REF] we have a lower bound on the quality of the simplices. To be precise, we have that for any simplex τ in T ,

t(τ ) ≥ 4 5 √ d 1 - 8 c t(T ) 2 t(T ), ( 16
)
as a consequence of Lemma 4.6, the fact that if you perturb the vertices by cL the edge lengths are perturbed by 2 c (that is 2 c = c 0 ), and the fact that if σ ⊂ τ , then t(σ ) ≥ t(τ ). So we have established [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[END_REF].

We now give the scheme where the vertices are perturbed inductively. Suppose that the vertices ṽ1 , . . . , ṽi-1 have been determined, and thus the complex Ti-1 with these vertices. A simplex { ṽ j 1 , . . . , ṽ j k } lies in Ti-1 if and only if {v j 1 , . . . , v j k } lies in T . We shall now find ṽi and thus Ti so that for any τ k ∈ Ti of dimension k ≤ dn -1, ( 14) is satisfied. We distinguish two cases:

Case 1: d(v i , M ) ≥ 3L/2.
In this case we choose ṽi = v i . The inequality ( 14) is established as follows: Because L < (1 + 2 c)L, which means that for any point x in the star of ṽi = v i we have d(x, ṽi

(= v i )) < (1 + 2 c)L. By the triangle inequality we see that d(x, M ) ≥ d(v i , M ) -d(x, ṽi (= v i )) ≥ (1/2 -2 c)L.
That is, any simplex in T with vertex ṽi = v i is at least distance (1/2 -2 c)L from the manifold. Thanks to [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF] we have that (1/2 -2 c)L > 5L/12. This means that d(τ k , M ) > 5L/12 for any simplex in the star of ṽi = v i . This lower bound is much larger than

α k L < L/18 k . Case 2: d(v i , M ) < 3L/2. Let p be a point in M such that d(v i , p) < 3L/2. Let τ 0 (= ∅), τ 1 , . . . , τ ν (18)
123 be the simplices of Ti-1 such that the joins τ j = τ j * ṽi are simplices of T , and dim(τ j * ṽi ) ≤ dn -1 (and thus dim τ j ≤ dn -2), with 0 ≤ j ≤ ν. We note that ν ≤ N ≤d-n-1 , with N ≤k as defined in [START_REF] Allgower | Numerical Continuation Methods: An Introduction[END_REF]. We now consider the span, denoted by span(τ j , T p M ), for all 0 ≤ j ≤ ν. Note that the dimension of span(τ j ,

T p M ) is at most (d -n -2) + n + 1 = d -1.
We now claim the following:

Lemma 5. [START_REF] Behr | Simplex space-time meshes in finite element simulations[END_REF] We can pick ṽi such that it lies sufficiently far from each span(τ j , T p M ), that is

d( ṽi , span(τ j , T p M )) ≥ ρ 1 cL, ( 19
)
while it is not too far from v i , that is | ṽi -v i | ≤ cL.
Proof The argument is volumetric. Let us first introduce the notation U (X , r ) for the set of all points x ∈ R d such that d(x, X ) ≤ r , where X is any subset of R d . By definition of ρ 1 , see 'Balls and exclusion volumes' in Sect. 5.1, and because the dimension of span(τ j ,

T p M ) is at most d -1, we have that vol B(v i , cL) ∩ U (span(τ j , T p M ), ρ 1 cL) ≤ vol B d (r ) 2N ≤d-n-1 . It now follows that vol ⎛ ⎝ B(v i , cL) 1≤ j≤ν U (span(τ j , T p M ), ρ 1 cL) ⎞ ⎠ ≥ vol B(v i , cL) - 0≤ j≤ν vol B(v i , cL) ∩ U (span(τ j , T p M ), ρ 1 cL) > vol B(v i , cL) - 0≤ j≤ν vol B(v i , cL) 2N ≤d-n-1 = 1 - ν + 1 2N ≤d-n-1 vol B(v i , cL) ≥ vol B(v i , cL) 2 ,
where we used that ν ≤ N ≤d-n-1 in the last line, by definition, as mentioned in the description of Case 2. Because the volume is positive we know there exists a point ṽi that satisfies

d( ṽi , span(τ j , T p M )) > ρ 1 cL, (20)
for all 1 ≤ j ≤ ν.

The following lemma completes Case 2:

Lemma 5. [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF] The triangulation T is (dn -1)-skeleton safe, in particular ( 14) is satisfied.

Proof We first make use of the induction3 hypothesis d(τ j , M ) > α k-1 L to find a bound on the distance from τ j to the tangent space T p M , then bound the distance from ṽi * τ j = τ j to T p M based on this. For this argument to work, we have to assume that τ j is not the empty set, that is j = 0. This case is handled separately at the end.

If we combine 1. the induction hypothesis d(τ j , M ) > α k-1 L, 2. the fact that the ball in the tangent space B T p M ( p, r ), centred at p of radius 6L = r , satisfies

B T p M ( p, r ) ⊂ U M , 1 -1 - r rch M 2 rch M ,
thanks to Lemma 3.5, we find that

d(τ j , B T p M ( p, r )) > α k-1 L -1 -1 - r rch M 2 rch M .
This can be simplified:

d(τ j , B T p M ( p, r )) > α k-1 L -1 -1 - r rch M 2 rch M > α k-1 L - α 4+2n d-n 6(n + 1) 2 ζ 2n L (using (11)) > α k-1 L - α d-n L 3 (using (13)) (21)
≥ 2 3 α k-1 L (because α k-1 > α k ).
Because d(v i , p) < 3L/2, L < L + 2 cL, and c < 1/24, see [START_REF] Bendich | Inferring local homology from sampled stratified spaces[END_REF], we have the very coarse bound that

d(τ j , p) ≤ 4L, (22)
by the triangle inequality. We thus find that

d(τ j , T p M \ B T p M ( p, r )) > 2L.
This means that (21) holds for the entire tangent space, that is,

d(τ j , T p M ) > 2 3 α k-1 L. ( 23
)
Lemma B.2, with A 1 = T p M and A 2 = span(τ j , T p M ), now gives

d(τ j , T p M ) ≥ d(τ j , T p M )d(v i , span(τ j , T p M )) L + 2 cL .
This can again be simplified:

d(τ j , T p M ) ≥ d(τ j , T p M )d(v i , span(τ j , T p M )) L + 2 cL > (2α k-1 L/3)ρ 1 cL L + 2 cL
(thanks to ( 23) and ( 19))

> 2L/3 4L/3 α k-1 ρ 1 cL (because c ≤ 1/24) (24)
= 4 3 α k L.
(using the relation ( 8) for α k )

Similarly to [START_REF] Cheng | Dynamic skin triangulation[END_REF], we have that

d(τ j , p) ≤ 4L < 6L.
We can go from the distance from τ j to the tangent space, as given in [START_REF] Cheng | Delaunay Mesh Generation[END_REF], to the distance to the manifold as follows. Because of Corollary 3.4 we can localise the results and Lemma 3.5 allows us to estimate the difference in distance to the manifold and the tangent space. This gives

d(τ j , M ) > 4 3 α k L - ⎛ ⎝ 1 -1 - 6L rch M 2 ⎞ ⎠ rch M .
This can be again simplified: 11) and ( 13))

d(τ j , M ) > 4 3 α k L - ⎛ ⎝ 1 -1 - 6L rch M 2 ⎞ ⎠ rch M > 4 3 α k L - α d-n L 3 (using (
≥ α k L. (because α k ≥ α d-n if k ≤ d -n-1 by (8))
This completes the proof for the case where j = 0 or τ j is non-empty. For j = 0, (20) and Lemma 3.5 yield

d(τ j , M ) > ρ 1 cL - ⎛ ⎝ 1 -1 - 6L rch M
We simplify: 11) and ( 13))

d(τ j , M ) > ρ 1 cL - ⎛ ⎝ 1 -1 - 6L rch M 2 ⎞ ⎠ rch M > ρ 1 cL - α d-n L 3 (using (
> α 1 L.

(by definition of ( 8))

We emphasise that in the perturbation of the points it suffices to look at the tangent spaces at specific points, making this constructive proof an algorithm.

6 Constructing the Triangulation of M Section 6.1 gives geometric consequences of the perturbation we discussed in the previous section. Most importantly we shall see that a simplex σ in T intersects M if and only if it intersects the tangent space T p M of M at a nearby point p close to σ , see Lemma 6.2. Here we again rely on Appendix B. The triangulation K of M is defined in Sect. 6.2.

The Geometry of the Intersection of Simplices in T and M (the Ambient Triangulation and the Manifold)

In this section, we discuss the geometry of simplices in T in relation to M . We follow [52, Sect. IV.19], with the usual exceptions of the use of Coxeter triangulations, the thickness, and the reach to quantify the results. The proofs also differ in a fair number of places from the original. For any p ∈ M we first establish a lower bound on the distance between T p M and simplices in the (dn -1)-skeleton of T that are close to p.

Lemma 6.1 Let p ∈ M and suppose that τ

k ∈ T , with k ≤ d -n -1, be such that τ k ⊂ B( p, 6L). Then d(τ k , T p M ) > 2 3 α k L.
The following proof differs from Whitney's proof.

Proof of Lemma 6.1 Because τ k ⊂ B( p, 6L), the point in T p M that is closest to τ lies in T p M ∩ B( p, 6L) = B T p M ( p, 6L). We now see that

123 d(τ k , T p M ) ≥ d(τ k , B T p M ( p, 6L)) (first sentence of the proof) > d(τ k , M ) - ⎛ ⎝ 1 -1 - 6L rch M 2 ⎞ ⎠ rch M (Lemma 3.5) > α k L - α d-n L 3 (d(τ k , M ) > α k L and (13)) > 2 3 α k L, ( α k ≥ α d-n for k ≤ d -n)
which completes the proof.

We can now examine the relation between intersections with the manifold and nearby tangent spaces.

Lemma 6.2 Suppose that M intersects τ k ∈ T . Let p ∈ M be such that τ k ⊂ B( p, 6(L), then T p M intersects τ k .
Proof Let p ∈ M ∩ τ k . Lemma 3.2 (and ( 11), ( 13)) gives p ∈ π -1 p (B T p M ( p, 6L)), where we use the notation of Definition 3.3. Lemma 3.5 implies

d( p, T p M ) ≤ ⎛ ⎝ 1 -1 - 6L rch M 2 ⎞ ⎠ rch M < α d-n L 3 .
Let τ ⊂ τ k be the face of smallest dimension such that d( τ , T p M ) ≤ 2α d-n L/3. This face exists thanks to the triangle inequality. By Lemma 6.1 we have dim τ ≥ dn. Lemma B.1 implies that τ intersects T p M . The reason for this is the following; τ is the simplex of the smallest dimension such that d( τ ,

T p M ) ≤ 2α k L/3, meaning in particular that d( τ , T p M ) < d(∂ τ , T p M ). Because τ is a face of τ k , clearly T p M intersects τ k .
We can now bound the angle between simplices and tangent spaces. In this case the proof identical to original, and included for completeness.

Lemma 6.3 Suppose that M intersects τ k ∈ T and τ

k has dimension d -n, that is k = d -n. Let p ∈ M be such that τ k ⊂ B( p, 6L), then sin ∠(aff τ k , T p M ) ≥ 2d(T p M , ∂τ k ) L + 2 cL ≥ 4α k L/3 L + 2 cL ≥ 16 13 α k .
Proof This is an immediate consequence of Lemma B.1, [START_REF] Boissonnat | Sampling and meshing submanifolds in high dimension[END_REF], and the previous lemmas.

Below we investigate the relation between intersections of tangent spaces and simplices, and intersections between the manifold and simplices. We combine two statements of [START_REF] Whitney | Geometric Integration Theory[END_REF]Sect. IV.19] in the following lemma. The proof differs from the original by Whitney.

Lemma 6.4 If p ∈ M , τ k ∈ T , τ k ⊂ B( p, 6L), and moreover T p M intersects τ k , then k ≥ dn and M intersects τ k . If k = dn this point is unique, which in particular means that every simplex of dimension dn contains at most one point of M .

Proof Let τ be a face of smallest dimension of τ k such that d( τ , T p M ) ≤ 2α n L/3. Now Lemmas B.1 and 6.2 give that τ and T p M have a unique point p in common and the dimension of τ is dn.

Thanks to Lemma 3.2, M can be written as the graph of a function f , in a neighbourhood of at most size rch M . We note that f :

T p M R n → N p M R d-n
, where here we think of the tangent and normal spaces as embedded in R d . Using the identification of T p M with R n , we now define

F : R × R n → R d , (λ,x) → (x, λ f (x)).
Note that F(0, • ) gives a parametrisation of T p M . Similarly, we can define G : R d-n → R d to be a linear (orthonormal) parametrisation of aff τ . We now consider the difference of the two functions

F -G : R × R n × R d-n = R × R d → R d .
Thanks to Lemma 6.3 we have that sin ∠(aff τ , T p M ) ≥ 16 13 α d-n .

Lemma 3.1 and (13) give that for any q ∈ B( p, 6L)

sin ∠(T p M , T q M ) 2 ≤ 6L 2 rch M ≤ 6 2 • α 2 d-n 54 = α 2 d-n 18 .
It is clear that this also gives an upper bound on the angle between T p M and the graph of F(λ, •) (denoted by graph F(λ, • )) for all λ ∈ [0, 1], due to linearity of the inner product. Because the upper bound on the angle between the tangent spaces is much smaller than the lower bound on ∠(aff τ k , T p M ), aff τ and the tangent space to the graph T q graph F(λ, • ) span R d , for any λ ∈ [0, 1] and q ∈ B( p, 6L). The implicit function theorem and the fact that τ and T p M have a unique point p in common now give that the intersection pλ between graph F(λ, • ) ∩ B( p, 6L) and aff τ exists and is unique, for all λ ∈ [0, 1].

We can now use Lemmas 3.2, 3.5, and 6.3, to bound | p -pλ |. The distance from the manifold to the tangent space is bounded from above by

α 3+2n d-n 3(n + 1) ζ 2n L < α 2 d-n L 3 ,
due to [START_REF] Boissonnat | Meshing of surfaces[END_REF] and [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF]. The same bound holds for graph F(λ, • ). We also have that sin ∠(aff τ , T p M ) ≥ 16α d-n /13. Combining these observations gives

| p -pλ | ≤ α 2 d-n L/3 16α d-n /13 ≤ α d-n L 3 .
This distance bound is smaller than the distance bound for p and the boundary of τ , due to Lemma 6.1. This means that pλ ∈ τ , and in particular that M intersects τ k .

Finally, we study the faces of a simplex that intersects M . This is essential for the barycentric subdivision in part 2 of the algorithm. The proof is identical to the original, but added for completeness. Lemma 6.5 If M intersects τ = {v 0 , . . . , v r } ∈ T , then for each v i ∈ τ , there exists some (dn)-face τ of τ such that v i ∈ τ and τ intersects M . Proof Take p ∈ M ∩ τ . Let τ k be a face of the smallest dimension of τ , with v i ∈ τ k , that intersects T p M . Now assume that k > dn. Let us write τ k-1 for the face of τ k opposite v i . Because the dimension of τ k ∩ T p M is at least 1, the intersection of T p M and τ k-1 is non-empty.

Similarly to the first argument in the proof of Lemma 6.4, we see that T p M intersects some (dn)-face of τ k-1 . Thanks to Lemma 6.3, the angle between this (d -n)-face and T p M is bounded from below. Due to Lemma 6.1, the intersection lies in the interior of the (dn)-face. The angle bound and the fact that the intersection lies in the interior gives that any simplex in T that contains this (dn)-face has points in the interior that lie in T p M . In particular, the interior of τ k contains part of T p M . Because both the interior of τ k and τ k-1 contain points of T p M , linearity gives that T p M must intersect ∂ τ k \ τ k-1 . From this contradiction of the assumption, we conclude that k = dn. Lemma 6.4 finally says that M intersects τ k , because T p M does.

The Triangulation of M : The Complex K

The construction of the complex follows [START_REF] Whitney | Geometric Integration Theory[END_REF]Sect. IV.20].

In each simplex τ of T that intersects M , we choose a point v(τ ) and construct a complex K with these points as vertices. The construction goes via barycentric subdivision of general polytopes or even CW-complexes, see for example [START_REF] Lundell | The Topology of CW Complexes[END_REF]Thm. 1.7

of Chapter III]. For each sequence τ

0 ⊂ τ 1 ⊂ . . . ⊂ τ k of distinct simplices in T such that τ 0 intersects M , σ k = {v(τ 0 ), . . . , v(τ k )} ( 25
)
will be a simplex of K . The definition of v(τ ) depends on the dimension of τ :

-If τ is a simplex of dimension dn, then there is an unique point of intersection with M , due to Lemma 6.4. We define v(τ ) to be this unique point. -If τ has dimension greater than dn, then we consider the faces τ d-n 1 , . . . , τ d-n j of τ of dimension dn that intersect M . These faces exist thanks to Lemma 6.5. We now define v(τ ) as follows:

v(τ ) = v(τ d-n 1 ) + • • • + v(τ d-n j ) j . ( 26
)
Remark 6. [START_REF] Behr | Simplex space-time meshes in finite element simulations[END_REF] We stress that thanks to Lemma 6.4, choosing the point v(τ d-n ) to be the point of intersection with T p M , assuming p is sufficiently close, locally gives the same combinatorial structure as intersections with M . We also stress that for the combinatorial structure it does not really matter where M intersects a simplex of T , as long as it does.

We can now state the following bound on the altitudes of the simplices we constructed in this manner. Lemma 6.7 Let σ n be a top dimensional simplex as defined in [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF], then

min alt σ n > ζ n (α d-n-1 ) n L,
where min alt denotes the minimal altitude or height, and we used the notation ζ as defined in [START_REF] Boissonnat | Geometric and Topological Inference[END_REF].

Proof This inequality relies on estimates on the barycentric coordinates and Lemma 4.6. We first establish a bound on the barycentric coordinates of v(τ d-n i ) for some

(d -n)-dimensional simplex τ d-n i ∈ T that intersects M . By Lemma 6.1, v(τ d-n i ) lies at least a distance 2α d-n-1 L/3 from the boundary ∂τ d-n i
, and the longest edge is at most L + 2 cL. This means that all the barycentric coordinates λ l with respect to (the vertices of) τ d-n i are at least ) satisfy the bound [START_REF] Von Danwitz | Simplex space-time meshes in compressible flow simulations[END_REF], while the other n coordinates are zero. This also means that for the barycentric coordinates with respect to

λ l (τ d-n i ) > 2 3 α d-n-1 L L + 2 c L = 2 3 α d-n-1 1 1 + 2 c . ( 27
τ d of v(τ k ) = v(τ d-n 1 ) + • • • + v(τ d-n j ) j ,
for k > dn, we have that:

k + 1 of the coordinates λ l satisfy

λ l > 2 3 j • α d-n-1 1 + 2 c .
-The other dk coordinates are zero.

Note that j ≤ d d-n . This means that 5 Both triangles are right angled. We stress that the projection of p onto aff σ may be quite a distance from σ itself

d(v(τ k ), ∂τ k ) ≥ 2α d-n-1 p d min aff L h σ ′ σ ′ ′ Fig.
We now have

d(v(τ k ), ∂τ k ) ≥ 2α d-n-1 t( T ) L 3 d d -n • (1 + 2 c)
(by definition of thickness) 2 (by the estimate ( 16))

≥ 2α d-n-1 t(T ) L 3 d d -n • (1 + 2 c) • 4 5 √ d 1 - 8 c t(T )
≥ 8α d-n-1 t(T ) L 15 √ d d d -n • (1 + 2 c) 1 - 8 c t(T ) 2 . ( 28
)
Using this estimate and the fact that σ n is defined through a sequence

τ d-n 0 ⊂ τ d-n+1 1 ⊂ . . . ⊂ τ d n ,
we can give a lower bound on the minimal altitude of the simplex. We are going to use the following easy observation on the minimal altitude simplices. Suppose that:

-The simplex σ is the join of a point p and the simplex σ .

-d( p, aff σ ) ≥ d min . -min alt σ ≥ h . -The maximum edge length of σ is L(σ ).
Then min alt σ ≥ h d min /L(σ ), as can be established by simple trigonometric arguments; see Fig. 5. Applying this result n times gives that

min alt σ n ≥ h(σ n-1 )d min (σ n ) L(σ n ) ≥ h(σ n-2 )d min (σ n )d min (σ n-1 ) L(σ n ) 2 ≥ d min (σ n ) . . . d min (σ 0 ) L(σ n ) n , ( 29
)
where we indicated the dimensions explicitly.

Plugging (28) (using the definition of the simplex σ n as in ( 25)) into [START_REF] Dey | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF] gives that min alt σ n is lower bounded as follows:

min alt σ n > ⎛ ⎜ ⎜ ⎝ 8α d-n-1 t(T ) 15 √ d d d -n • (1 + 2 c) 1 - 8 c t(T ) 2 ⎞ ⎟ ⎟ ⎠ n L = ζ n α n d-n-1 L,
which completes the proof.

The Triangulation Proof

Given the triangulation T , we want to prove that the intersection of M ∩ τ d is homeomorphic to the triangulated polytope described in Sect. 6.2. This immediately gives a global homeomorphism between the triangulation and the manifold. The homeomorphism we discuss in this section differs greatly from Whitney's own approach. Firstly, he used the closest point projection as a map (which does not respect simplices, meaning that the point in the complex K (as defined in the previous section) and its projection may lie in different simplices of T . Secondly, to prove that this map is a homeomorphism, he uses what has become known as Whitney's lemma in much the same way as in [START_REF] Boissonnat | Local criteria for triangulation of manifolds[END_REF].

The great advantage of our approach to the homeomorphism proof is that it is extremely explicit and it is elementary in the sense that it does not rely on topological results. We also need precise bounds on the angles, which do not require deep theory, but are quite intricate.

Because we work with an ambient triangulation of type à and we do not perturb too much, the simplices of T are Delaunay. The homeomorphism from M ∩ τ d to the triangulated polytope K ∩ τ d , with K as defined in Sect. 6.2 and τ d ∈ T , gives that the intersection of any simplex in T with M is a topological ball of the appropriate dimension. This may remind the reader of the closed ball property of Edelsbrunner and Shah [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]. We stress that the homeomorphism we construct is explicit.

Overview of the homeomorphism proof The proof consists of three steps:

-For each d-simplex τ ∈ T we provide a 'tubular neighbourhood' for K ∩ τ adapted to τ . By this we mean that, for each point p in K ∩ τ , we designate a 'normal' space N p that has dimension equal to the codimension of M and K , and is transversal to K ∩ τ . Moreover, these directions shall be chosen in a sufficiently controlled and smooth way, so that every point x in τ that is sufficiently close to K has a unique point p on K ∩ τ such that x -p ∈ N p. -We give conditions that enforce that the 'normal' spaces N p intersect M transversely. More precisely, we prove that the angle between Ñ p and N q M , for any q ∈ M ∩ τ , is upper bounded by a quantity strictly less than 90 degrees. -We conclude that the projection along N p gives a homeomorphism from M to K .

Constructing the Tubular Neighbourhood

We now give the construction of a 'tubular neighbourhood' of K . We refer to Fig. 6 for a pictorial overview of the construction. We use two results from the previous sections:

-The normal space is almost constant, see Lemma to N p M is a (linear) bijection, for any p that is sufficiently close to τ d-n k , with k ∈ {1, . . . , j}. We will denote the inverse of this map by π -1

τ d-n k .
We can now define the 'normal spaces' for the complex K . We first do this for the vertices v(τ ), where τ has dimension dn (these vertices lie on M ), secondly for general vertices of K (these vertices do not necessarily lie on M ), and finally, using barycentric coordinates, for arbitrary points in K .

We start, as mentioned, with the vertices that are associated to a simplex τ = τ d-n ∈ T of dimension d -n. We stress that there is one face of

τ of dimension d -n, namely τ itself, so τ = τ d-n = τ d-n 1 . For the vertex v(τ ) = v(τ d-n ) = v(τ d-n 1 ) we choose the normal space N v(τ d-n 1 ) to be aff τ d-n 1 .
For v(τ ) such that the dimension of τ is greater than dn we make the following construction, which is reminiscent of the construction of v(τ ) in Sect. 6.2. Let p ∈ M be such that τ ⊂ B( p, 6L). For now p is arbitrary, we will specify this later. We consider the faces

τ d-n 1 , . . . , τ d-n j of τ of dimension d -n that intersect M . Now consider the orthogonal projection map π v(τ d-n k ), p : aff τ d-n k → N p M . For any w ∈ N p M we define N τ, p (w) = 1 j π -1 v(τ d-n 1 ), p (w) + • • • + π -1 v(τ d-n j ), p (w) . ( 30
)
To construct the normal space at v(τ ), we pick p = π M (v(τ )) and define the normal space as

N v(τ ) = span N τ,π M (v(τ )) (w). Let σ n = {v(τ d-n 0 ), . . . , v(τ d n )
} be a simplex of K . Now choose a point p ∈ M as before. For any point p in σ n with barycentric coordinates λ = (λ 0 , . . . , λ n ), and any w ∈ N p M , we define

N p, p (w) = λ 0 N τ d-n 0 , p (w) + • • • + λ n N τ d n , p (w).
We now set p = π M ( p). By defining N p = span N p,π M ( p) (w), we get affine spaces for each point in each σ n ∈ K .

Remark 7.1 By construction, these spaces are consistent on the faces of simplices in K as well as with the boundaries of the d-dimensional simplices in T .

Fig. 6 The tubular neighbourhood

The tubular neighbourhood is defined as the set of all points in R d that that lie in a unique N p, with p ∈ K .

The Size of the Tubular Neighbourhoods and the Homeomorphism

In this section, we establish the size of the neighbourhood of K as defined by N p. The following angle estimate is an essential part of the estimate of the size of the neighbourhood of the triangulation K . Lemma 7.2 Suppose that p ∈ M , τ d ⊂ B( p, 6L), and σ n ∈ K are such that σ n ⊂ τ d , where we regard σ n and τ d as subsets of R d . Then the angle between T p M and aff σ n is bounded as follows:

sin ∠(aff σ n , T p M ) ≤ α 4+n d-n 6(n + 1) ζ n .
Proof By Lemma 4.7, we have

sin ∠(aff σ n , T p M ) ≤ (n + 1)d max min alt σ n ,
where d max denotes the maximum distance of the vertices of σ n to T p M . Lemma 6.7 gives us the following bound:

min alt σ n > α n d-n-1 ζ n L.
Finally, d max is bounded thanks to [START_REF] Boissonnat | Meshing of surfaces[END_REF]. Combining these results yields sin ∠(aff

σ n , T p M ) ≤ (n + 1)α 4+2n d-n ζ 2n L/(6(n + 1) 2 ) α n d-n-1 ζ n L ≤ α 4+n d-n ζ n 6(n + 1)
,

123
because α d-n-1 < α d-n and L ≥ L (there are unperturbed simplices in T ).

With this we can give a bound on the size of the neighbourhood of K .

Lemma 7.3 Let p, q ∈ σ n , with barycentric coordinates λ = (λ 0 , . . . , λ n ), λ = (λ 0 , . . . , λ n ), respectively. Suppose that N p and N q are defined as in Sect. 7.1. Suppose now that the intersection between p +N p and q +N q is non-empty. Here p +N p and q + N q denote the affine spaces that go through p, q and are parallel to N p, N q , respectively.

If x ∈ p + N p ∩ q + N q , then d(x, aff σ n ) ≥ 15 13 2 α 4 d-n n + 1 ζ n α n d-n-1 L. (31)
Because, by construction, the N p agree on the faces of the n-dimensional simplices in K , this provides a tubular neighbourhood for K of the size indicated in the right hand side of [START_REF] Dey | Localized Delaunay refinement for piecewise-smooth complexes[END_REF].

Proof The main idea of the proof of this lemma is the following: Given two points p, q ∈ σ n ⊂ K , the 'normal' spaces N p and N q are not intersecting too close to K if the angle between N p and N q is not too large compared to the distance between p and q (and the angle between N p and aff σ is not too small). The proof consists of several steps. Step 0 gives some very rough estimates, mainly on the angles between the various 'normal' spaces of K that we construct and N p M . Steps 1, 2, and 3 work from these very naïve bounds to fairly sharp estimates on ∠(N p, N q ). In the fourth and final step the bound on ∠(N p, N q ) is used to give a lower bound on the size of the tubular neighbourhood.

Step 0: preliminary estimates. Lemma 6.3 gives that for each

τ d-n sin ∠(aff τ d-n , T p M ) ≥ 16 13 α d-n or cos ∠(aff τ d-n , N p M ) ≥ 16 13 α d-n , so that for u ∈ N p M of unit length, cos ∠ π -1 v(τ d-n k ), p (u), u ≥ 16 13 α d-n , with π v(τ d-n k ), p : aff τ d-n k → N p M the orthogonal projection map. This means that π -1 v(τ d-n k ), p (u) ≤ 13/(16α d-n ).
Together with the triangle inequality this yields that

|N τ, p (u)|, |N p(u )| ≤ 13 16α d-n , (32)
for any u ∈ N p M and u ∈ N π M ( p) M of unit length. By construction, the component of N τ, p (u) in the u direction is u and the component of N p(u ) in the u direction is also u . This in turn gives us that

∠(N τ, p (u), u), ∠(N p(u ), u ) ≤ arccos 16α d-n 13 . (33)
Thus

∠ span u∈N p M N τ, p (u), N p M , ∠(N p, N π M ( p) M ) ≤ arccos 16α d-n 13 . ( 34
)
If we want to compare the two different normal spaces N p M and N q M , with | p-q| ≤ 4L, we again invoke Lemma 3.1 and ( 13) to see that

sin ∠(N p M , N q M ) 2 ≤ 2L rch M < α 4+2n d-n 29(n + 1) 2 ζ n .
Using ( 9) and the fact that ζ is small, we can further simplify:

sin ∠(N p M , N q M ) 2 < α d-n 18 3 • 29 = α d-n 169128 .
The triangle inequality for angles (or points on the sphere) now implies that

∠ span u∈N p M N τ, p (u), span u∈N q M N τ,q (u) , ∠(N p, N q ) < 2 arccos 16α d-n 13 + arcsin α d-n 18 3 • 29 .
Overview steps 1, 2, and 3: angle estimates. Having established some preliminary estimates, we will tighten this result for ∠(N p, N q ). The angle between these two terms is determined by both p and p in N p, p (u). We will examine the effects of both separately.

Step 1: Bounding ∠(N p, p (u), N q, p (u)). We start by fixing p and varying p. We now consider

N p, p (u) = λ 0 N τ d-n 0 , p (u) + • • • + λ n N τ d n , p (u) and N q, p (u) = λ 0 N τ d-n 0 , p (u) + • • • + λ n N τ d n , p (u).
We are now going to estimate the angle between these vectors and thus the angle between span u N p, p (u) and span u N q, p (u) in terms of the barycentric coordinates.

The u components of N p, p (u) and N q, p (u) are u by construction as we mentioned before. We are going to compare this with the length of N p, p (u) and N q, p (u), and the length of their difference. For estimates on these lengths we need to introduce the following notation:

≤ 13 √ n + 1 16 D d-n |O -O | 1 ≥ 1 T 0 I 0
Fig. 7 The worst case for the angle between the vectors N p (e j ) and N q (e j ). We write φ 0 for an upper bound on ∠(N p (e j ), N q (e j )). Moreover θ 0 ≥ arcsin (16α d-n /10). The length or bound on the length of two of the edges is also indicated in the figure With this notation, we can now derive the following bound:

λ 0 N τ d-n 0 , p (u) + • • • + λ n N τ d n , p (u) -λ 0 N τ d-n 0 , p (u) + • • • + λ n N τ d n , p (u) = (λ 0 -λ 0 ) N τ d-n 0 , p (u) + • • • + (λ n -λ n ) N τ d n , p (u) = N τ d-n 0 , p (u). . .N τ d n , p (u) ⎛ ⎜ ⎝ λ 0 -λ 0 . . . λ n -λ n ⎞ ⎟ ⎠ ≤ N τ d-n 0 , p (u). . .N τ d n , p (u) 2 |λ -λ | (35)
≤ N τ d-n 0 , p (u). . .N τ d n , p (u) F |λ -λ | (because • 2 ≤ • F ) = N τ d-n 0 , p (u) 2 + • • • + N τ d n , p (u) 2 • |λ -λ | (by definition of • F ) ≤ 13 √ n + 1 16α d-n |λ -λ |.
(by [START_REF] Dey | Normal and feature approximations from noisy point clouds[END_REF])

We now turn our attention to the triangle with edges N p, p (u), N q, p (u), and N p, p (u) -N q, p (u), as depicted in Fig. 7. We apply the sine rule to this triangle, using [START_REF] Doi | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[END_REF] and [START_REF] Dyer | Riemannian simplices and triangulations[END_REF], to find sin ∠(N p, p (u), N q, p (u))

≤ sin φ 0 = 13 √ n + 1 16α d-n • |λ -λ | • 13 16α d-n = 13 16α d-n
Note that this can be tightened a fair bit at the cost of complicating the bound. We conclude from [START_REF] Edelsbrunner | Triangulating topological spaces[END_REF]

that sin ∠ span u N p, p (u), span u N q, p (u) ≤ 13 16α d-n 2 √ n + 1 • |λ -λ |. ( 37
)
Step 2: bounding ∠ span u∈N p M N p, p (u), span u∈N q M N p,q (u) . We now want to bound the angle between span u∈N p M N p, p (u) and span u∈N q M N p,q (u) based on the distance between the points p and q in M . We use the fact that p and q are such that τ ⊂ B( p, 6L), B(q, 6L), so the conditions of Lemma 6.3 hold. This also means that p and q are close, so the angle between N q M and N p M is very small. This gives that the projection π N q M →N p M induces a (linear) bijection from N q M to N p M , so that the inverse π -1 N q M →N p M makes sense. Having established this map, we see that

∠ span u∈N p M N p, p (u), span u∈N q M N p,q (u) ≤ sup u∈N p M ∠ N p, p (u), N p,q π -1 N q M →N p M (u) . (38)
To bound this angle, we look at the individual terms in [START_REF] Dey | Localized Delaunay refinement for volumes[END_REF], that is π -1

v(τ d-n k ), p (u)
and π -1

v(τ d-n 1
),q π -1 N q M →N p M (u) . See Fig. 8 for an illustration. We can write π -1 N q M →N p M (u) = u + wq,p , with wq,p ∈ T p M , | wq,p | ≤ tan ∠(N p M , N q M ), and

|u + wq,p | ≤ 1 cos ∠(N p M , N q M ) . Similarly, π -1 v(τ d-n k ), p
(u) can be written as u + wk,p , with wk,p ∈ T p M , and

| wk,p | ≤ tan ∠(aff τ d-n k , N p M ). Likewise, π -1 v(τ d-n k
),q (u + wq,p ) can be written as u + wq,p + wk,q , with wk,q ∈ T q M , and

| wk,q | ≤ tan ∠(aff τ d-n k , N q M ) cos ∠(N p M , N q M ) .
123

aff k π -1 ( k d -n ), q π -1 π -1 ( k d -n ), p (u + q,p ) - (u) u N p N p N q N q (u) = u + q,p - → w w M M M M τ τ τ v v
Fig. 8 Note that the two-dimensional nature of the figure is slightly misleading

The distance from π -1

v(τ d-n k
),q (u + wq,p ) to the translation of T p M that goes through (u + wq,p ) = π -1 N q M →N p M (u) are both contained in aff τ d-n k , we have that

u is at most tan ∠(aff τ d-n k , N q M ) sin ∠(T p M , T q M ) cos ∠(N p M , N q M ) = tan ∠(aff τ d-n k , N q M ) sin ∠(N p M , N q M ) cos ∠(N p M , N q M ) .
π -1 v(τ d-n k ), p (u) -π -1 v(τ d-n k ),q (u + wq,p ) = π -1 v(τ d-n k ), p (u) -π -1 v(τ d-n k ),q π -1 N q M →N p M (u) ≤ tan ∠(aff τ d-n k , N q M ) sin ∠(N p M , N q M ) cos ∠(N p M , N q M ) sin ∠(aff τ d-n k , T p M ) (39) = tan ∠(aff τ d-n k , N q M ) sin ∠(N p M , N q M ) cos ∠(N p M , N q M ) cos ∠(aff τ d-n k , N p M ) = tan ∠(aff τ d-n k , N q M ) tan ∠(N p M , N q M ) cos ∠(aff τ d-n k , N p M ) .
Lemma 6.3 gives us that sin ∠(aff

τ d-n k , T p M ), sin ∠(aff τ d-n k , T q M ) ≥ 16α d-n /13, so cos ∠(aff τ d-n k , T q M ) ≥ 16α d-n /13 and tan ∠(aff τ d-n k , N q M ) ≤ 1 -(16α d-n /13) 2 16α d-n /13 . By Lemma 3.1 sin ∠(T p M , T q M ) 2 ≤ | p -q| 2 rch M , so that, using the identity tan (2 arcsin x) = 2x √ 1 -x 2 /(1 -2x 2 ), tan ∠(N p M , N q M ) ≤ | p -q| rch M 1 - | p -q| 2 4(rch M ) 2 1 - | p -q| 2 2(rch M ) 2 .
This means that (39) yields

π -1 v(τ d-n k ), p (u) -π -1 v(τ d-n k ),q π -1 N q M →N p M (u) ≤ 1 - 16 13 α d-n 2 16 13 α d-n 2 • | p -q| rch M 1 - | p -q| 2 4(rch M ) 2 1 - | p -q| 2 2(rch M ) 2 . ( 40
)
By the triangle inequality (applied to the terms in the sum in the definition [START_REF] Dey | Localized Delaunay refinement for volumes[END_REF]), this gives N p, p (u) -N p,q π -1 N q M →N p M (u) ≤ the right hand side of [START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF].

Because the u component of N p, p (u)) is u and its length is at least 1, we find that also sin ∠ N q, p (u), N q,q π -1 N q M →N p M (u) ≤ the right hand side of [START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF].

Because p and q are very close (in fact, they lie in an L neighbourhood of τ , which is small due to ( 13))

1 - | p -q| 2 4(rch M ) 2 1 - | p -q| 2 2(rch M ) 2 ≤ 2,
as can verified using the fact that √ 1x 2 /(1 -2x 2 ) is monotone increasing for sufficiently small x. This means we can simplify the result further: sin ∠ N q, p (u), N q,q π -1

N q M →N p M (u) ≤ 10 16α d-n 2 | p -q| rch M .
Thanks to [START_REF] Hirsch | Differential Topology[END_REF], we now have

∠ span u∈N p M N q, p (u), span u∈N q M N q,q (u) ≤ arcsin 13 16α d-n 2 | p -q| rch M . ( 41
)
Step 3: Combining into a bound on ∠(N p, N q ). Combining ( 41), [START_REF] Federer | Curvature measures[END_REF], and the triangle inequality, we see that

∠ span u∈N p M N p, p (u), span u∈N q M N q,q (u) ≤ ∠ span u∈N p M N p, p (u), span u∈N p M N q, p (u) + ∠ span u∈N p M N q, p (u), span u∈N q M N q,q (u) ≤ arcsin 13 16α d-n 2 √ n + 1 • |λ -λ | + arcsin 13 16α d-n 2 | p -q| rch M .
Because we need estimates on ∠(N p, N q ) we have to set p = π M ( p) and q = π M ( q). To estimate the distance between the two points, we first note that | p -q| ≤ |λ -λ | L, because p and q have barycentric coordinates λ and λ . Thanks to [START_REF] Federer | Curvature measures[END_REF]Theorem 4.8 (8)], we have that if

d(x, M ), d(y, M ) ≤ (rch M )/2, then |π M (x) -π M (y)| ≤ 2|x -y|. This means that ∠(N p, N q ) = ∠ span u∈N π M ( p) M N p,π M ( p) (u), span u∈N π M ( q) M N q,π M ( q) (u) ≤ arcsin 13 16α d-n 2 √ n + 1 • |λ -λ | + arcsin 13 16α d-n 2 2|λ -λ | L rch M . Since sin (arcsin x + arcsin y) = y √ 1 -x 2 + x 1 -y 2 ≤ x + y, we see that sin ∠(N p, N q ) ≤ 13 16α d-n 2 √ n + 1 • |λ -λ | + 13 16α d-n 2 2|λ -λ | L rch M .
Because of ( 13) and ( 7),

13 16α d-n 2 2|λ -λ | L rch M ≤ 13 16α d-n 2 4|λ -λ | α 4+2n d-n 54(n + 1) 2 ≤ 6 100 |λ -λ |α 2 d-n .
Thanks to [START_REF] Berger | State of the art in surface reconstruction from point clouds[END_REF], the first term in the following sum is by far the larger one:

sin ∠(N p, N q ) ≤ 13 16α d-n 2 √ n + 1 • |λ -λ | + 6 100 |λ -λ |α 2 d-n .
We finally arrive at the following simple, but weaker bound:

sin ∠(N p, N q ) ≤ 13 15α d-n 2 √ n + 1 • |λ -λ |. (42)
Step 4: From angles to a lower bound on the neighbourhood size. We now consider the triangle p qx and we estimate | px| and | qx|. Recall that in the statement of the lemma we defined x as the point where the normal spaces N p and N q first intersect. The estimate will use:

1. the sine rule; 2. the fact that the distance between p and q is at least |λ -λ | min alt σ/ √ n, thanks to [10, Lem. 5.12]; 3. Lemma 6.7 to bound min alt σ ; 4. inequality [START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF], which gives a bound on the angle ∠ px q, namely φ 0 . 5. Lemma 7.2 gives that

∠(N p M , (aff σ n ) ⊥ ) ≤ arcsin α 4+n d-n ζ n 6(n + 1) ≤ arcsin α 4 d-n 6 ,
where (aff σ n ) ⊥ denotes the space perpendicular to aff σ n . Because e j ∈ N p M , combining this with Lemma 6.3 and the triangle inequality for angles yields

∠(N p (e j ), (aff σ n ) ⊥ ) ≤ ∠(N p (e j ), e j ) + ∠(N p M , (aff σ n ) ⊥ ) ≤ arccos 16α d-n 13 + arcsin α 4 d-n 6 . ( 43
)
We need a lower bound on sin ∠ p qx and sin ∠ q px, that is sin ∠(N p (e j ), aff σ n ) = cos ∠(N p (e j ), (aff σ n ) ⊥ ).

We also recall the trigonometric identity cos (arccos a

+ arcsin b) = a 1 -b 2 -b 1 -a 2 .
123

Using [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] (as α k ≤ 1/18 k by ( 9))

≥ α d-n .
This completes the fifth point.

The considerations we summed up yield

| px|, | qx| ≥ α d-n |λ -λ | • min alt σ/ √ n (10/(15α d-n )) 2 • √ n + 1 • |λ -λ | ≥ α d-n (15α d-n /10) 2 • min alt σ n + 1 ≥ α d-n (15α d-n /13) 2 n + 1 (ζ α d-n-1 ) n L.
Using [START_REF] Oudot | Meshing volumes bounded by smooth surfaces[END_REF] again yields that the distance from x to aff σ n is bounded from below by

d(x, aff σ n ) ≥ α 2 d-n (15α d-n /13) 2 n + 1 (ζ α d-n-1 ) n L.
This completes the proof.

n-dimensional simplex σ n ⊂ K ∩ τ d . Note that v(τ d ) ∈ K ∩ τ d . Thanks to Lemma 7.2, sin ∠(aff σ n , T v M ) ≤ α 4+n d-n 6(n + 1) ζ n .
From this bound we conclude that

d H T v M ∩ B(v, 2L), aff σ n ∩ B(v, 2L) ≤ 2 α 4+n d-n 6(n + 1) ζ n L,
where d H denotes the Hausdorff distance. Because of Lemma 3.5 and (11), we have that

d H T v M ∩ B(v, 2L), π -1 v (B T v M (v, 2L)) ≤ α 4+2n d-n 6(n + 1) 2 ζ 2n L,
where B T v M (v, 2L) denotes the ball in T v M with centre v and radius 2L. This gives us

d H aff σ n ∩ B(v, 2L), π -1 v (B T v M (v, 2L)) ≤ 2 α 4+n d-n 6(n + 1) ζ n L + α 4+2n d-n 6(n + 1) 2 ζ 2n L (by the triangle inequality) ≤ α 4+n d-n ζ n n + 1 L. Because M ∩ τ ⊂ π -1 v (B T v M (v, 2L
)) and the distance between M ∩ τ and aff σ n is small compared to the size of the neighbourhood of K given in Lemma 7.3, that is

α 4+n d-n ζ n n + 1 L ≤ (15/13) 2 α 4 d-n n + 1 ζ n α n d-n-1 L, M ∩ τ is contained in this neighbourhood of K .
Having established that M lies in the tubular neighbourhood around K , it is meaningful to speak about the projection from M to K along the direction N . Because we also have that the projection from M to K in the direction N (as defined in Sect. 7.1) is transversal (because π/2 minus the angle between N p and N p M , see [START_REF] Duistermaat | Multidimensional Real Analysis. II. Integration[END_REF], is much bigger than the variation of the tangent/normal space as bounded by Lemma 3.1 and ( 12)), we see that M ∩ τ d is homeomorphic to K ∩ τ d . By construction the projection map is compatible on the boundaries of τ d , so we also immediately have an explicit homeomorphism between M and K . Moreover, this homeomorphism is piecewise smooth and not just continuous. This completes the proof of Theorem 1.1. We emphasise that along the way we have also given bounds on -the Hausdorff distance between M and K , see Lemmas 7.4 and 7.3, -the quality of simplices, see Lemma 6.7, -the variation of the tangent spaces, see Lemma 3.1, (34), and (12). material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix A: Notation

In the following table we give an overview of the notation used in this paper and compare it to Whitney's notation.

Overview of the Most Important Bounds

We recall here for the reader's convenience the most important bounds and constants used in the paper. The constant ρ1 > 0 (depending only on d) is defined as follows: For any two parallel (d -1)-hyperplanes whose distance apart is less than 2 ρ1 r , the intersection of the slab between the two hyperplanes with the ball B d (r ) is denoted by S . Now, ρ1 is the largest number such that the volume (vol) of any S satisfies vol S ≤ vol B d (r ) 2N ≤d-n-1 ,

where N ≤k is an upper bound on the total number of faces of dimension less or equal to k that contain a given vertex, see (4). α 1 and α k have been defined by a recursion relation as follows: The notation is illustrated in Fig. 9.

α 1 = 4 3 ρ 1 c, 2 3 α k-1 cρ 1 = α k , ( 8
Proof Choose p ∈ σ and q ∈ A 0 such that

| p -q| = d(A 0 , σ ).
Now suppose that there is a vector v = 0 that lies in the intersection of aff σ and A 0 . Then there exists some c ∈ R such that p + cv ∈ ∂σ . Because v lies in the intersection of aff σ and A 0 , we have that q + cv ∈ A 0 . As translation leaves distances invariant,

d(A 0 , σ ) = |p -q| = |( p + cv) -(q + cv)| ≥ d(A 0 , ∂σ ),
which clearly contradicts the assumption. This means we can conclude that there is no such v and therefore s + n = d.

Because there is no v in the intersection of aff σ and A 0 , there is a unique point p in this intersection. We will now show that p ∈ σ . We will assume that p / ∈ σ . This means in particular that q = p. Because d(A 0 , σ ) < d(A 0 , ∂σ ), pq is normal to aff σ and p ∈ σ \ ∂σ . Now consider the line from q to p, which lies in A 0 . The distance from a point on this line to σ decreases (at least at first) as you go from q toward p. This contradicts the definition of q. We conclude that p ∈ σ . Now suppose that l 0 is a line in A 0 that goes though p. In order to derive a contradiction, we assume that sin φ < 2d(A 0 , ∂σ ) L(σ ) ,

where sin φ denotes the angle between l 0 and aff σ . Denote by π aff σ (l 0 ) the orthogonal/closest point projection on aff σ of l 0 . Because p ∈ σ , π aff(σ ) (l 0 ) intersects ∂σ at a point q and we may assume that | p -q| ≤ L(σ )/2 so that l 0 contains a point of distance

L(σ ) sin φ 2 < L(σ ) 2 • 2d(A 0 , ∂σ ) L(σ ) = d(A 0 , ∂σ )
from ∂σ , a contradiction. Because l 0 was an arbitrary line in A 0 the result now follows.

The following is a variation on Lemma 14b from [52, App. II.14]. The proof presented here is almost identical to the original.

Lemma B.2

Let A 1 and A 2 be two affine subspaces in R d , with A 1 ⊂ A 2 . Let τ be a simplex in A 2 , and let v be a point in R d \ τ . Define J to be the join of τ and v. Then d(J , A 1 ) ≥ d(τ, A 1 ) d(v, A 2 ) L(J ) , [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF] where the distances between sets are defined as d(B, C) = inf x∈B,y∈C |x -y| and L(J ) denotes the longest edge length of an edge in J .

The notation is illustrated in Fig. 10.

Proof Let us suppose that (45) is false. Let J c be the truncated cone that consists of all half lines that start at a point of τ and pass through v. Then we may choose p J ∈ J c and a 1 ∈ A 1 so that

| p J -a 1 | = d(J c , A 1 );
by the definition of J c and the hypothesis we also have Now suppose that p J lies on the half line that starts at w ∈ τ and goes through v. Because τ ⊂ A 2 , we see that d(v, A 2 ) ≤ L e (J ). This means that [START_REF] Rennie | On Stirling numbers of the second kind[END_REF] gives that d(J c , A 1 ) < d(τ, A 1 ), so p J = w. We now immediately see that the line segment a 1 p J is orthogonal to the line that goes through w and v, which extends the half line we mentioned above. Let now be the line that goes through a 1 and w, and π (v) ∈ the point that is closest to v. It follows that π (v)w is perpendicular to . Because a 1 is nearer to p J than w, a 1 and π (v) are on the same side of w in . This means, because two of the angles are the same (and thus the third), that the triangles p J wa 1 and π (v)wv are similar. We now have that

d(J c , A 1 ) ≤ d(J , A 1 ) < d(τ, A 1 ) d(v, A 2 ) L(J ) . ( 46
)
d(J c , A 1 ) = |p J -a 1 | = |a 1 -w| • |v -π (v)| |v -w| ≥ d(τ, A 1 ) d(v, A 2 ) L(J ) ,
contradicting the hypothesis and thus proving the lemma.

Fig. 1

 1 Fig.[START_REF] Aamari | Estimating the reach of a manifold[END_REF] The two parts of the algorithm: part 1, where we perturb the vertices of the ambient triangulation, is depicted on top. Part 2, where the triangulation is constructed from the points of intersection of M and the edges, is depicted below

Definition 4 . 1 A

 41 monohedral 2 triangulation is called a Coxeter triangulation if all its d-simplices can be obtained by consecutive orthogonal reflections through facets of the d-simplices in the triangulation and the affine hulls of facets entirely consist of facets of d-simplices in the triangulation.

Corollary 4 .

 4 10 ([42, Corr. 3.15]) The number of (d + 1k)-faces of the permutahedron is k! S(d + 1, k).

) 1 ,

 1 Let τ d now be a top dimensional simplex in T that intersects M . Let τ d-n . . . , τ d-n j be the faces of τ d that intersect M . This means that dn + 1 barycentric coordinates with respect to τ d of any v(τ d-n i

1 ,

 1 3.1, near a simplex τ ∈ T , because it is small. So T M and N M near p are well approximated by T p M and N p M . -The angles between the normal space and faces τ d-n . . . , τ d-n j of τ of dimension dn that intersect M are bounded from below by Lemma 6.3. As a consequence, the orthogonal projection map π aff τ d-n k →N p M ≡ π τ d-n k from the affine hull aff τ d-n k

), p the point π - 1 v 1 v

 11 By definition of the projection map π v(τ d-nk lies in the translation of T p M that goes through u. Because, also by definition, π -

Lemma 7 . 4

 74 Suppose that τ d ∈ T and M ∩ τ d = ∅. Then, M ∩ τ d lies in the tubular neighbourhood of K ∩ τ d as defined in Sect. 7.1 (whose size is lower bounded byLemma 7.3).Proof Consider v(τ d ) ⊂ K ∩ τ d ,where we use the definition[START_REF] Coxeter | Discrete groups generated by reflections[END_REF], and choose an arbitrary

) 123 and thus α k = 2 k+1 ρ k 1 ck / 3 kFig. 9 Lemma B. 1

 1391 Fig. 9 An illustration of the notation of Lemma B.1

Fig. 10

 10 Fig. 10 Notation for the proof of Lemma B.2

  The 'normal' space of K at p, that is span N p (e i ) π M Simplices. We have tried to reserve τ for T or T and σ for K . However, for arbitrary simplices (such as in Appendix B) we use arbitrary choices. Subscripts are used for indices and superscripts for the dimension.

	Notation			Definition	Whitney's notation (if relevant)
	N v(τ ) (e i )		See (30)
	N p
				Closest point projection on M	π *
	π T p M π -1 p			Orthogonal projection on the tangent spaces T p M See Definition 3.3
	π aff τ d-n k π τ d-n k rch M	→N p M	=	The orthogonal projection map from the affine hull aff τ d-n to N p M i The reach to the manifold M
	ρ1			Volume fraction of the part of a ball inside a	ρ 1
	Notation ρ 1			Definition slab Lower bound on ρ1 , see (5)	Whitney's notation (if relevant)
	S A i aff B d (c, r ) T M , T p M T		Affine subspaces Slab between two hyperplanes intersected with a ball The affine hull A ball in R d with centre c and radius r , if we do not need to emphasise the centre or The tangent bundle and the tangent space at p The ambient Coxeter triangulation of type Ã	P, P and Q Q U r (c) P p Lis the ambient triangulation, but
	T			radius or they are to be determined, these are suppressed from the notation Perturbed ambient triangulation	is not a Coxeter triangulation L *
	τ, σ B T is t(σ ) Thickness of σ at most r 2 cL Perturbation radius of the vertices of T U (X , r ) A neighbourhood of radius r of a set X c Normalised perturbation radius d Ambient dimension (R d ) v i Vertices of T v * i Vertices of T	Same U r (X ) ρ ρ * m p i p * i
	d( • , • )			Euclidean distance between sets
	d M ( • , • ) δ		Distance on M Protection
				The sampling density as in an ( , μ)-net (the
				circumradius of the simplices in the
				Coxeter triangulation)
	K			Triangulation of M	K
	L( • )			Longest edge length	δ is the longest edge length of the
					ambient triangulation L
	L			L = L(T )
	L			L = L( T )
	λ			barycentric coordinates
	M			The manifold	M
	μ			Separation as in an ( , μ)-net (the shortest
				edge length in T for Coxeter
				triangulations)
	μ 0			The normalised separation, that is μ = μ 0
	n			Dimension of M	n
	N M , N p M		The normal bundle and normal space at p
	N ≤k			An upper bound on the total number of faces	Whitney does not distinguish
				of dimension less or equal to k that contain	dimensions and uses N as an
				a given vertex	upper bound (no value given)

p M (c, r ) A ball in T p M , using the same conventions as for

B d (c, r ) C(T p M , r 1 , r 2 )

Open cylinder given by all points that project orthogonally onto an open ball of radius r 1 in T p M and whose distance to this ball

IST Austria, Klosterneuburg, Austria

Discrete & Computational Geometry (2021) 66:386-434

The join of a simplex and a vertex is the convex hull of the vertices of the original simplex as well as the new vertex. Generally, the join of two subsets A, B ⊂ R d is defined as A * B = {λa + μb | a ∈ Ab ∈ B}, where λ, μ ∈ R, λ, μ ≥ 0, and μ + λ = 1, see for example[START_REF] Rourke | Introduction to Piecewise-Linear Topology[END_REF] Chap. 1].

A triangulation of R d is called monohedral if all its d-simplices are congruent.

In particular τ j ⊂ Ti .

⎞ ⎠ rch M .

d dn • (1 + 2 c) min alt τ d .

√ n + 1 • |λ -λ |.(36)

Acknowledgements

We thank the reviewers for their comments, which helped to improve the exposition.

Funding Information Open access funding provided by the Institute of Science and Technology (IST Austria).

Kenneth Clarkson

This work has been funded by the European Research Council under the European Union's ERC Grant Agreement Number 339025 GUDHI (Algorithmic Foundations of Geometric Understanding in Higher Dimensions). The third author also received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 754411.

Appendix B: Some Properties of Affine Spaces

In this appendix, we discuss two variants of lemmas from [START_REF] Whitney | Geometric Integration Theory[END_REF]App. II.14] that are essential in the building of the triangulation, see Sect. 6.1 in particular. Both lemmas are due to Whitney. However, in both cases, the statement is different, because we prefer to work directly with angles and use the thickness as our quality measure. In the

Appendix C: Remark on the C 1,1 case

We now first discuss a simpler version of Lemma 3.2 before going in to the C 1,1 setting. The result in this case is weaker, but can be easily extended to the C 1,1 setting as we shall see below. The following consequence of Lemma 3.1 is a stronger version of [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]Lem. 5.4]:

Proof Let q ∈ M be such that | p -q| ≤ r , then the differential of the projection map π T p M at q is non-degenerate, because, by Lemma 3.1, the angle ∠(T p M , T q M ) is less than π/2. Because M ∩ B( p, r ) is a topological ball of the right dimension by [18, Prop. 1], the result now follows.

Similarly to Lemma 3.1 we have for C 1,1 manifolds that:

,

This lemma gives us a corollary, which is the equivalent of Corollary C.1:

and p ∈ M , then for all r < (rch M )/3, the projection π T p M onto the tangent space T p M , restricted to M ∩ B( p, r ), is a diffeomorphism onto its image.

These are in fact all the fundamental results that are needed to be able to extend to the C 1,1 setting.

Assuming the manifold is C 1,1 would lead to minor changes in the calculations in the proof of Lemma 6.4 and would in theory influence the final conclusion in Sect. 7.2. However, because we have a significant margin in the difference between π/2 and the angle between N p and N p M , we would not need to change the constants in Sect. 7.2. Because we use the projection on the manifold, which is only Lipschitz, the map is a homeomorphism which is no longer piecewise smooth, but just Lipschitz. The rest of proofs hold verbatim.

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.