J.-D Boissonnat
email: jean-daniel.boissonnat@inria.fr

S Kachanovich
email: siargey.kachanovich@inria.fr

Mathijs Wintraecken
email: m.h.m.j.wintraecken@gmail.com

Triangulating submanifolds: An elementary and quantified version of Whitney's method

Keywords: Triangulations, Manifolds, Coxeter triangulations

We quantize 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.

Introduction

The question whether every C 1 manifold admits a triangulation was of great importance to topologists in the first half of the 20th century. This question was answered in the affirmative by Cairns [START_REF] Cairns | On the triangulation of regular loci[END_REF], see also Whitehead [START_REF] Whitehead | On C 1 -complexes[END_REF]. However the first proofs were complicated and not very geometric, let alone algorithmic. It was Whitney [START_REF] Whitney | Geometric Integration Theory[END_REF]ter IV], who eventually gave an insightful geometric constructive proof. Here, we will be reproving Theorem 12A of [START_REF] Whitney | Geometric Integration Theory[END_REF]Section IV.12], in a more quantitative/algorithmic fashion for C 2 manifolds: Theorem 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]Remark 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 generalize 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 implies that any C 1 manifolds 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 (see [START_REF] Whitney | Geometric Integration Theory[END_REF]Section 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] Cazals | Delaunay triangulation based surface reconstruction[END_REF][START_REF] Boissonnat | Meshing of surfaces[END_REF][START_REF] Cheng | Computer and information science series[END_REF][START_REF] Berger | State of the Art in Surface Reconstruction from Point Clouds[END_REF][START_REF] Dey | Curve and Surface Reconstruction; Algorithms with Mathematical Analysis[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 construct 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'll recall in detail in Section 4. 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 ndimensional convex polytope, more precisely the intersection is piecewise smoothly homeomorphic to a polytope. Proving this is the core of the homeomorphism proof in Section 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, which given criteria for a homeomorphism, but not for a 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 neigh-bourhood, 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] William | Marching cubes: A high resolution 3d surface construction algorithm[END_REF][START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF][START_REF] Doi | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[END_REF][START_REF] Allgower | An algorithm for piecewise-linear approximation of an implicitly defined manifold[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] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF][START_REF] Boissonnat | Isotopic implicit surface meshing[END_REF] in two and three dimensions) no topological correctness (homeomorphism) is proven 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]Theorem 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 with 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, such as in [START_REF] Lundell | The Topology of CW Complexes[END_REF][START_REF] Edelsbrunner | Triangulating topological spaces[END_REF].

The part of the algorithm described in this paper that constructs the triangulation (see part 2 of the algorithm in Section 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 Section 2. [START_REF] Aamari | Estimating the reach of a manifold[END_REF] is not yet complete at the moment of writing.

The algorithm and overview

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 Whitneys 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'll prove in Section 7). The second part is where the triangulation is constructed and is based on barycentric subdivision of polytopes. precise definition) and consists of two steps. In these two steps we have to carefully choose a significant 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 | Geometric and Topological Inference[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 Section 4 for the definitions) of the ambient triangulation as well as the longest edge length , separation, and dimension, see [START_REF] Brown | Sheaf-Theoretic Stratification Learning[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 found in [START_REF] Boissonnat | Manifold reconstruction using tangential Delaunay complexes[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 list, 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 ) ≥ 3 2 L, as well as the simplices with these vertices (see Case 1 of Section 5.2). For a vertex v i such that d(v i , M ) < 3 2 L (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 join 1 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 , such that it lies sufficiently far from the union of these spans, but also not to 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] Cheng | Dynamic skin triangulation[END_REF]. The existence of such a point can by proven by volume estimates and is shown in Lemma 21. The fact that such a perturbation ensures that the dn -1-skeleton is sufficiently far away from the manifold is non-trivial and proven in Lemma 22. We note that for a curve in 2 dimensions, as depicted in Figure 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 1 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]Chapter 1]. 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] Tamal | 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 [START_REF] Behr | Simplex spacetime meshes in finite element simulations[END_REF] below. The number of simplices in turn is bounded by (d -n) d d (d-n) , see Lemma 17. 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.

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] Tamal | Delaunay meshing of piecewise smooth complexes without expensive predicates[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'll see in Section 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 Figure 2. Because of 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 solutions to numerical partial differential equations, in particular space time methods [START_REF] Shewchuk | What is a good linear finite element? -interpolation, conditioning, anisotropy, and quality measures[END_REF][START_REF] Behr | Simplex spacetime meshes in finite element simulations[END_REF][START_REF] Max Von Danwitz | Simplex space-time meshes in compressible flow simulations[END_REF]. This also serves as a first step in generalizing the work on the triangulation of general stratifolds in three dimensions [START_REF] Oudot | Meshing Volumes Bounded by Smooth Surfaces[END_REF][START_REF] Rineau | Meshing Volumes bounded by Piecewise Smooth Surfaces[END_REF][START_REF] Dey | Localized delaunay refinement for piecewise-smooth complexes[END_REF][START_REF] Tamal | Localized delaunay refinement for volumes[END_REF][START_REF] Tamal | Delaunay meshing of piecewise smooth complexes without expensive predicates[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].

Fig. 2 The same triangulation as depicted in Figure 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 Section 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 Section 3 we recall some results on the geometry of manifolds of poistive reach. Coxeter triangulations, Delaunay protection, and the combinatorial stability of a triangulation under perturbations is the topic of Section 4.

In Section 5, we both give the details of the perturbation of the vertices and some geometric consequences for the triangulation. In Section 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 Section 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, NM is the normal bundle and N p M a normal space. Distances on the 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, Smale and Weinberger [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]:

Lemma 2 (Lemma 6 and Corollary 3 of [START_REF] Boissonnat | The Reach, Metric Distortion, Geodesic Convexity and the Variation of Tangent Spaces[END_REF]) 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| 2rch(M )
.

In Lemma 3 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, Smale, and Weinberger [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].

p T p M M C(T p M , r 1 , r 2 ) r 1 r 2
Fig. 3 The cylinder C(T p M , r 1 , r 2 ), with the manifold and tangent space.

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 Figure 3 for an illustration. We now have: Lemma 3 Suppose that M is C 2 and p ∈ M , then for all r < rch(M ), the projection π 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 2, we'll 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 minimizing geodesic γ on M with length parametrized 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 2) as long as < 1 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 )). (2)
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 4

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 Lemma IV.8a of [START_REF] Whitney | Geometric Integration Theory[END_REF]: Corollary 5 Suppose that M is C 2 and p inM , 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. (3)
Proof Lemma 3 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 6 (Distance to Manifold, Lemma 11 of [START_REF] Boissonnat | The Reach, Metric Distortion, Geodesic Convexity and the Variation of Tangent Spaces[END_REF]) 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 4. 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 4 families for which there is one member for (almost) every dimension d. We will focus on one such family, Ãd , which is Delaunay protected.

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 8 Let P = {(x i ) ∈ R d+1 | ∑ i x i = 0} and consider the d-simplex with ver- tices 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 9

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 dsimplices 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]Section 6]. We shall denote the triangulation of this type by T .

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 Ã2 Ã3

Fig. 4 The vertex sets of the Coxeter triangulations in dimensions two and three are the triangular lattice and the body-centred cubic lattice, respectively.

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 [START_REF] Whitney | Geometric Integration Theory[END_REF]Appendix Section II.16]. The argument for stability of triangulations for à type Coxeter triangulations is much simpler, because it is a Delaunay triangulation and δ -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 normalized 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. Theorem 4.14 of [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] immediately gives: Corollary 10 The triangulation T is (combinatorially) stable under a cL-perturbation as long as

cL ≤ t(T )µ 0 18d δ . (4)
We claim the following concerning the behaviour of c Lemma 11

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

Proof Choudhary et al. [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF]Appendix B] provide explicit values of all the quantities mentioned in Corollary 10 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 ε = d(d + 2) 12(d + 1) δ (σ ) = √ d 2 + 2d + 24 - √ d 2 + 2d 12(d + 1) . (5)
The value of µ easily follows from the general expression of edge lengths (see [25, Appendix B, Ãd , item 5]) and is equal to µ = d d+1 . From ( 5), we get that

µ 0 = µ ε = 12 d+2 . The bound in (4) is therefore c ≤ t(σ )µ 0 18d δ L =      √ 2 √ d 2 +2d+24- √ d 2 +2d 9d 3/2 √ d+2(d+1) 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 + 2(d + 1) ∼ √ 32 3d 4 ,
where we used that

√ 1 + x 1 + 1 2 x 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 Section IV.14 and Section IV.15 of [START_REF] Whitney | Geometric Integration Theory[END_REF].

We remind ourselves of the following:

Lemma 12 (Thickness under distortion [34, Lemma 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(σ ).
We can now state a variation of Whitney's angle bound result, see [START_REF] Whitney | Geometric Integration Theory[END_REF]Section IV.15].

Lemma 13 (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 on 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 Section 5. We will now give an explicit bound on this number.

In general the (dk)-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 14 (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 ourselves of the following definition, see [START_REF] Abramowitz | 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]: 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

Definition 15 Let S(
S(d, k) = 1 k! k ∑ j=0 (-1) j k j (k -j) d .
Corollary 16 (Corollary 3.15 of [START_REF] Maes | On the permutahedron and the quadratic placement problem[END_REF]) The number of

(d + 1 -k)-faces of the per- mutahedron is k!S(d + 1, k).
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

N ≤k = 2 + k ∑ j=1 j!S(d + 1, j), (6)
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 ( 20)), and have a strict inequality. We now claim the following:

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

≤ j ≤ d -1 1 2 ( j 2 + j + 2) j d-j-1 -1 ≤ S(d, j) ≤ 1 2 d j j d-j .
Furthermore, Stirlings 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 Section 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 quantized. 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 Section 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 Section 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 the constants. This subsection corresponds to Section IV.18 of [START_REF] Whitney | Geometric Integration Theory[END_REF].

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] Behr | Simplex spacetime meshes in finite element simulations[END_REF].

A precise bound on ρ1 can be given, see Remark 1 below. We will use an easier bound ρ 1 , at the cost of weakening the result:

Lemma 18 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. ( 7
)
Note that

ρ 1 ∼ 1 √ dN ≤d-n-1 .
Proof We can bound the volume of the slab S by the 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 rvol(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 [33, page 622]. Note that the inequality is strict because ρ 1 > 0. We see that therefore ρ 1 may be chosen to be

ρ 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. (7)
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 √ dN ≤d-n-1 .
Remark 1 Because of symmetry, the largest volume S can attain is when both delimiting hyperplanes are equidistant to the center 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 δ 18dL , 1 24 t(T ) 2 . ( 8
)
The reasons for this particular choice will be discussed after [START_REF] Brown | Sheaf-Theoretic Stratification Learning[END_REF] below. We stress that ( 8) 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 . (9)
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 , (10)
that is α k = 2 k+1 3 k ρ k 1 ck
. 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 , (11)
because c ≤ 1 24 , as we have seen in ( 9). ρ 1 is also very small, as a direct consequence of Lemma 18. Furthermore we notice that α k < α k-1 . To make sure the formulae do 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 ). (12)
Note that ζ depends on both the ambient and intrinsic dimension and the perturbation parameter c. Because c ≤ 1 24 t(T ) 2 and t(T ) ≤ 1 we see that ζ ≤ 1. We set the coarseness of the ambient triangulation by demanding that L satisfies

  1 -1 - 6L(T ) rch(M ) 2   rch(M ) = (α d-n ) 4+2n 6(n + 1) 2 ζ 2n L, (13)
or equivalently

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

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

Remark 2

We have to choose the right hand side in (13) 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 28.

(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 19 ((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, (16)
for all faces τ k in T , with k ≤ dn -1, and

L < 13 12 L (17)
t( T ) ≥ 4 5 √ d 1 - 8 c t(T ) 2 t(T ). (18)

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, Section 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 δ , 1 24 t(T ) 2 L , (19)
using the notation of Section 4. With this bound we have that ( 17) 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 20 T has the same combinatorial structure as T . Moreover, (18) is satisfied.

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

By the third condition of [START_REF] Brown | Sheaf-Theoretic Stratification Learning[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 ), (18)
as a consequence of Lemma 12, the fact that if you perturb the vertices by cL the edge lengths are perturbed 2 c (that is 2 c = c 0 ) and the fact that if σ ⊂ τ, then t(σ ) ≥ t(τ). So we have established [START_REF] Boissonnat | The Reach, Metric Distortion, Geodesic Convexity and the Variation of Tangent Spaces[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, ( 16) is satisfied. We distinguish two cases:

Case 1 (d(v i , M ) ≥ 3 2 L).
In this case we choose ṽi = v i . The inequality ( 16) 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] Berger | State of the Art in Surface Reconstruction from Point Clouds[END_REF] we have that ( 1 2 -2 c)L > 5 12 L. This means that d(τ k , M ) > 5 12 L for any simplex in the star of ṽi = v i . This lower bound is much larger than

α k L < 1 18 k L. Case 2 (d(v i , M ) < 3 2 L). Let p be a point in M such that d(v i , p) < 3 2 L. Let τ 0 (= / 0), τ 1 , . . . , τ ν (20)
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] Behr | Simplex spacetime meshes in finite element simulations[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 21 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, ( 21
)
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 Section 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 = vol(B(v i , cL)) 1 - ν + 1 2N ≤d-n-1 ≥ 1 2 vol(B(v i , cL)),
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, (22)
for all 1 ≤ j ≤ ν.

The following lemma completes Case 2:

Lemma 22 The triangulation T is (d -n -1)-skeleton safe, in particular (16) 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 on 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 6, 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 - (α d-n ) 4+2n 6(n + 1) 2 ζ 2n L (using (13)) > α k-1 L - 1 3 α d-n L (using (15)) ≥ 2 3 α k-1 L. (because α k-1 > α k ) (23)
Because d(v i , p) < 3 2 L and L < L + 2 cL, and c < 1 24 , see (9), we have the very coarse bound that

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

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

d(τ j , T p M ) > 2 3 α k-1 L. ( 25
)
Lemma 33, 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 [START_REF] Choudhary | Coxeter triangulations have good quality[END_REF] and ( 21))

d(τ j , T p M ) ≥ d(τ j , T p M )d(v i , span(τ j , T p M )) L + 2 cL > 2 3 α k-1 L ρ 1 cL L + 2 cL (thanks to
> 2 3 L 4 3 L α k-1 ρ 1 cL (because c ≤ 1 24 ) = 4 3 α k L.
(using the relation [START_REF] Boissonnat | Meshing of surfaces[END_REF]

for α k ) (26)
Similarly to [START_REF] Cheng | Computer and information science series[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] Coxeter | Discrete groups generated by reflections[END_REF], to the distance to the manifold as follows. Because of Corollary 5 we can localize the results and Lemma 6 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 13) and ( 15))

d(τ j , M ) > 4 3 α k L -   1 -1 - 6L rch(M ) 2   rch(M ) > 4 3 α k L - 1 3 α d-n L (using (
≥ α k L. (because α k ≥ α d-n if k ≤ d -n -1 by (10))
This completes the proof for the case where j = 0 or τ j is non-empty. For j = 0, [START_REF] Cheng | Dynamic skin triangulation[END_REF] and Lemma 6 yield

d(τ j , M ) > ρ 1 cL -   1 -1 - 6L rch(M ) 2   rch(M ).
We simplify 13) and ( 15))

d(τ j , M ) > ρ 1 cL -   1 -1 - 6L rch(M ) 2   rch(M ). > ρ 1 cL - 1 3 α d-n L (using (
> α 1 L.
(by definition of ( 10))

We emphasize 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 24. Here we again rely on Appendix B. The triangulation K of M is defined in Section 6.2.

The geometry of the intersection of simplices in T and M

In this section, we discuss the geometry of simplices in T in relation to M . We follow [52, Section 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 23 Let p ∈ M and suppose that

τ k ∈ T , with k ≤ d -n -1, 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 23 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

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 6) > α k L - 1 3 α d-n L, (d(τ k , M ) > α k L and (15)) > 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 24 Suppose that M intersects τ k ∈ T . Let p ∈ M , such that τ k ⊂ B(p, 6(L), then T p M intersects τ k .

Proof Let p ∈ M ∩τ k . Lemma 3 (and ( 13), ( 15)) gives us that p ∈ π -1 p (B T p M (p, 6L)), where we use the notation of Definition 4. Lemma 6 gives that

d( p, T p M ) ≤   1 -1 - 6L rch(M ) 2   rch(M ) < 1 3 α d-n L.
Let τ ⊂ τ k be the face of smallest dimension such that d( τ, T p M ) ≤ 2 3 α d-n L. This face exists thanks to the triangle inequality. By Lemma 23 we have the dim( τ) ≥ dn. Lemma 32 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 3 α k L 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 25 Suppose that M intersects τ k ∈ T and τ k has dimension d -n, that is k = d -n. Let p ∈ M , such that τ k ⊂ B(p, 6L), then sin ∠(aff(τ k ), T p M ) ≥ 2d(T p M , ∂ τ k ) L + 2 cL ≥ 4 3 α k L L + 2 cL ≥ 16 13 α k .

Proof

This is an immediate consequence of Lemma 32,[START_REF] Brown | Sheaf-Theoretic Stratification Learning[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 Section IV.19 in the following lemma. The proof differs from the original by Whitney. Proof Let τ be a face of smallest dimension of τ k such that d( τ, T p M ) ≤ 2 3 α n L. Now Lemma 32 and Lemma 24 give that τ and T p M have a unique point p in common and the dimension of τ is dn.

Thanks to Lemma 3, 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 really 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 parametrization of T p M . Similarly, we can define G : R d-n → R d to be a linear (orthonormal) parametrization 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 25 we have that sin ∠(aff( τ), T p M ) ≥ 16 13 α d-n .
Lemma 2, and [START_REF] Boissonnat | Local Criteria for Triangulation of Manifolds[END_REF] give that for any q ∈ B(p, 6L)

sin ∠(T p M , T q M ) 2 ≤ 6L 2rch(M ) ≤ 6 2 • (α d-n ) 2 54 = 1 18 (α d-n ) 2 .
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, 6, and 25, to bound | ppλ |. The distance from the manifold to the tangent space is bounded from above by

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

| p -pλ | ≤ 1 3 (α d-n ) 2 L 16 13 α d-n ≤ 1 3 α d-n L.
This distance bound is smaller than the distance bound of p to the boundary of τ, due to Lemma 23. 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 27 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 26, we see that T p M intersects some (dn)-face of τk-1 . Thanks to Lemma 25, the angle between this (dn)-face and T p M is bounded from below. Due to Lemma 23, 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 26 finally says that M intersects τk , because T p M does.

6.2

The triangulation of M : The complex K

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

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]Theorem 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 )} (27)
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 26. 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 27. We now define v(τ) as follows:

v(τ) = 1 j v(τ d-n 1 ) + • • • + v(τ d-n j ) . ( 28
)
Remark 3 We stress that thanks to Lemma 26, 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 28 Let σ n be a top dimensional simplex as defined in [START_REF] Dey | Curve and Surface Reconstruction; Algorithms with Mathematical Analysis[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 | Delaunay stability via perturbations[END_REF].

Proof This inequality relies on estimates on the barycentric coordinates, and Lemma 12. 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 23, v(τ d-n i ) lies at least a distance 2 3 α d-n-1 L 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 ( 29), 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 . ( 29
τ d of v(τ k ) = 1 j v(τ d-n 1 ) + • • • + v(τ d-n j ) , for k > d -n, we have that -k + 1 of the coordinates λ l satisfy λ l > 2 3 1 j 1 1 + 2 c α d-n-1 .
-The other dk coordinates are zero.

Note that j ≤ d d-n . This means that

d(v(τ k ), ∂ τ k ) ≥ 2 3 d d-n • (1 + 2 c) α d-n-1 min alt(τ d ).
We now have

d(v(τ k ), ∂ τ k ) ≥ 2 3 d d-n • (1 + 2 c) α d-n-1 t( T ) L (by definition of the thickness) ≥ 2 3 d d-n • (1 + 2 c) α d-n-1 4 5 √ d 1 - 8 c t(T ) 2 t(T ) L (by the estimate (18)) ≥ 8 15 √ d d d-n • (1 + 2 c) 1 - 8 c t(T ) 2 α d-n-1 t(T ) L ( 30
)
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 the min alt(σ ) ≥ h d min L(σ ) , as can be established by simple trigonometric arguments, as illustrated in Figure 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 , ( 31
)
where we indicated the dimensions explicitly. Plugging (30) (using the definition of the simplex σ n as in ( 27)) into [START_REF] Dey | Localized delaunay refinement for piecewise-smooth complexes[END_REF] gives that min alt(σ n ) is lower bounded as follows:

min alt(σ n ) > 8 15 √ d d d-n • (1 + 2 c) 1 - 8 c t(T ) 2 α d-n-1 t(T ) n L = ζ n (α d-n-1 ) n 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 Section 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 Section 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 homeomorphism proof. The homeomorphism 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 xp ∈ 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 the 'tubular neighbourhood' of K. We use two results from the previous sections:

-The normal space is almost constant, see Lemma 2, near a simplex τ ∈ T , because it is small. So T M and NM near p are well approximated by T p M and N p M . -The angles between the normal space and faces τ d-n , 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.

Fig. 6 The tubular neighbourhood.

We start, as mentioned, with the vertices that are associated to a simplex τ = τ d-n ∈ T of dimension dn. We stress that there is but one face of dimension dn of τ, 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 Section 6.2. Let p ∈ M be such that τ ⊂ B(p, 6L). For now p is arbitrary, we'll 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) . (32)
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). (33)
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 4 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 .

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 29 Suppose that p ∈ M , τ d ⊂ B(p, 6L) and σ n ∈ K 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 ) ≤ (α d-n ) 4+n 6(n + 1) ζ n .
Proof By Lemma 13, 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 28 gives us the following bound

min alt(σ n ) > (α d-n-1 ) n ζ n L.
Finally, d max is bounded thanks to [START_REF] Boissonnat | Geometric and Topological Inference[END_REF]. Combining these results yields

sin ∠(aff(σ n ), T p M ) ≤ (n + 1) (α d-n ) 4+2n 6(n+1) 2 ζ 2n L (α d-n-1 ) n ζ n L ≤ (α d-n ) 4+n 6(n + 1) ζ n .
(because α d-n-1 < α d-n , and L ≥ L because there are unperturbed simplices in T )

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

Lemma 30 Let p, q ∈ σ n , with barycentric coordinates λ = (λ 0 , . . . , λ n ), λ = (λ 0 , . . . , λ n ) respectively. Suppose that N p and N q be as defined in Section 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 (α d-n ) 4 n + 1 ζ n (α d-n-1 ) n L. (34)
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 (34).

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 do not intersect 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 naive bounds to fairly sharpe 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 25 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 |π -

| ≤ 13 16α d-n 1 v(τ d-n k ),p (u)
. Together with the triangle inequality this yields that

|N τ,p (u)|, |N p(u )| ≤ 13 16 α d-n , (35)
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 13 α d-n . (36)
Thus ∠( span

u∈N p M (N τ,p (u)), N p M ), ∠(N p, N π M ( p) M ) ≤ arccos 16 13 α d-n . (37)
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 2, and [START_REF] Boissonnat | Local Criteria for Triangulation of Manifolds[END_REF]

to see that sin ∠(N p M , N q M ) 2 ≤ 2L rch(M ) < (α d-n ) 4+2n 29(n + 1) 2 ζ n .
Using [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF] and 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 13 α d-n + 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: With this notation, we can now derive the following bound:

-(N τ d-n 0 ,p (u) . . . N τ d n ,p (
|λ 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 |λ -λ | ≤ (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 the Frobenius norm)

≤ 13 √ n + 1 16 α d-n |λ -λ |.
(by ( 35))

≤ 13 √ n+1 16 α d-n |λ -λ | 1 ≥ 1 θ 0 θ 0 φ 0 (38)
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 10 α d-n ). The length or bound on the length of two of the edges is also indicated in the figure .   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 Figure 7. We apply the sine rule to this triangle, using [START_REF] Federer | Curvature measures[END_REF], and [START_REF] Kachanovich | Meshing submanifolds using Coxeter 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 1 = 13 16 α d-n 2 √ n + 1|λ -λ |. (39)
Note that this can be tightened a fair bit at the cost of complicating the bound.

We conclude from (39) that sin ∠(span

u (N p,p (u)), span u (N q,p (u))) ≤ 13 16 α d-n 2 √ n + 1|λ -λ |. (40)
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)), span u∈N q M (N p,q (u)) based on the distance between the points p and q in M . We use that p and q are such that τ ⊂ B(p, 6L), B(q, 6L) so that the conditions of Lemma 25 hold. This also means that p and q are close, which means that 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)) (41)
To bound this angle, we look at the individual terms in [START_REF] Doi | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[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 Figure 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 ) .
The distance from π -1 v(τ d-n k ),q (u + wq,p ) to the translation of T p M that goes through u ), we have that, .

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 ) = 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 ) . (42)
By Lemma 2 sin ∠(T p M , T q M ) 2 ≤ |p -q| 2rch(M ) , so that, using 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 4rch(M ) 2 1 -|p-q| 2 2rch(M ) 2 .
This means that (42) 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 4rch(M ) 2 1 -|p-q| 2 2rch(M ) 2
By the triangle inequality (on the terms in the sum in the definition [START_REF] Doi | An efficient method of triangulating equi-valued surfaces by using tetrahedral cells[END_REF]), this gives

|N p,p (u)) -N p,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 4rch(M ) 2 1 -|p-q| 2 2rch(M ) 2
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))) ≤ 1 -16 13 α d-n 2 16 13 α d-n 2 |p-q| rch(M ) 1 -|p-q| 2 4rch(M ) 2 1 -|p-q| 2 2rch(M ) 2
Because p and q are very chose (in fact they lie in an L neighbourhood of τ, which is small due to ( 15))

1 -|p-q| 2 4rch(M ) 2 1 -|p-q| 2 2rch(M ) 2 ≤ 2,
as can verified using that

√ 1-x 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] Maes | On the permutahedron and the quadratic placement problem[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 ) . (43)
Step 3: Combining into a bound on ∠(N p, N q). Combining ( 43), [START_REF] Lundell | The Topology of CW Complexes[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 ) ≤ 1 2 rch(M ), 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)) = √ 1 -x 2 y + 1 -y 2 x ≤ 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 ( 15) and ( 9),

13 16α d-n 2 2|λ -λ | L rch(M ) ≤ 13 16α d-n 2 4|λ -λ | (α d-n ) 4+2n 54(n + 1) 2 ≤ 6 100 |λ -λ |(α d-n ) 2 .
Thanks to [START_REF] Boissonnat | The stability of Delaunay triangulations[END_REF], the first term in the following sum is by far the largest

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

2 √ n + 1|λ -λ |, (44)
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|. We remind ourselves 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 following 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 Lemma 5.12 of [START_REF] Boissonnat | Geometric and Topological Inference[END_REF] 3. Lemma 28 to bound min alt(σ ) 4. Equation [START_REF] Plantinga | Isotopic approximation of implicit curves and surfaces[END_REF], which gives a bound on the angle ∠ px q, namely φ 0 5. Lemma 29 gives that

∠(aff(σ n ) ⊥ , N p M ) ≤ arcsin (α d-n ) 4+n 6(n + 1) ζ n ≤ arcsin (α d-n ) 4 6 , (45)
where aff(σ n ) ⊥ denotes the space perpendicular to aff(σ n ). Because e j ∈ N p M , combining this with Lemma 25 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 13 α d-n + arcsin (α d-n ) 4 6 . (46)
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 remind ourselves of the trigonometric identity cos (arccos(a) + arcsin(

b)) = a 1 -b 2 -b 1 -a 2 .
Using (46) now gives sin ∠(N p (e j ), aff(σ n )) ≥ cos arccos 16 13

α d-n + arcsin (α d-n ) 4 6 = 16 13 α d-n 1 - (α d-n ) 8 6 2 - (α d-n ) 4 6 1 - 16 13 α d-n 2 ≥ 16 13 α d-n - 16 13 α d-n (α d-n ) 8 6 2 - (α d-n ) 4 6 ≥ 16 13 α d-n - (α d-n ) 18 3 (because α k ≤ 1 18 k by (11)) ≥ α d-n . (47)
This completes the fifth point.

The considerations we summed up yield:

| px|, | qx| ≥ (|λ -λ | min alt(σ )/ √ n)α d-n 10 15 α d-n 2 √ n + 1|λ -λ | ≥ min alt(σ )α d-n 15 10 α d-n 2 n + 1 ≥ α d-n 15 13 α d-n 2 n + 1 (ζ α d-n-1 ) n L
Using [START_REF] Rourke | Introduction to piecewise-linear topology[END_REF] again yields that the distance from x to aff(σ n ) is bounded from below by

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

Lemma 31 Suppose τ d ∈ T and that M ∩ τ d = / 0. Then, M ∩ τ d lies in the tubular neighbourhood of K ∩ τ d as defined in Section 7.1 (whose size is lower bounded by Lemma 30).

Proof Consider v(τ d ) ⊂ K ∩ τ d , where we use the definition [START_REF] Tamal | Delaunay meshing of piecewise smooth complexes without expensive predicates[END_REF] and choose an arbitrary n dimensional simplex

σ n ⊂ K ∩ τ d . Note that v(τ d ) ∈ K ∩ τ d . Thanks to Lemma 29, sin ∠(aff(σ n ), T v M ) ≤ (α d-n ) 4+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 (α d-n ) 4+n 6(n + 1) ζ n L,
where d H denotes the Hausdorff distance. Because of Lemma 6 and ( 13), we have that

d H (T v M ∩ B(v, 2L), π -1 v (B T v M (v, 2L))) ≤ (α d-n ) 4+2n 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 (α d-n ) 4+n 6(n + 1) ζ n L + (α d-n ) 4+2n 6(n + 1) 2 ζ 2n L (by the triangle inequality) ≤ (α d-n ) 4+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 30, that is

(α d-n ) 4+n ζ n n + 1 L ≤ 15 13 2 (α d-n ) 4 n + 1 ζ n (α d-n-1 ) n L, (48)
M ∩ τ is contained in this neighbourhood of K.

Having established that M lies in the tubular neighbourhood around K we can sensibly 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 Section 7.1) is transversal (Because π/2 minus the angle between N p and N p M , see [START_REF] Hirsch | Differential Topology[END_REF], is much bigger than the variation of the tangent/normal space as bounded by Lemma 2 and ( 14)) 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. We emphasize that along the way we have also given bounds on the Hausdorff distance between M and K, Lemmas 31 and 30 the quality of simplices, see Lemma 28 the variation of the tangent spaces, see Lemma 2 and equations [START_REF] Hirsch | Differential Topology[END_REF], and [START_REF] Boissonnat | Isotopic implicit surface meshing[END_REF].

A Notation

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

Notation definition

Whitney's notation (if relevant) A i

Affine subspaces P, P and Q aff

The affine hull B d (c, r)

A ball with in R d of dimension with centre c and radius r, if we do not need to emphasize the centre or radius or they are to be determined, these are suppressed from the notation 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 [START_REF] Behr | Simplex spacetime meshes in finite element simulations[END_REF]. α 1 and α k have been defined by a recursion relation as follows

α 1 = 4 3 ρ 1 c 2 3 α k-1 cρ 1 = α k , (10)
We often use

L rch(M ) < (α d-n ) 4+2n 54(n + 1) 2 ζ n < (α d-n ) 2 54 , (α d-n ) 4+2n 6(n + 1) 2 ζ 2n < (α d-n ) 2 3 ≤ α d-n 3 . (15)
The normalized perturbation radius c satisfies

|v i -ṽi | ≤ cL = min t(T )µ 0 18d δ , 1 24 t(T ) 2 L , (19)
from which it follows that c ≤ 1 24 .

B Some properties of affine spaces

In this appendix, we discuss two variants of lemmas from Appendix Section II.14 of [START_REF] Whitney | Geometric Integration Theory[END_REF], that are essential in the building of the triangulation, see Section 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 first case, the proof we provide differs significantly from the original. The first lemma will allow us to prove that if T p M intersects a simplex τ ∈ T and p and τ are not too far from each other then M intersects τ and vice versa. The second result is essential in proving that the perturbation of the vertices as described in Section 2.1 part 1, gives a triangulation for which the low dimensional simplices are sufficiently far away from the manifold. We start with a variation on Lemma 14a of Appendix Section II.14 of [START_REF] Whitney | Geometric Integration Theory[END_REF].

Lemma 32 Let σ be an s-simplex and A 0 an affine n-dimensional subspace in R d . Assume that s + n ≥ d and d(A 0 , σ ) < d(A 0 , ∂ σ ).

Then s + n = d, A 0 intersects σ in a single point and sin ∠(aff(σ ), A 0 ) ≥ 2d(A 0 , ∂ σ )/L(σ ). C Remark on the C 1,1 case

We now first discuss a simpler version of Lemma 3 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 2 of is a stronger version of Lemma 5.4 of [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF]:

Corollary 34 Suppose M is C 2 and p ∈ M , then for all 0 < r < rch(M ) √ 2 , the projection π TpM onto the tangent space T p M , restricted to M ∩ B(p, r) is a diffeomorphism onto its image.

Proof Let q ∈ M such that |p -q| ≤ r, then the differential of the projection map π TpM at q is nondegenerate, because by Lemma 2, 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 Proposition 1 of [START_REF] Boissonnat | The Reach, Metric Distortion, Geodesic Convexity and the Variation of Tangent Spaces[END_REF], the result now follows.

Similarly to Lemma 2 we have for C 1,1 manifolds that: Lemma 35 (Theorem 3 of [START_REF] Boissonnat | The Reach, Metric Distortion, Geodesic Convexity and the Variation of Tangent Spaces[END_REF]) Now suppose that M has positive reach, that is M is at least C 1,1 , and let |p -q| ≤ rch(M )/3, then

sin ∠(T p M , T q M ) 2 ≤ 1 - √ 1 -α 2 α 2 4 -( α 2 2 + 1 - √ 1 -α 2 ) 2
, where α = |p -q|/rch(M ).

This lemma gives us a corollary, which is the equivalent of Corollary 34:

Corollary 36 Suppose M is C 1,1 and p ∈ M , then for all r < rch(M )

3

, the projection π TpM 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 26 and would in theory influence the final conclusion in Section 7.2. However, because we have a significant margin in the difference between π/2 minus the angle between N p and N p M we would not need to change the constants in Section 7.2. Because we use the projection on the manifold, which is only Lipschitz, the map is homeomorphism is no longer piecewise smooth, but just Lipschitz. The rest of proofs hold verbatim.

Fig. 1

 1 Fig.1The 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 7 A

 7 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.

Lemma 26

 26 If p ∈ M , τ k ∈ T and τ 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 .

) 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

Fig. 5

 5 Fig.5Both triangles are right angled. We stress that the projection of p onto aff(σ ) may be quite a distance from σ itself.

1 ,

 1 . . . , τ d-n j of τ of dimension dn that intersect M are bounded from below by Lemma 25. 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 ) to N p M is a (linear) bijection, for any p that is sufficiently close to τ d-n k

  u)) denotes the matrix whose columns are N τ d-n 0 ,p (u), . . . , N τ d n ,p (u) -• 2 denotes the operator 2-norm, -• F the Frobenius norm.

1 v 1 v

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

Fig. 8

 8 Fig. 8 Note that the two dimensional nature of the figure is slightly misleading.

  wq,p ) = π -1 N q M →N p M (u) are both contained in aff(τ d-n k

Lemma 25 gives

 25 us that sin ∠(aff(τ d-n k ), T p M ), sin ∠(aff(τ d-n k ), T q M ) ≥16 13 α d-n , so that cos ∠(aff(τ d-n k ), T q M ) ≥16 13 α d-n and tan ∠(aff(τ d-

1 pk

 1 U r (c) B TpM (c, r)A ball in T p M , using the same conventions as forB 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 is at most r 2 cL Perturbation radius of the vertices of T ρ cNormalized perturbation radiusρ * d Ambient dimension (R d ) m 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) Separation as in an (ε, µ)-net (the shortest edge length in T for Coxeter triangulations) µ 0The normalized separation, that isµ = µ 0 ε n Dimension of M n NM , N p MThe normal bundle and normal space at p N ≤k An upper bound on the total number of faces of dimension less or equal to k that contain a given vertex. Whitney does not distinguish dimensions and uses N as an upper bound. (no value given) N v(τ) (e i ) See (32) N p The 'normal' space of K at p, that is span(N p(e i )) π M Closest point projection on M π * π TpM Orthogonal projection on the tangent spaces T p M π -See The orthogonal projection map from the affine hull aff(τ d-n i ) to N p M . rch(M ) The reach to the manifold M ρ1 Volume fraction of the part of a ball inside a slab ρ 1 ρ 1 Lower bound on ρ1 , see (7) S Slab between two hyperplanes intersected with a ball Q T M , T p M The tangent bundle and the tangent space at p P p T The ambient Coxeter triangulation of type à L is the ambient triangulation, but is not a Coxeter triangulation T Perturbed ambient triangulation L * τ, σ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.Same t(σ )Thickness of σ U(X, r)A neighbourhood of radius r of a set X U r

and thus α k = 2 k+1 3 k ρ k 1 ck.

 1 In particular, we have the bound ) 2 t(T ).

Fig. 9

 9 Fig. 9 An illustration of the notation of Lemma 32.

1 Fig. 10

 110 Fig. 10 Notation for the proof of Lemma 33.

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

In particular τ j ⊂ Ti .

Acknowledgments

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

This work has been funded by the European Research Council under the European Unions 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.

J.-D. Boissonnat , INRIA Sophia-Antipolis

Proof of Lemma 32 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 . Clearly translation leaves distances invariant, so

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'll now show that p ∈ σ . I'll 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

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 of Appendix Section II.14 of [START_REF] Whitney | Geometric Integration Theory[END_REF]. The proof presented here is almost identical to the original. Lemma 33 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

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.

Proof of Lemma 33 Let us suppose that ( 49) 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

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] Wendel | Note on the gamma function[END_REF] gives that d(J c , A 1 ) < d(τ, A 1 ), so that 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

contradicting the hypothesis and thus proving the lemma.