The Topological Correctness of PL Approximations of Isomanifolds

Isomanifolds are the generalization of isosurfaces to arbitrary dimension and codimension, i.e. manifolds defined as the zero set of some multivariate multivalued smooth function f:Rd→Rd-n\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f: {\mathbb {R}}^d\rightarrow {\mathbb {R}}^{d-n}$$\end{document}. A natural (and efficient) way to approximate an isomanifold is to consider its piecewise-linear (PL) approximation based on a triangulation T\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}$$\end{document} of the ambient space Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^d$$\end{document}. In this paper, we give conditions under which the PL approximation of an isomanifold is topologically equivalent to the isomanifold. The conditions are easy to satisfy in the sense that they can always be met by taking a sufficiently fine and thick triangulation T\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {T}$$\end{document}. This contrasts with previous results on the triangulation of manifolds where, in arbitrary dimensions, delicate perturbations are needed to guarantee topological correctness, which leads to strong limitations in practice. We further give a bound on the Fréchet distance between the original isomanifold and its PL approximation. Finally, we show analogous results for the PL approximation of an isomanifold with boundary.


Introduction
Isosurfacing Given a surface represented in R 3 as the zero set of a function f : R 3 → R, the goal of isosurfacing is to find a piecewise-linear (PL) approximation of the surface. This question naturally extends to higher dimensions and codimensions, in which case the generalized surface is called an isomanifold. Isosurfaces play a crucial role in medical imaging, computer graphics and geometry processing [1]. Higher-dimensional isomanifolds are also of fundamental importance in many fields such as statistics [2], dynamical systems [3], econometrics or mechanics [1].

Marching Algorithms
The standard algorithmic solution to the isosurfacing problem is to use some marching algorithm. This approach was initiated by Lorensen and Cline, with their marching cube algorithm [4]. Many variants of the algorithm have been introduced, see, for example, [5][6][7][8], and the overview [1]. The approach, however, is always the same: one first subdivides the ambient space into cubes (in which case the algorithm is called a marching cube algorithm), or simplices [5,6,9] (in which case the algorithm is called a marching tetrahedra algorithm). One starts with a cube or a simplex (cell) in which a part of the zero set of the function is contained, and finds a piecewise-linear approximation of the zero set in that cell. One then propagates or marches to adjacent cells that also contain the zero set and approximates the zero set in that cell. This process can be continued until all cells that intersect the zero set have been visited.
For the marching cube algorithm [4], one also has to decide how to approximate the zero set inside a cube. As observed by Dürst, there is in general no canonical way how to do this due to ambiguous configurations, see [10] for an extensive discussion in the three-dimensional setting. For the marching tetrahedra algorithm, there is a canonical way to construct a piecewise-linear approximation of the zero set, as we will discuss below. However, the result of the algorithm is still not necessarily topologically correct.

Guarantees for Isosurfacing
For the marching simplex algorithm [5] in arbitrary dimensions, bounds have been given on the one-sided Hausdorff distance between the zero set of f and its PL approximation, and also on the difference between the gradient of f and the gradient of the PL approximation. It can be proven that the result of the algorithm is a manifold under appropriate assumptions [11,12].
An important requirement in the work of Allgower and Georg [12] is that the zero set avoids simplices that have dimension less than the codimension, see [12,Definition 12.2.2] and the text above [12,Theorem 15.4.1]. The idea to avoid these low-dimensional simplices originates with Whitney [13], with whom Allgower and George [11,12] were apparently unfamiliar. Very heavy perturbation schemes for the vertices of the ambient triangulation T are needed to ensure that the manifold stays sufficiently far from simplices in the ambient triangulation that have dimension less than the codimension of the manifold [13,14]. Various techniques have been developed to compute such perturbations with guarantees. They typically consist in perturbing the position of the sample points or in assigning weights to the points. Complexity bounds are then obtained using volume arguments. See, for example, [15][16][17][18]. However, these techniques suffer from several drawbacks. The constants in the complexity depend exponentially on the ambient dimension. Moreover, the analysis assumes that the probability of the simplices of dimension less than the codimension to intersect the manifold is zero, which is not true when dealing with finite precision. As a result, the actual implementations we are aware of fail to work well in practice except in very simple cases.
More complete correctness results have been achieved in three dimensions in the computational geometry community: Boissonnat, Cohen-Steiner and Vegter [19] base their proof on a combination of Morse theory and simplicial collapses. Vegter and Plantinga's proof [20] is in its philosophy closely related to normal surface theory, see, for example, [21], but relies rather heavily on case analysis. The results of [19,20] seem not extendable to higher dimensions.

Triangulating General Manifolds (Without Boundary)
The approximation of a manifold that is the zero set of a function is an example of the more general question of how to triangulate a manifold. It is known that C 1 manifolds are triangulable, see, for example, [13], and algorithms have been proposed recently to triangulate smooth manifolds [14,17,22,23]. However, all known methods use intricate perturbation schemes to guarantee the correctness of the triangulation algorithms when the intrinsic dimension of the manifold exceeds 2. As for the case of isomanifolds, perturbation schemes work fine in theory but the constants are miserable and the methods do not work in practice in high dimensions.

Manifolds with Boundary
In this paper, we also consider the piecewise-linear approximation of manifolds with boundary (that are given as a zero set) and briefly mention the extension to stratifolds. Apart from some Delaunay-based work on triangulations of stratifolds in three dimensions [24][25][26][27][28], we are not aware of similar results on manifolds with boundary. Significant effort also went in the detection of strata, in this case in arbitrary dimension, see, for example, [29][30][31].

Contribution
This paper contains three main results, the first two (Theorem 25 and Corollary 27) concern manifolds without boundary and the third (Theorem 48) manifolds with boundary.
We state here simplified versions of the statements to be fully described later on.

Isomanifolds (without boundary)
Let f : R d → R d−n be a smooth function and suppose that 0 is a regular value of f , meaning that at every point x such that f (x) = 0, the Jacobian of f is nondegenerate. Assume that T is a triangulation of R d . Define the function f P L as the linear interpolation of the values of f at the vertices if restricted to a single simplex σ ∈ T . Then, Theorem 25 (Ambient isotopy) The zero set of f P L is a manifold that is ambient isotopic to the zero set of f , provided that the triangulation T is sufficiently fine and thick.
We recall that the thickness of a simplex is the ratio of the height (smallest altitude) over the longest edge length, and it is a measure for the quality or how well shaped a simplex is.
Corollary 27 (Bound on the Fréchet distance) The Fréchet distance between f P L and f is of the order of D 2 , where D is the longest edge length of T .
We also give a variant of a result due to Allgower and George [11]: Proposition 10 The difference between the gradient of f and the gradient of its piecewise-linear approximation is of order d D inside each simplex of T . Isomanifolds with boundary Suppose that apart from f we are also given another function f ∂ : R d → R and f ∂,P L is defined similar to f P L . Write f i for the ith component of f . Let us further assume that the zero set is regular in the following sense: The gradients of f i span a (d − n)-dimensional space at each point of f −1 (0) and the gradients of , provided that the triangulation T is sufficiently fine and thick.
An important aspect of these results is that they hold under mild conditions: they simply ask for a sufficiently fine and thick triangulation T . In contrast to previous results on the triangulation of manifolds, no perturbations are needed to guarantee topological correctness.
Our method provides guarantees on the piecewise-linear (PL) approximation of isomanifolds, regarding the topology, the Fréchet distance and the approximation of the gradients (the latter was already known to Allgower and Georg [11]).
However, we stress that it does not give lower bounds on the quality of the linear pieces in the PL approximation. This is a clear difference with previous methods [13,14,16,23] whose output is a thick triangulation. Although this is an appealing property, it complicates the analysis further and requires unpractical perturbation schemes. Such perturbation techniques could be added to our method to improve the simplex quality (to some limited extent). However, they are not required to make the algorithm work and to obtain the guarantees mentioned above.
The techniques used in this paper are also different from many of the standard tools and do not rely on Delaunay triangulations [32,33], the closed ball property [15,34,35], Whitney's lemma [36] or collapses [37]. The current paper mainly relies on the non-smooth implicit function theorem [38] with some Morse theory.

Outline
The rest of this paper is subdivided as follows. In Sect. 2, we treat closed isomanifolds, i.e. compact manifolds without boundary. In Sect. 3, we treat isomanifolds with boundary. Extension to general isostratifolds is briefly discussed in Sect. 4. In the final section, we quantify the robustness of the method by studying how much the zero set of f changes if f is perturbed slightly in the C 1 -sense [39].
This paper is closely related to another paper where the data structure needed to efficiently propagate along the manifold is presented [40]. Altogether, these two papers show that one can construct PL approximations of isomanifolds in space and time polynomial in the resolution 1/D of the ambient triangulation and in the dimension d of the ambient space.

Isomanifolds (Without Boundary)
Let f : R d → R d−n be a smooth (C 2 suffices) function and suppose that 0 is a regular value of f , meaning that at every point x such that f (x) = 0, the Jacobian of f is non-degenerate. Then, the zero set of f is an n-dimensional manifold as a direct consequence of the implicit function theorem, see, for example, [41,Section 3.5]. We further assume that f −1 (0) is compact. As in [11] we consider a triangulation T of R d . The function f P L is the linear interpolation of the values of f at the vertices if restricted to a single simplex σ ∈ T , i.e.
where the λ v are the barycentric coordinates of x with respect to the vertices of σ . For any function g : We prove that under certain conditions there is an ambient isotopy from the zero set of f to the zero set of f P L . The proof will be using the piecewise-smooth map which interpolates between f and f P L and is based on the generalized implicit function theorem. We are, by definition, only interested in f −1 (0) and so can ignore points that are sufficiently far from this zero set. More precisely, we observe the following: if f i (x) is positive for all x in a geometric simplex σ , then so is f i P L (x) because f i P L (x) is a convex combination of the (positive) values at the vertices. This in turn implies that 1] as, for each τ , it is a convex combination of positive numbers. The same argument holds for negative values. So we see that The results will be expressed using constants defined in terms of f and the ambient triangulation T .

Definition 2 We define
D : the longest edge length of a simplex in T 0 (5) T : the smallest thickness of a simplex in T 0 , where -Gram(∇ f ) denotes the Gram matrix whose elements are ∇ f i ·∇ f j where · stands for the dot product, -λ min (x) denotes the smallest absolute value of the eigenvalues of Gram(∇ f (x)), 1 -Hes( f ) = (∂ k ∂ l f i ) k,l denotes the Hessian matrix of second order derivatives, -| · | denotes the Euclidean norm of a vector and · 2 the operator 2-norm of a matrix. 2 -The thickness is the ratio of the height (smallest altitude) over the longest edge length.
We will assume that γ max , λ min , α max , D, T ∈ (0, ∞). The constant λ min quantifies how close 0 is to not being a regular value of f . D is a measure of the size of the simplices of T . We will call δ = 1/D the resolution of T . The thickness is a quality measure of a simplex. A good choice for T is the Coxeter triangulation of type A d , see [42,43], or the related Freudenthal triangulations, see [3,[44][45][46], which can be defined for different values of D while keeping T constant (for a given dimension d).
Our results hold for any dimensions d and n. We are especially interested in the case where the ambient dimension d is large. We thus consider d and D as the two main parameters. For our bounds, we will give both exact and asymptotic expressions. The asymptotic expressions are given to emphasize the dependency on the two most important parameters d and D, and hold for any d, D and T such that d D < T and for any fixed positive γ max , λ min , α max . For Coxeter triangulations of typeÃ d (the only type of Coxeter triangulations mentioned in this paper), we have T > 2 √ d+2 , 1 Because a Gram matrix is a symmetric square matrix, its eigenvalues are well defined and real. 2 The operator norm is defined as A p = max x∈R n |Ax| p |x| p , with | · | p the p-norm on R n .  [43]. Therefore, the condition d D < T is satisfied for these triangulations when D < (2d √ d + 2) −1 . For convenience, exact expressions are gathered in Appendix A.

The Result
We are going to construct an ambient isotopy based on (1), see Fig. 1 for a pictorial overview. In fact, the map τ → {x | F P L (x, τ ) = 0} gives an ambient isotopy between the zero set of F P L (x, 0), which is identical to the smooth isosurface f −1 (0), and the zero set of F P L (x, 1), which is the PL approximation f −1 P L (0). The latter can be turned into a triangulation of the isosurface f −1 (0) by triangulating the non-simplicial cells using barycentric subdivision. We will also bound the Fréchet distance between f −1 (0) and f −1 P L (0) . Proving the ambient isotopy consists of three technical steps. The first two consume most of the space in the proof namely: is a smooth manifold, under certain conditions (Corollary 13). -Global step. We prove that F −1 P L (0) is a manifold, under certain conditions, using techniques from non-smooth analysis (Corollary 24).
A crucial ingredient will be the implicit function theorem and its non-smooth extension. Along the way, we shall also see that F −1 P L (0) is never tangent to the τ = c planes, where c is a constant. The gradient of (x, τ ), → τ in R d × R is (0, 1). Projecting this vector onto the tangent space of F −1 P L (0) gives the gradient of (x, τ ), → τ restricted to F −1 P L (0). Because of the non-tangency, this projection is nonzero. So the gradient field of the function (x, τ ), → τ restricted to F −1 P L (0), is piecewise-smooth (because F −1 P L (0) is piecewise-smooth) and never vanishes.
The third step is similar to a standard observation in Morse theory [47,48], with the exception that we now consider piecewise-smooth instead of smooth vector fields. We refer to Milnor [47] for an excellent introduction and to Lemma 2.4 and Theorem 3.1 in particular. Proof This is a straightforward consequence of the existence and uniqueness of the solution to a differential equation.
Bounds on the gradient of τ on the manifold give a bound on the Fréchet distance, which is defined in the following.

Preliminaries
The following elementary lemma will be useful.

Lemma 5
Proof The first statement follows from the fact that the supremum of the absolute value of a derivative of a function bounds the Lipschitz constant of the function. The second statement follows from standard bounds on matrix norm (see, for example, [49,Equation (2.3.11)]) together with (4). These bounds imply that for all k, l and i. Arguing as before, we deduce a bound on the Lipschitz constant of ∂ l f i : Bound (8) now follows.

The Implicit Function Theorem
The main technical tool to prove the existence of the ambient isotopy from f −1 (0) to f −1 P L (0) is the implicit function theorem which we recall now. To prove the existence of the isotopy from f −1 (0) to f −1 P L (0), we will apply the implicit function theorem (Theorem 6) to several functions g that are close to f and we will therefore need to prove that their Jacobians are of maximal rank. A matrix has maximal rank if and only if the Gram matrix of its columns has a nonzero determinant or, equivalently, nonzero eigenvalues. In our context, we will need lower bounds on the absolute values of the eigenvalues of the Gram matrices Gram(∇g), given the lower bound λ min on the absolute values of the eigenvalues of Gram(∇ f ).

Eigenvalues and Perturbations
We will follow the convention that the eigenvalues of the matrices we consider are sorted by increasing order of their absolute values, i.e. |λ i | ≤ |λ j | if i ≤ j. We first recall Weyl's perturbation theorem that bounds the difference between the ith eigenvalues of two symmetric matrices:

Lemma 7 (Weyl's bound, Corollary III.2.6 of [50]) Let A andÃ = A + E be two symmetric (or Hermitian) matrices and write λ i andλ i for the eigenvalues of A and A, respectively. Then,
where · p denotes the p-norm.
We further note that E 2 ≤ E F where · F denotes the Frobenius norm, see [49, (2.3.7)]. By definition of the Frobenius norm, we have that |E i j | ≤ e max , for all

Estimates for a Single Simplex
In this section, we concentrate on a single simplex σ and write f L for the linear function whose values on the vertices of σ coincide with f . In other words, f L is the linear extension of the interpolation of f . Note that f L coincides with f P L within the geometric simplex σ (but not necessarily outside).

Estimates on the Linear Approximation f L and Its Gradient
We need a simple estimate similar to Proposition 2.1 of Allgower and George [11].

Lemma 9
Let σ ⊂ T 0 and let f L be as described above. Then, for all x ∈ σ , We included a proof for completeness.

Proof
Let v k be a vertex of σ . Taylor's theorem, see, for example, [41,Theorem 2.8.4], yields that with Thanks to the bounds on R(v k ) and Cauchy-Schwarz, one has We will also be using an estimate similar to Proposition 2.2 of Allgower and George [11].

Proposition 10
Let σ ⊂ T 0 and let f L be as described above. Then for all x in the simplex σ .
We provide a proof for completeness. In this proof, we use the following variant of a common refinement of two sets:

Claim 11 Suppose we are given two finite sets
Proof The proof goes via intersection as indicated in Fig. 2. In the figure, the sets A, B, C are represented as a union of intervals of lengths a i , b j , c k , respectively. We observe the following: -In the figure, the set C is obtained by the common refinement of the intervals.
We thus have k max ≤ i max + j max with equality if there are no ι, ι such that -Any interval of A or B is the union of consecutive intervals of C, and the result follows. Fig. 2 The construction of the c k s from Claim 11 Proof of Proposition 10 We use the same notation as before and again use that with Subtracting Because f L is the linear interpolation of f , we have Let now u and v be two points of σ . Writing u and v in terms of their barycentric coordinates with respect to the vertices v k of σ , we have We can now apply Claim 11 to the two sets The refined set is denoted by {π l |l = 0, ..., L}, where L ≤ |K |. The refinement associates with each l a k + ∈ K + and a k − ∈ K − . We will write l + = k + (l) and l − = k − (l) to emphasize this dependence. With this notation, we can write By Claim 11, we have We now see that (11) and (13)).
Because the simplex σ contains a ball of radius the smallest altitude over d centred at its barycentre, that is, T D/d with T the thickness, the vector u − w can be chosen to be any vector of length less than T D/d. In particular, we can choose Plugging this choice into (14) gives We stress that the bound in Proposition 10 depends on the quality of the simplices in the ambient triangulation T but not on the shape of the cells of the PL approximation. This is fortunate since we know ambient triangulations of very good quality (e.g. Coxeter triangulations [43]), while we do not have control on the shapes of the cells of the PL approximation which depend on the way the isomanifold intersects T .

Applying the Implicit Function Theorem
Let σ be a simplex of T and let f L be the linear approximation defined above. We now define a homotopy F L : We intend to show that . This follows from the Implicit function theorem (Theorem 6) provided that, for any point such that F L (x, τ ) = 0 in this neighbourhood, the Jacobian is of maximal rank or, equivalently as recalled above, if and only if the Gram matrix of its columns has nonzero eigenvalues. The following lemma will provide lower bounds on the eigenvalues of this Gram matrix.
We denote by ∇ x F L or simply ∇ F L the gradient of the restriction of F L to the x variable and ∇ x,τ F L the gradient of F L . Note that For convenience, we will write 3 and write λ min and λ min for the smallest absolute values of the eigenvalues of G and G, respectively.
. The precise expression of e L is given in (19) and (18).
Proof Let, in addition to the notations of the lemma, G = Gram(∇ x F L ), λ min be the smallest absolute value of the eigenvalues of G , and write G i, j , G i, j and G i, j for the entries of G, G and G, respectively. Proposition 10 yields that The addition of the τ component gives a small extra contribution.
by Lemma 9 Applying Corollary 8, we obtain From (18) and assuming d D < T , we see that , see [43]. The following corollary follows directly from the previous lemma and the discussion before the lemma.
is a smooth manifold.

Transversality with Regard to the -Direction
We now prove that inside each σ × [0, 1] the gradient of τ on F L = 0 is smooth and does not vanish. We need the following straightforward lemma. We include a proof for completeness. Proof We see that

Lemma 14 Now suppose that
We also need to bound the angle of the vectors ∇ x,τ (F i L ) and the x plane, that is Lemma 15 Let Ξ be as above. We have In particular, the manifold F −1 is never tangent to the τ = c planes, where c is a constant, provided that the following transversality condition holds Proof By (16),  (20) and (22), which both holds for D/

The Non-smooth Implicit Function Theorem
For the global result, we need to recall some definitions and results from non-smooth analysis. We refer to [38] for an extensive introduction.

Applying the Non-smooth Implicit Function Theorem
We recall the definition of F P L : Further recall that the closed star of a vertex v in a simplicial complex is the closure of all simplices in the complex that contain v. We will also be using the following remark often.

Remark 21
Let v be a vertex in T , x 1 , x 2 ∈ star(v), then We now have Proof Because (8) and (23), Proposition 10 and Lemma 9).
We generalize Lemma 12 as follows.
. . , μ m are positive weights such that μ 1 + · · · + μ m = 1, and let Λ min be the smallest modulus of the eigenvalues of G. Then, where e P L = O(d 2 D/T ) (assuming d D < T ) and is precisely defined in (26). When T is a Coxeter triangulation of typeÃ d , e P L = O(d 5/2 D).
Proof Let G = Gram(∇ x,τ F L (x 0 , y 0 )) and let λ min be the smallest modulus of the eigenvalues of G. We claim that the elements of the two matrices G and G are pairwise close. Specifically, using the identity by Cauchy-Schwarz and the triangle inequality where g P L is given in Lemma 22. It remains to bound |∇ x,τ F i P L (x k , τ k )|: where we used Lemma 9 and Proposition 10. We conclude that Applying Corollary 8, we get which completes the proof of the lemma.
From the previous lemmas, we immediately have that, where e P L is defined in Lemma 23, the generalized implicit function theorem, Theorem 20, applies to F P L (x, τ ) = 0. In particular, The second technical step of the proof is now completed. The third step follows from an application of Lemma 3. The fact that F L (x, τ ) = 0 is a piecewise-smooth manifold and transversality, as proven in Lemma 15, gives that the gradient of τ is a piecewise-smooth vector field whose flow we can integrate to give an ambient isotopy from the zero set of f to that of f P L .
We summarize in a theorem:

Fréchet Distance
To bound the Fréchet distance, denoted by d F , between the zero sets of f (x) and f P L , it suffices to bound the angle that the gradient of τ (as restricted to F L (x, τ ) = 0) makes with the τ -direction (in R d+1 ). We write e τ for the unit vector in the τ direction (again in R d+1 ). For this, we will use the angle bound of Lemma 15, together with some estimates that are similar in spirit to those in [51,Lemma C.13].

Lemma 26 For any
where e P L = O(d 2 D/T ) is defined in (26) and θ = O(D 2 ) is defined in (21).
where Λ min is defined in Lemma 23. Proposition 10 and |∇( f i )| ≤ γ max give Lemma 15 states that By definition, Ξ is the space orthogonal to e τ (with e τ aligned with the τ direction), so that Hence, by definition of the cosine and (29), we see Using (28) and Cauchy-Schwarz, we then obtain The result now follows thanks to Lemma 23.
Let F P L 0 be the restriction of F P L to F −1 P L (0, τ ). Moreover, let ∇ τ F P L 0 be the gradient of τ restricted to F −1 P L (0), whenever it exists. We want to bound the angle of ∇ τ F P L 0 and the τ -direction. Because the isotopy is given by the gradient flow and we have a bound on the norm of the gradient, the Fréchet distance is bounded. Specifically, the bound is equal to the norm of the gradient since the time we follow the flow is 1.
There is one subtlety. Because the manifold is only piecewise-smooth, we need to take into account the points where ∇ τ F P L 0 is not uniquely defined. Because, for each simplex σ , F L extends to a neighbourhood of σ × [0, 1], there exists a limit of ∇ τ F P L 0 (x i , τ i ) for any sequence (x i , τ i ) that lies in int(σ ) × [0, 1], where int denotes the interior. This means that, if we bound ∇ τ F P L 0 for each simplex, we also bound its limits, where the limits are as just described.
Corollary 27 (Bound on the Fréchet distance) Suppose that the conditions of Theorem 25 are satisfied. Then, Proof Let, as before, Ξ = R d ⊂ R d+1 be the space spanned by the d basis vectors corresponding to the x-directions. Lemma 26 gives, for w ∈ span i (∇ x,τ (F i )), Since the tangent space to F L = 0 is normal to span i (∇ x,τ (F i L )), the same bound holds for sin ∠(∇ τ F P L 0 , e τ ). This means that, as τ ∈ [0, 1], the distance between the begin and the end points of the gradient flow, and thus the Fréchet distance, is bounded by tan ∠(∇ τ F P L 0 , e τ ), that is, The asymptotic dependence follows because sin(arctan(x)) = x √ 1+x 2 , and tan(arcsin(x)) = x √ 1−x 2 .

Isomanifolds with Boundary
We will now consider isomanifolds with boundary. By this, we mean that on top of the function f : R d → R d−n , we will have another function f ∂ : R d → R and the set we consider is ). This is a manifold with boundary if the gradients of f i span a (d − n)-dimensional space at each point of f −1 (0) and the gradients of f i and f ∂ span a (d − n + 1)-dimensional space at each point of , as a consequence of the submersion theorem. We will again write f P L for the PL interpolation of f . Similarly, we write f ∂,P L for the PL interpolation of f ∂ .
We prove that, under certain conditions, there is an isotopy from The conditions are very similar to the conditions we have before but, of course, we need to include bounds on the gradient of f ∂,P L .

Overview of the Proof
We will again construct an isotopy, but in this case it will consist of two steps, see Fig. 3 for a pictorial overview.
-In the first step, we isotope the part of f −1 (0) that is far from f −1 ∂ (0) to its piecewise-linear approximation, while leaving the part of f −1 (0) that is close to f −1 ∂ (0) smooth. We will denote the result by M 1 = (F P L,1 (·, 1)) −1 (0). -In the second step, we consider a (small) tubular neighbourhood around f −1 ∂ (0) as restricted to M 1 by looking at all f −1 ∂ ( ) for | | sufficiently small. 4 We then isotope M 1 ∩ f −1 ∂ ( ) to its piecewise-linear approximation. Again, the isotopy is chosen in such a way that, for relatively large, it leaves M 1 ∩ f −1 ∂ ( ) invariant (for the points such that M 1 is already piecewise-linear). This gives an isotopy of a tubular neighbourhood of M 1 ∩ f −1 ∂ (0) to its piecewiselinear approximation.
We will first partition the manifold in two parts using a smooth bump function φ : R → [0, 1], defined so that φ(y) = 0 in a neighbourhood of zero and φ(y) = 1 if |y| > y 0 , for some y 0 > 0. Such bump functions can be easily constructed, see, for example, [39, Section 2.2]. We will be using the function φ i ( f i ) 2 + f 2 ∂ often. In fact, because it is used so often, it will be convenient to introduce the following shorthand The first step will be using the zero set of the following function: on which we will apply the same gradient flow argument as before.
The resulting set M 1 is the same zero set of f P L as before if we stay sufficiently far away from ∂ M and the isotopy leaves the manifold invariant close to ∂ M. In particular, ∂ M 1 = ∂ M (Fig. 3).
In the second step, we define an isotopy that will act only on a small neighbourhood of ∂ M. Consider the sets B 1 ( ) = M 1 ∩ f −1 ∂ ( ) and, for each of them, define the function where ψ : R → [0, 1] is now a smooth bump function that is 1 in a sufficiently large neighbourhood of zero (somewhat larger than y 0 ) and zero outside some compact set. Using the result for isomanifolds (with some modifications), we can prove that each individual set B 1 ( ) is isotopic to f −1 P L (0) ∩ f −1 ∂,P L ( ) for small while, for sufficiently large , it leaves the set invariant.

Step 1
The proof closely follows the proof for the case without boundary in Sect. 2. The main technical difficulty will be to provide bounds that serve as the counterparts to Lemma 22 for both steps in the proof. To be able to do so, we first need to discuss bounds on the bump functions φ and ψ.

Inside a Single Simplex
Similar to Corollary 13, we need a condition that ensures that the zero set of F i P L,1 (x, τ ) restricted to σ × [0, 1]is a smooth manifold. In fact, similar to (15), we define where φ is as defined above. Observe that F i L,1 (x, τ ) can be extended to a neighbourhood of σ × [0, 1].

Remark 29
For the constants, it is better if y 0 can be chosen as large as possible, but we need y 1 to be quite a bit larger than y 0 . In turn, we cannot choose y 1 arbitrarily large because this would mean that the gradient field ∇ f ∂ | f −1 (0) (seen as restricted on f −1 (0)) would never vanish. The latter is in general impossible thanks to the hairy ball theorem [52].
We introduce the following definition that complements Definition 2:

Definition 30
where G B 1 = Gram(∇ x,τ F L,1 ) and λ min (A) denotes, as before, the smallest absolute value of the eigenvalues of matrix A.
We have then the analogue of Lemma 12:

Lemma 31
We have, (41) and (38). For Coxeter triangula- Proof We start with an estimate on the individual ∇ x,τ F i L,1 (x, τ ). As noted before, We now write G i, j and G B 1 i, j for the elements (i, j) of G and G B 1 , respectively. Expanding yields We now see by Cauchy-Schwarz, the triangle inequality and Eq. (38) that and e B 1 t is defined in (38). The following corollary is then the analogue of Corollary 13: holds, then F −1 L,1 (0) is a smooth manifold inside an neighbourhood of σ × [0, 1].

Transversality with Regard to the -Direction
We note that, similar to Lemma 15, we have Lemma 33 Using the notation of Lemma 15, we have

Proof
The proof is identical to the proof of Lemma 15 with the replacement of 4d Dα max T in the denominator by e B 1 t The latter constant is a consequence of (38).

Now, similar to Corollary 16, we find that
Corollary 34 (Transversality with respect to τ for Step 1) Assume that both the regularity condition (42) and the transversality condition

Global Result
We now have to prove that F −1 P L,1 (0) is a manifold. For this, we again employ the generalized implicit function theorem. But first of all, we need the following bound, which is similar to Lemma 22. Note that ∇ x,τ F i (x 0 , τ 0 ) is well defined as soon as x 0 lies in the interior of a d-simplex in T .

Lemma 35
Assuming that the gradients are well defined, we have (44). When T is a Coxeter triangulation of typeÃ d , g B 1 P L = O(d 3/2 D). Proof By expansion, we see that (23), (8), and (38) twice) This completes the proof.
Suppose x 0 , x 1 , . . . , x m ∈ star(v), τ 0 , . . . , τ m ∈ [0, 1] and that ∇ x,τ F i P L,1 (x i , τ i ) is well defined for all i. Further, assume that μ 1 , . . . , μ m are positive weights such that μ 1 + · · · + μ m = 1. We write Proof The proof is more or less the same as the proof of Lemma 23, but with more complicated bounds. We assume that x 0 ∈ star(v) and τ 0 ∈ [0, 1] are such that ∇ x,τ F i (x 0 , τ 0 ) is well defined (i.e. x 0 lies in the interior of a d-simplex of T ). Lemma 31 gives that Using ∇( f i ) ≤ γ max and (38), we note that We want to use Weyl's bound to determine a bound on the smallest absolute value of the eigenvalues of G B 1 . Writing G B 1 i, j and G B 1 i, j for element (i, j) of matrices G B 1 (x 0 , τ 0 ) and G B 1 , respectively, we show that G B 1 i, j and G B 1 i, j (x 0 , τ 0 ) are pairwise close (compare to (40)).
by Lemma 35 and (45). (46) Using the result of Lemma 31 and invoking Corollary 8 once more gives Lemma 36 immediately yields that holds, then the generalized implicit function theorem, Theorem 20, applies to F P L,1 (x, τ ) = 0. In particular, F −1 P L,1 (0) is a manifold.
We stress again that, inside the set , the zero set of F P L,1 (x, 1) coincides with the zero set of f P L (x).

Step 2
Before we can proceed, we have to specify the bump function ψ. We suppose that (the constants have not been optimized) In particular,

Remark 39
We stress that the choice of y 1 and y 2 for the function ψ is different from the choice we made for φ, in the first step of the proof.
First, we stress that the zero set of F P L,2, (x, 1) coincides with the zero set of Secondly, we now claim the following:

Lemma 40
The zero set of F P L,2, (x, 1) is a subset of the zero set of f P L (x), for each .
Proof We focus on the first d − n coordinates of (34). We see that where we used that We can further rewrite (49), where we used the same argument as before.
The technical result that remains to be proven is the counterpart of Theorem 25 for F P L,2, (x, τ ) and for each sufficiently small . To be precise, it suffices for ≤ 2y 0 . We remark that it is likely that this bound on can be improved.
We again follow the same path to prove this result. That is, we first concentrate on a single simplex and prove that inside that simplex the zero set of F P L,2, is a smooth manifold on which the gradient of τ as restricted to the manifold does not vanish. We then prove that the zero set of F P L,2, is globally a manifold.

Assumptions and Notations
Because we are now faced with both f (x) and f ∂ (x), we need to introduce a bound on how far the gradients of all the entrees of these functions are from being collinear. We write Before we were only interested in the set T 0 . Similarly here, we sometimes concentrate on a neighbourhood of the zero set of both f ∂ and f . Therefore, we write T B for all σ ∈ T such that ( l ( f l ) 2 We also write G B = Gram(∇ f B ) and λ B min for the minimal absolute value of the eigenvalues of G B , where the minimization is over all simplices in the set T B ∩ T 0 . The restriction to the set T B ∩ T 0 is important, because if the minimization would be just over T 0 , G B would generically be 0 as a consequence of the hairy ball theorem.
We note that by taking gradients the constant drops from the expression, so that the properties we now define are independent of . For the lengths of the gradients of f B , we define, for all 1 ≤ i ≤ d − n + 1. Similar to α max , we define α B max as the bound on the operator 2-norm of all Hessians of f B , that is, We stress that that α max ≤ α B max . We use the same notation for the ambient triangulation T , the lower bound on the thickness of the simplices T and upper bound on the longest edge length D. We also need to introduce a bound on the differential of the bump function ψ. Similar to (35) because we picked y 1 = 101 100 y 0 and y 2 = 2y 0 , for ψ, see Definition 38 and Remark 39.

Inside a Single Simplex
Similar to Lemma 31, we now give a condition that ensures that the zero set of F P L,2, (x, τ ) is smooth inside σ × [0, 1]. In fact, similar to (15), we define which can be extended to a neighbourhood of σ × [0, 1]. We also write G B 2 = Gram(∇ x,τ F P L,2, andλ B 2 min for the smallest absolute value of the eigenvalues of G B 2 . Proof We start with an estimate on the individual ∇ x,τ F i L,2, (x, τ ). We will write (v, w) i for the ith coordinate of the composed vector (v, w). We now see that
is a smooth manifold inside an small neighbourhood of σ × [0, 1] provided that the following regularity condition holds where e B 2 L is defined in Lemma 41.

Transversality with Regard to the -Direction
Once more, similar to Lemma 15 we have

Lemma 43
Let Ξ be as in Lemma 15. We have

Global Result
We now have to prove that F −1 P L,2, (0) is a manifold, for all sufficiently small . For this we first need the following bound, which is similar to the one in Lemma 35.

Lemma 45
Let v be a vertex in T , x 1 , x 2 ∈ star(v), and τ 1 , τ 2 ∈ [0, 1], such that ∇ x,τ F i P L,2, (x 1 , τ 1 ) and ∇ x,τ F i P L,2, (x 2 , τ 2 ) are well defined, then Proof The proof follows the same steps as the proof of Lemma 35. By expansion, we see that Suppose x 0 , x 1 , . . . , x m ∈ star(v), τ 0 , . . . , τ m ∈ [0, 1], and that, for all i, ∇ x,τ F i P L,2 (x i , τ i ) is well defined. Further assume that μ 1 , . . . , μ m are positive weights such that μ 1 + · · · + μ m = 1. We write Proof The proof is more or less the same as the proof of Lemma 36. Let x 0 ∈ star(v) and τ 0 ∈ [0, 1], be such that ∇ x,τ F i P L,2, (x 0 , τ 0 ) is well defined. Note that it is sufficient for x 0 to lie in the interior of a d-simplex in T . Lemma 41 gives that Using ∇( f i B ) ≤ γ B max and (54), we get We want to bound the smallest absolute value of the eigenvalues of G B 2 = Gram(∇ x,τ F L,2, )). We proceed similarly to (46) (with F P L,1 replaced by F i P L,2, ). Let G B 2 be as in (60), and denote by G B 2 i, j and G B 2 i, j the (i, j) elements of G B 2 and G B 2 , respectively, and by Λ B 2 min the smallest absolute value of the eigenvalues of G B 2 . (62)).
Thanks to Corollary 8 and Lemma 41, we have that Lemma 46 immediately yields that

Corollary 47 (The generalized implicit function theorem in
Step 2) If the regularity condition the generalized implicit function theorem, Theorem 20, applies to F P L,1 (x, τ ) = 0. In particular, F −1 P L,1 (0) is a manifold.

Fréchet Distance
The bounds on the Fréchet distance can be achieved in the same way as before.

Theorem 49
Suppose that the conditions of Theorem 48 are satisfied. Then, Proof We apply the same argument as in Lemma 26 and Corollary 27, for both steps of the proof. This yields the sum of two terms that are of the same form as (31). For the first step, we need the following substitutions: , as a consequence of Lemma 33.
λ min − e P L is replaced by λ min − e B 1 P L , as a consequence of Lemma 36. -γ max + 4d Dα max T is replaced by γ max + 2D 2 α max + e B 1 t , as a consequence of (38).
For the second step, we need the following substitutions: , as a consequence of Lemma 43.
-√ λ min − e P L is replaced by λ min − e B 2 P L , as a consequence of Lemma 46. -γ max + 4d Dα max T is replaced by γ B max + e B 2 t as a consequence of (54).
This yields

Isostratifolds
There is no obstruction that prevents us from extending the approach above to isostratifolds. By isostratifolds, we mean stratifolds that are given by the zero sets of functions and inequalities. For example, suppose that we want to find a PL approximation of the unit sphere centred at 0 in R 3 including the PL approximations of the intersections of the sphere with slightly deformed x = 0, y = 0, and z = 0-planes, as depicted in Fig. 4. This would also give PL approximations of the respective 'octants' of the sphere.
We could follow the same procedure as for a manifold with boundary to give precise bounds on the longest edge length D of the ambient triangulation that ensure that the PL approximation is correct. However, this would mean that we have to introduce an extra bump function for each stratum as well as an extra isotopy. Even though this should be relatively straightforward, finding the precise constants involved would become prohibitively lengthy.

Robustness
Suppose that f and f δ are smooth functions and moreover f δ is small in terms of the C 1 -topology. Thanks to the implicit function theorem, we know that if 0 is a regular value of f , the zero set of f and the zero set of the slightly perturbed function f + f δ are isotopic. We now give quantitative conditions that guarantee that f −1 (0) and ( f + f δ ) −1 (0) are ambient isotopic. Let f , f δ : R d → R d−n , and writẽ α max = max where τ ∈ [0, 1]. We, again, first establish that the zero set of this function is a (n + 1)-dimensional manifold. Secondly, we will see that the gradient of τ restricted to this (n + 1)-dimensional manifold never vanishes. As we have seen in the previous sections, this suffices to establish the isotopy F −1 (0) from f −1 (0) to ( f + f δ ) −1 (0), by Lemma 3.
We find that Acknowledgements First and foremost, we acknowledge Siargey Kachanovich for discussions. We thank Herbert Edelsbrunner and all members of his group, all former and current members of the Datashape team (formerly known as Geometrica), and André Lieutier for encouragement. We further thank the reviewers of Foundations of Computational Mathematics and the reviewers and program committee of the Symposium on Computational Geometry for their feedback, which improved the exposition.
Funding Open access funding provided by Institute of Science and Technology (IST Austria).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Notations and Overview of Constants
For notations, we followed the following rules: 1. Greek letters, except τ and , are for constants related to functions.

Bounds on Eigenvalues
Bounds on Angles