Anisotropic triangulations via discrete Riemannian Voronoi diagrams

The construction of anisotropic triangulations is desirable for various applications, such as the numerical solving of partial differential equations and the representation of surfaces in graphics. To solve this notoriously difficult problem in a practical way, we introduce the discrete Riemannian Voronoi diagram, a discrete structure that approximates the Riemannian Voronoi diagram. This structure has been implemented and was shown to lead to good triangulations in $\mathbb{R}^2$ and on surfaces embedded in $\mathbb{R}^3$ as detailed in our experimental companion paper. In this paper, we study theoretical aspects of our structure. Given a finite set of points $\cal P$ in a domain $\Omega$ equipped with a Riemannian metric, we compare the discrete Riemannian Voronoi diagram of $\cal P$ to its Riemannian Voronoi diagram. Both diagrams have dual structures called the discrete Riemannian Delaunay and the Riemannian Delaunay complex. We provide conditions that guarantee that these dual structures are identical. It then follows from previous results that the discrete Riemannian Delaunay complex can be embedded in $\Omega$ under sufficient conditions, leading to an anisotropic triangulation with curved simplices. Furthermore, we show that, under similar conditions, the simplices of this triangulation can be straightened.


Introduction
Anisotropic triangulations are triangulations whose elements are elongated along prescribed directions. Anisotropic triangulations are known to be well suited when solving PDE's [12,22,27]. They can also significantly enhance the accuracy of a surface representation if the anisotropy of the triangulation conforms to the curvature of the surface [18].
Many methods to generate anisotropic triangulations are based on the notion of Riemannian metric and create triangulations whose elements adapt locally to the size and anisotropy prescribed by the local geometry. The numerous theoretical and practical results [1] of the Euclidean Voronoi diagram and its dual structure, the Delaunay triangulation, have pushed authors to try and extend these well-established concepts to the anisotropic setting. Labelle and Shewchuk [20] and Du and Wang [14] independently introduced two anisotropic Voronoi diagrams whose anisotropic distances are based on a discrete approximation of the Riemannian metric field. Contrary to their Euclidean counterpart, the fact that the dual of these anisotropic Voronoi diagrams is an embedded triangulation is not immediate, and, despite their strong theoretical foundations, the anisotropic Voronoi diagrams of Labelle and Shewchuk and Du and Wang have only been proven to yield, under certain conditions, a good triangulation in a two-dimensional setting [8,9,11,14,20].
Both these anisotropic Voronoi diagrams can be considered as an approximation of the exact Riemannian Voronoi diagram, whose cells are defined as V g (p i ) = {x ∈ Ω | d g (p i , x) ≤ d g (p j , x), ∀p j ∈ P\p i }, where d g (p, q) denotes the geodesic distance. Their main advantage is to ease the computation of the anisotropic diagrams. However, their theoretical and practical results are rather limited. The exact Riemannian Voronoi diagram comes with the benefit of providing a more favorable theoretical framework and recent works have provided sufficient conditions for a point set to be an embedded Riemannian Delaunay complex [2,16,21]. We approach the Riemannian Voronoi diagram and its dual Riemannian Delaunay complex with a focus on both practicality and theoretical robustness. We introduce the discrete Riemannian Voronoi diagram, a discrete approximation of the (exact) Riemannian Voronoi diagram. Experimental results, presented in our companion paper [26], have shown that this approach leads to good anisotropic triangulations for two-dimensional domains and surfaces, see Figure 1.
We introduce in this paper the theoretical side of this work, showing that our approach is theoretically sound in all dimensions. We prove that, under sufficient conditions, the discrete Riemannian Voronoi diagram has the same combinatorial structure as the (exact) Riemannian Voronoi diagram and that the dual discrete Riemannian Delaunay complex can be embedded as a triangulation of the point set, with either curved or straight simplices. Discrete Voronoi diagrams have been independently studied, although in a two-dimensional isotropic setting by Cao et al. [10].

Riemannian geometry
In the main part of the text we consider an (open) domain Ω in R n endowed with a Riemannian metric g, which we shall discuss below. We assume that the metric g is Lipschitz continuous. The structures of interest will be built from a finite set of points P, which we call sites.

Riemannian metric
A Riemannian metric field g, defined over Ω, associates a metric g(p) = G p to any point p of the domain. This means that for any v, w ∈ R n we associate an inner product v, w g = v t g(p)w, in a way that smoothly depends on p. Using a Riemannian metric, we can associate lengths to Figure 1: Left, the discrete Riemannian Voronoi diagram (colored cells with bisectors in white) and its dual complex (in black) realized with straight simplices of a two-dimensional domain endowed with a hyperbolic shock-based metric field. Right, the discrete Riemannian Voronoi diagram and the dual complex realized with curved simplices of the "chair" surface endowed with a curvature-based metric field [26].
curves and define the geodesic distance d g as the minimizer of the lengths of all curves between two points. When the map g : p → G is constant, the metric field is said to be uniform. In this case, the distance between two points x and y in Ω is d G Most traditional geometrical objects can be generalized using the geodesic distance. For example, the geodesic (closed) ball centered on p ∈ Ω and of radius r is given by B g (p, r) = {x ∈ Ω | d g (p, x) ≤ r}. In the following, we assume that Ω ⊂ R n is endowed with a Lipschitz continuous metric field g.
We define the metric distortion between two distance functions d g (x, y) and d g (x, y) to be the function ψ(g, g ) such that for all x, y in a small-enough neighborhood we have: 1/ψ(g, g ) d g (x, y) ≤ d g (x, y) ≤ ψ(g, g ) d g (x, y). Observe that ψ(g, g ) ≥ 1 and ψ(g, g ) = 1 when g = g . Our definition generalizes the concept of distortion between two metrics g(p) and g(q), as defined by Labelle and Shewchuk [20] (see Appendix B).

Geodesy
Let v ∈ R n . From the unique geodesic γ satisfying γ(0) = p with initial tangent vectorγ = v, one defines the exponential map through exp(v) = γ (1). The injectivity radius at a point p of Ω is the largest radius for which the exponential map at p restricted to a ball of that radius is a diffeomorphism. The injectivity radius ι Ω of Ω is defined as the infimum of the injectivity radii at all points. For any p ∈ Ω and for a two-dimensional linear subspace H of the tangent space at p, we define the sectional curvature K at p for H as the Gaussian curvature at p of the surface exp p (H).
In the theoretical studies of our algorithm, we will assume that the injectivity radius of Ω is strictly positive and its sectional curvatures are bounded.

Power protected nets
Controlling the quality of the Delaunay and Voronoi structures will be essential in our proofs. For this purpose, we use the notions of net and of power protection.

Anisotropic triangulations via discrete Riemannian Voronoi diagrams 5
Power protection of point sets Power protection of simplices is a concept formally introduced by Boissonnat, Dyer and Ghosh [2]. Let σ be a simplex whose vertices belong to P, and let B g (σ) = B g (c, r) denote a circumscribing ball of σ where r = d g (c, p) for any vertex p of σ. We call c the circumcenter of σ and r its circumradius.
For 0 ≤ δ ≤ r, we associate to B g (σ) the dilated ball B +δ g (σ) = B(c, √ r 2 + δ 2 ). We say that σ is δ-power protected if B +δ g (σ) does not contain any point of P \ Vert(σ) where Vert(σ) denotes the vertex set of σ. The ball B +δ g is the power protected ball of σ. Finally, a point set P is δ-power protected if the Delaunay ball of its simplices are δ-power protected.
Nets To ensure that the simplices of the structures that we shall consider are well shaped, we will need to control the density and the sparsity of the point set. The concept of net conveys these requirements through sampling and separation parameters.
The sampling parameter is used to control the density of a point set: if Ω is a bounded domain, P is said to be an ε-sample set for Ω with respect to a metric field g if d g (x, P) < ε, for all x ∈ Ω. The sparsity of a point set is controlled by the separation parameter: the set P is said to be µ-separated with respect to a metric field g if d g (p, q) ≥ µ for all p, q ∈ P. If P is an ε-sample that is µ-separated, we say that P is an (ε, µ)-net.

Riemannian Delaunay triangulations
Given a metric field g, the Riemannian Voronoi diagram of a point set P, denoted by Vor g (P), is the Voronoi diagram built using the geodesic distance d g . Formally, it is a partition of the domain in Riemannian Voronoi cells The Riemannian Delaunay complex of P is an abstract simplicial complex, defined as the nerve of the Riemannian Voronoi diagram, that is the set of simplices Del g (P) = {σ | Vert(σ) ∈ P, ∩ p∈σ V g (p) = 0}. There is a straightforward duality between the diagram and the complex, and between their respective elements.
In this paper, we will consider both abstract simplices and complexes, as well as their geometric realization in R n with vertex set P. We now introduce two realizations of a simplex that will be useful, one curved and the other one straight.
The straight realization of a n-simplex σ with vertices in P is the convex hull of its vertices. We denote it by σ. In other words, The curved realization, notedσ is based on the notion of Riemannian center of mass [19,15]. Let y be a point ofσ with barycentric coordinate λ p (y), p ∈ σ. We can associate the energy functional E y (x) = 1 2 p∈σ λ p (y)d g (x, p) 2 . We then define the curved realization of σ as The edges ofσ are geodesic arcs between the vertices. Such a curved realization is well defined provided that the vertices of σ lie in a sufficiently small ball according to the following theorem of Karcher [19].
Theorem 3.1 (Karcher). Let the sectional curvatures K of Ω be bounded, that is Λ − ≤ K ≤ Λ + . Let us consider the function E y on B ρ , a geodesic ball of radius ρ that contains the set {p i }.

The discrete Riemannian Voronoi Diagram
To define the discrete Riemannian Voronoi diagram of P, we need to give a unique color to each site of P and to color the vertices of the canvas accordingly. Specifically, each canvas vertex is colored with the color of its closest site. Definition 4.1 (Discrete Riemannian Voronoi diagram). Given a metric field g, we associate to each site p i its discrete cell V d g (p i ) defined as the union of all canvas simplices with at least one vertex of the color of p i . We call the set of these cells the discrete Riemannian Voronoi diagram of P, and denote it by Vor d g (P).
Observe that contrary to typical Voronoi diagrams, our discrete Riemannian Voronoi diagram is not a partition of the canvas. Indeed, there is a one canvas simplex-thick overlapping since each canvas simplex σ C belongs to all the Voronoi cells whose sites' colors appear in the vertices of σ C . This is intentional and allows for a straightforward definition of the complex induced by this diagram, as shown below.

The discrete Riemannian Delaunay complex
We define the discrete Riemannian Delaunay complex as the set of simplices Del d g (P) = {σ | Vert(σ) ∈ P, ∩ p∈σ V d g (p) = 0}. Using a triangulation as canvas offers a very intuitive way to construct the discrete complex since each canvas k-simplex σ of C has k + 1 vertices {v 0 , . . . , v k } with respective colors {c 0 , . . . , c k } corresponding to the sites {p c0 , . . . , p c k } ∈ P. Due to the way discrete Voronoi cells overlap, a canvas simplex σ C belongs to each discrete Voronoi cell whose color appears in the vertices of σ. Therefore, the intersection of the discrete Voronoi cells {V d g (p i )} whose colors appear in the vertices of σ is non-empty and the simplex σ with vertices {p i } thus belongs to the discrete Riemannian Delaunay complex. In that case, we say that the canvas simplex σ C witnesses (or is a witness of) σ. For example, if the vertices of a canvas 3simplex τ C have colors yellow-blue-blue-yellow, then the intersection of the discrete Voronoi cells of the sites p yellow and p blue is non-empty and the one-simplex σ with vertices p yellow and p blue belongs to the discrete Riemannian Delaunay complex. The canvas simplex τ C thus witnesses the (abstract, for now) edge between p yellow and p blue . Figure 2 illustrates a canvas painted with discrete Voronoi cells, and the witnesses of the discrete Riemannian Delaunay complex. Figure 2: A canvas (black edges) and a discrete Riemannian Voronoi diagram drawn on it. The canvas simplices colored in red are witnesses of Voronoi vertices. The canvas simplices colored in grey are witnesses of Voronoi edges. Canvas simplices whose vertices all have the same color are colored with that color.
In other words, if a canvas simplex σ C witnesses a simplex σ, then for each face τ of σ, there exists a face τ C of σ C that witnesses τ . As we assume that there is no boundary, the complex is pure and it is sufficient to only consider canvas n-simplices whose vertices have all different colors to build Del d g (P). Similarly to the definition of curved and straight Riemannian Delaunay complexes, we can define their discrete counterparts we respectively denote by Del d g (P) and Del d g (P). We will now exhibit conditions such that these complexes are well-defined and embedded in Ω.

Equivalence between the discrete and the exact structures
We first give conditions such that Vor d g (P) and Vor g (P) have the same combinatorial structure, or, equivalently, that the dual Delaunay complexes Del g (P) and Del d g (P) are identical. Under these conditions, the fact that Del d g (P) is embedded in Ω will immediately follow from the fact that the exact Riemannian Delaunay complex Del g (P) is embedded (see Sections 3.1 and 3.2). It thus remains to exhibit conditions under which Del d g (P) and Del g (P) are identical. Requirements will be needed on both the set of sites in terms of density, sparsity and protection, and on the density of the canvas. The central idea in our analysis is that power protection of P will imply a lower bound on the distance separating two non-adjacent Voronoi objects (and in particular two Voronoi vertices). From this lower bound, we will obtain an upper bound on the size on the cells of the canvas so that the combinatorial structure of the discrete diagram is the same as that of the exact one. The density of the canvas is expressed by e C , the length of its longest edge.

9
The main result of this paper is the following theorem.
Theorem 5.1. Assume that P is a δ-power protected (ε, µ)-net in Ω with respect to g. Assume further that ε is sufficiently small and δ is sufficiently large compared to the distortion between g(p) and g in an ε-neighborhood of p. Let {λ i } be the eigenvalues of g(p) and 0 a value that depends on ε and δ (Precise bounds for ε, δ and l 0 are given in the proof).
The rest of the paper will be devoted to the proof of this theorem. Our analysis is divided into two parts. We first consider in Section 6 the most basic case of a domain of R n endowed with the Euclidean metric field. The result is given by Theorem 6.1. The assumptions are then relaxed and we consider the case of an arbitrary metric field over Ω in Section 7. As we shall see, the Euclidean case already contains most of the difficulties that arise during the proof and the extension to more complex settings will be deduced from the Euclidean case by bounding the distortion.

Equality of the Riemannian Delaunay complexes in the Euclidean setting
In this section, we restrict ourselves to the case where the metric field is the Euclidean metric g E . To simplify matters, we initially assume that geodesic distances are computed exactly on the canvas. The following theorem gives sufficient conditions to have equality of the complexes. Theorem 6.1. Assume that P is a δ-power protected (ε, µ)-net of Ω with respect to the Euclidean metric field g E . Denote by C the canvas, a triangulation with maximal edge length e C . If e C < min µ/16, δ 2 /64ε , then Del d E (P) = Del E (P). We shall now prove Theorem 6.1 by enforcing the two following conditions which, combined, give the equality between the discrete Riemannian Delaunay complex and the Riemannian Delaunay complex: (1) for every Voronoi vertex in the Riemannian Voronoi diagram v = ∩ {pi} V g (p i ), there exists at least one canvas simplex with the corresponding colors {c pi }; (2) no canvas simplex witnesses a simplex that does not belong to the Riemannian Delaunay complex (equivalently, no canvas simplex has vertices whose colors are those of non-adjacent Riemannian Voronoi cells).
Condition (2) is a consequence of the separation of Voronoi objects, which in turn follows from power protection. The separation of Voronoi objects has previously been studied, for example by Boissonnat et al. [2]. Although the philosophy is the same, our setting is slightly more difficult and the results using power protection are new and use a more geometrical approach (see Appendix C).

Sperner's lemma
Rephrasing Condition (1), we seek requirements on the density of the canvas C and on the nature of the point set P such that there exists at least one canvas n-simplex of C that has exactly the colors c 0 , . . . , c d of the vertices p 0 , . . . , p d of a simplex σ, for all σ ∈ Del g (P). To prove the existence of such a canvas simplex, we employ Sperner's lemma [28], which is a discrete analog of Brouwer's fixed point theorem. We recall this result in Theorem 6.2 and illustrate it in a twodimensional setting (inset).

Theorem 6.2 (Sperner's lemma).
Let σ = (p 0 , . . . , p n ) be an n-simplex and let T σ denote a triangulation of the simplex. Let each vertex v ∈ T σ be colored such that the following conditions are satisfied: • The vertices p i of σ all have different colors.
• If a vertex p lies on a k-face (p i0 , . . . p i k ) of σ, then p has the same color as one of the vertices of the face, that is p ij .
Then, there exists an odd number of simplices in T σ whose vertices are colored with all n + 1 colors. In particular, there must be at least one.
We shall apply Sperner's lemma to the canvas C and show that for every Voronoi vertex v in the Riemannian Voronoi diagram, we can find a subset C v of the canvas that fulfills the assumptions of Sperner's lemma, hence obtaining the existence of a canvas simplex in C v (and therefore in C) that witnesses σ v . Concretely, the subset C v is obtained in two steps: -We first apply a barycentric subdivision of the Riemannian Voronoi cells incident to v.
From the resulting set of simplices, we extract a triangulation T v composed of the simplices incident to v (Section 6.2).
-We then construct the subset C v by overlaying the border of T v and the canvas (Section 6.3).
We then show that if the canvas simplices are small enough -in terms of edge length -then C v is the triangulation of a simplex that satisfies the assumptions of Sperner's lemma.
The construction of C v is detailed in the following sections and illustrated in Figure 3: starting from a colored canvas (left), we subdivide the incident Voronoi cells of v to obtain T v (middle), and deduce the set of canvas simplices C v which forms a triangulation that satisfies the hypotheses of Sperner's lemma, thus giving the existence of a canvas simplex (in green, right) that witnesses the Voronoi vertex within the union of the simplices, and therefore in the canvas.
v Figure 3: Illustration of the construction of C v . The Riemannian Voronoi diagram is drawn with thick orange edges and the sites are colored squares. The canvas is drawn with thin gray edges and colored circular vertices. The middle frame shows the subdivision of the incident Voronoi cells with think black edges and the triangulation T v is drawn in yellow. On the right frame, the set of simplices C v is colored in purple (simplices that do not belong to C) and in dark yellow (simplices that belong to C).

The triangulation T v
For a given Voronoi vertex v in the Euclidean Voronoi diagram Vor E (P) of the domain Ω, the initial triangulation T v is obtained by applying a combinatorial barycentric subdivision of the Voronoi cells of Vor E (P) that are incident to v: to each Voronoi cell V incident to v, we associate to each face F of V a point c F in F which is not necessarily the geometric barycenter. We randomly associate to c F the color of any of the sites whose Voronoi cells intersect to give F . For example, in a two-dimensional setting, if the face F is a Voronoi edge that is the intersection of V red and V blue , then c F is colored either red or blue. Then, the subdivision of V is computed by associating to all possible sequences of faces  Figure 4 in dimension 3. As shall be proven in Lemma 6.3, T v can be used to define a combinatorial simplex that satisfies the assumptions of Sperner's lemma. T v as a triangulation of an n-simplex By construction, the triangulation T v is a triangulation of the (Euclidean) Delaunay simplex σ v dual of v as follows. We first perform the standard barycentric subdivision on this Delaunay simplex σ v . We then map the barycenter of a k-face τ of σ v to the point c Fi on the Voronoi face F i , where F i is the Voronoi dual of the k-face τ . This gives a piecewise linear homeomorphism from the Delaunay simplex σ v to the triangulation T v . We call the image of this map the simplex σ S and refer to the images of the faces of the Delaunay simplex as the faces of σ S . We can now apply Sperner's lemma. Lemma 6.3. Let P be a δ-power protected (ε, µ)-net. Let v be a Voronoi vertex in the Euclidean Voronoi diagram, Vor E (P), and let Σ v be defined as above. The simplex σ S and the triangulation T v satisfy the assumptions of Sperner's lemma in dimension n.
Proof. By the piecewise linear map that we have described above, T v is a triangulation of the simplex σ S . Because by construction the vertices c Fi lie on the Voronoi duals F i of the corresponding Delaunay face τ , c Fi has the one of the colors of of the Delaunay vertices of τ . Therefore, σ S satisfies the assumptions of Sperner's lemma and there exists an n-simplex in T v that witnesses v and its corresponding simplex σ v in Del g (P).

Building the triangulation C v
Let p i be the vertices of the k-face τ S of σ S . In this section we shall assume not only that τ S is contained in the union of the Voronoi cells of V (p i ), but in fact that τ S is a distance 8e C removed from the boundary of ∪V (p i ), where e C is the longest edge length of a simplex in the canvas. We will now construct a triangulation C v of σ S such that: • σ S and its triangulation C v satisfy the conditions of Sperner's lemma, • the simplices of C v that have no vertex that lies on the boundary ∂σ S are simplices of the canvas C.
The construction goes as follows. We first intersect the canvas C with σ S and consider the canvas simplices σ C,i such that the intersection of σ S and σ C,i is non-empty. These simplices σ C,i can be subdivided into two sets, namely those that lie entirely in the interior of σ S , which we denote by σ int C,i , and those that intersect the boundary, denoted by σ ∂ C,i . The simplices σ int C,i are added to the set C v . We intersect the simplices σ ∂ C,i with σ S and triangulate the intersection. Note that σ ∂ C,i ∩ σ S is a convex polyhedron and thus triangulating it is not a difficult task. The vertices of the simplices in the triangulation of σ ∂ C,i ∩ σ S are colored according to which Voronoi cell they belong to. Finally, the simplices in the triangulation of σ ∂ C,i ∩ σ S are added to the set C v .
Since T v is a triangulation of σ S , the set C v is by construction also a triangulation of σ S . This triangulation trivially gives a triangulation of the faces τ S . Because we assume that τ S is contained in the union of its Voronoi cells, with a margin of 8e C we now can draw two important conclusions: • The vertices of the triangulation of each face τ S have the colors of the vertices p i of τ S .
• None of the simplices in the triangulation of σ ∂ C,i ∩ σ S can have n + 1 colors, because every such simplex must be close to one face τ S , which means that it must be contained in the union of the Voronoi cells V (p i ) of the vertices of τ S .
We can now invoke Sperner's lemma; C v is a triangulation of the simplex σ S whose every face has been colored with the appropriate colors (since σ S triangulated by T v satisfies the assumptions of Sperner's lemma, see Lemma 6.3). This means that there is a simplex C v that is colored with n + 1 colors. Because of our second observation above, the simplex with these n + 1 colors must lie in the interior of σ S and is thus a canvas simplex.
We summarize by the following lemma: Lemma 6.4. If every face τ S of σ S with vertices p i is at distance 8e C from the boundary of the union of its Voronoi cells ∂(∪V (p i )), then there exists a canvas simplex in C v such that it is colored with the same vertices as the vertices of σ S .
The key task that we now face is to guarantee that faces τ S indeed lie well inside of the union of the appropriate Voronoi regions. This requires first and foremost power protection. Indeed, if a point set is power protected, the distance between a Voronoi vertex c and the Voronoi faces that are not incident to c, which we will refer to from now on as foreign Voronoi faces, can be bounded, as shown in the following Lemma: Lemma 6.5. Suppose that c is the circumcenter of a δ-power protected simplex σ of a Delaunay triangulation built from an ε-sample, then all foreign Voronoi faces are at least δ 2 /8ε far from c.
The proof of this Lemma is given in the appendix (Section C.2). In almost all cases, this result gives us the distance bound we require: we can assume that vertices {c F0 , c F1 , . . . , c Fn−1 , c Fn } which we used to construct T v , are well placed, meaning that there is a minimum distance between these vertices and foreign Voronoi objects. However it can still occur that foreign Voronoi objects are close to a face τ S of σ S . This occurs even in two dimensions, where a Voronoi vertex v can be very close to a face τ S because of obtuse angles, as illustrated in Figure 5. Thanks to power protection, we know that v is removed from foreign Voronoi objects. This means that we can deform σ S (in a piecewise linear manner) in a neighborhood of v such that the distance between v and all the faces of the deformed σ S is lower bounded.

Anisotropic triangulations via discrete Riemannian Voronoi diagrams
In general the deformation of σ S is performed by "radially pushing" simplices away from the foreign Voronoi faces of v with a ball of radius r = min µ/16, δ 2 /64ε . The value µ/16 is chosen so that we do not move any vertex of σ v (the dual of v): indeed, P is µ-separated and thus d E (p i , p j ) > µ. The value δ 2 /64ε is chosen so that σ S and its deformation stay isotopic (no "pinching" can happen), using Lemma 6.5. In fact it is advisable to use a piecewise linear version of "radial pushing", to ensure that the deformation of σ S is a polyhedron. This guarantees that we can triangulate the intersection, see Chapter 2 of Rourke and Sanderson [25]. After this deformation we can follow the steps we have given above to arrive at a well-colored simplex. Lemma 6.6. Let P be a δ-power protected (ε, µ)-net. Let v be a Voronoi vertex of the Euclidean Voronoi diagram Vor E (P), and T v as defined above. If the length e C of the longest canvas edge is bounded as follows: e C < r = min µ/16, δ 2 /64ε , then there exists a canvas simplex that witnesses v and the corresponding simplex σ v in Del E (P).
Conclusion So far, we have only proven that Del g (P) ⊆ Del d g (P). The other inclusion, which corresponds to Condition (2) mentioned above, is much simpler: as long as a canvas edge is shorter than the smallest distance between a Voronoi vertex and a foreign face of the Riemannian Voronoi diagram, then no canvas simplex can witness a simplex that is not in Del g (P). Such a bound is already given by Lemma 6.5 and thus, if e C < δ 2 /8ε then Del d g (P) ⊆ Del g (P). Observe that this requirement is weaker than the condition imposed in Lemma 6.6 and it was thus already satisfied. It follows that Del d g (P) = Del g (P) if e C < min µ/16, δ 2 /64ε , which concludes the proof of Theorem 6.1.

Remark 6.7.
Assuming that the point set is a δ-power protected (ε, µ)-net might seem like a strong assumption. However, it should be observed that any non-degenerate point set can be seen as a δ-power protected (ε, µ)-net, for a sufficiently large value of ε and sufficiently small values of δ and µ. Our results are therefore always applicable but the necessary canvas density increases as the quality of the point set worsens (Lemma 6.6). In our companion practical paper [26,Section 7], we showed how to generate δ-power protected (ε, µ)-nets for given values of ε, µ and δ.

Extension to more complex settings
In the previous section, we have placed ourselves in the setting of an (open) domain endowed with the Euclidean metric field. To prove Theorem 5.1, we need to generalize Theorem 6.1 to more general metrics, which will be done in the two following subsections.
The common path to prove Del d g (P) = Del g (P) in all settings is to assume that P is a power protected net with respect to the metric field. We then use the stability of entities under small metric perturbations to take us back to the now solved case of the domain Ω endowed with an Euclidean metric field. Separation and stability of Delaunay and Voronoi objects has previously been studied by Boissonnat et al. [2,5], but our work lives in a slightly more complicated setting. Moreover, our proofs are generally more geometrical and sometimes simpler. For completeness, the extensions of these results to our context are detailed in Appendix C for separation, and in Appendix E for stability.
We now detail the different intermediary settings. For completeness, the full proofs are included in the appendices.

Uniform metric field
We first consider the rather easy case of a non-Euclidean but uniform (constant) metric field over an (open) domain. The square root of a metric gives a linear transformation between the base space where distances are considered in the metric and a metric space where the Euclidean distance is used (see Appendix B.1.1). Additionally, we show that a δ-power protected (ε, µ)-net with respect to the uniform metric is, after transformation, still a δ-power protected (ε, µ)-net but with respect to the Euclidean setting (Lemma E.1 in Appendix E), bringing us back to the setting we have solved in Section 6. Bounds on the power protection, sampling and separation coefficients, and on the canvas edge length can then be obtained from the result for the Euclidean setting, using Theorem 6.6. These bounds can be transported back to the case of uniform metric fields by scaling these values according to the smallest eigenvalue of the metric (Theorem H.1 in Appendix H).

Arbitrary metric field
The case of an arbitrary metric field over Ω is handled by observing that an arbitrary metric field is locally well-approximated by a uniform metric field. It is then a matter of controlling the distortion.
We first show that, for any point p ∈ Ω, density separation and power protection are locally preserved in a neighborhood U p around p when the metric field g is approximated by the constant metric field g = g(p) (Lemmas E.2 and E. 16 in Appendix E): if P is a δ-power protected (ε, µ)net with respect to g, then P is a δ -power protected (ε , µ )-net with respect to g . Previous results can now be applied to obtain conditions on δ , ε , µ and on the (local) maximal length of the canvas such that Del d g (P) = Del g (P) (Lemma H.2 in Appendix H). These local triangulations can then be stitched together to form a triangulation embedded in Ω. The (global) bound on the maximal canvas edge length is given by the minimum of the local bounds, each computed through the results of the previous sections. This ends the proof of Theorem 5.1.
Once the equality between the complexes is obtained, conditions giving the embeddability of the discrete Karcher Delaunay triangulation and the discrete straight Delaunay triangulation are given by previous results that we have established in Sections 3.1 and 3.2 respectively. 15

Extensions of the main result
Approximate geodesic computations Approximate geodesic distance computations can be incorporated in the analysis of the previous section by observing that computing inaccurately geodesic distances in a domain Ω endowed with a metric field g can be seen as computing exactly geodesic distances in Ω with respect to a metric field g that is close to g (Section H.3 in Appendix H).

General manifolds
The previous section may also be generalized to an arbitrary smooth nmanifold M embedded in R m . We shall assume that, apart from the metric induced by the embedding of the domain in Euclidean space, there is a second metric g defined on M. Let π p : M → T p M be the orthogonal projection of points of M on the tangent space T p M at p. For a sufficiently small neighborhood U p ⊂ T p M, π p is a local diffeomorphism (see Niyogi [23]).
Denote by P Tp the point set {π p (p i ), p i ∈ P} and P Up the restriction of P Tp to U p . Assuming that the conditions of Niyogi et al. [23] are satisfied (which are simple density constraints on ε compared to the reach of the manifold), the pullback of the metric with the inverse projection . This implies immediately that if P is a δ-power protected (ε, µ)-net on M with respect to g then P Up is a δ-power protected (ε, µ)-net on U p . We have thus a metric on a subset of a n-dimensional space, in this case the tangent space, giving us a setting that we have already solved. It is left to translate the sizing field requirement from the tangent plane to the manifold M itself. Note that the transformation π p is completely independent of g. Boissonnat et al. [2,Lemma 3.7] give bounds on the metric distortion of the projection on the tangent space. This result allows to carry the canvas sizing field requirement from the tangent space to M.

Implementation
The construction of the discrete Riemannian Voronoi diagram and of the discrete Riemannian Delaunay complex has been implemented for n = 2, 3 and for surfaces of R 3 . An in-depth description of our structure and its construction as well as an empirical study can be found in our practical paper [26]. We simply make a few observations here.
The theoretical bounds on the canvas edge length provided by Theorems 5.1 and 6.1 are far from tight and thankfully do not need to be honored in practice. A canvas whose edge length are about a tenth of the distance between two seeds suffices. This creates nevertheless unnecessarily dense canvasses since the density does not in fact need to be equal everywhere at all points and even in all directions. This issue is resolved by the use of anisotropic canvasses.
Our analysis was based on the assumption that all canvas vertices are painted with the color of the closest site. In our implementation, we color the canvas using a multiple-front vector Dijkstra algorithm [7], which empirically does not suffer from the same convergence issues as the traditional Dijkstra algorithm, starting from all the sites. It should be noted that any geodesic distance computation method can be used, as long as it converges to the exact geodesic distance when the canvas becomes denser. The Riemannian Delaunay complex is built on the fly during the construction of the discrete Riemannian Voronoi diagram: when a canvas simplex is first fully colored, its combinatorial information is extracted and the corresponding simplex is added to Del g (P).

Overview of the appendices
The topics of the appendices are as follows: A We review basic definitions related to simplices and complexes.
B We describe our generalization of the notion of distortion between metrics due to Labelle and Shewchuk.
C We discuss the separation of Voronoi objects. The main differences between this appendix and [2] are in our definition of metric distortion, the use of power protection, and the more geometrical nature of the proofs. F We prove that the Delaunay simplices can be straightened, under sufficient conditions. The proof of the stability of the center of mass, which forms the core of this appendix, is also remarkable in the sense that it generalizes trivially to a far more general setting.
G We illustrate a degenerate case of Section 6.3.
H This appendix gives the proofs for our main result, Theorem 5.1, in the general setting of an arbitrary metric. We naturally rely heavily on Appendices C, D, E, and H.
The references for the Appendix can be found at the end of the Appendix.

A Simplices and complexes
The purpose of this section is to offer the precise definitions of concepts and notions related to simplicial complexes. The following definitions live within the context of abstract simplices and complexes.
A simplex σ is a non-empty finite set. The dimension of σ is given by dim σ = #(σ) − 1, and a j-simplex refers to a simplex of dimension j. The elements of σ are called the vertices of σ. The set of vertices of σ is noted Vert(σ).
If a simplex τ is a subset of σ , we say it is a face of σ , and we write τ ≤ σ. A 1-dimensional face is called an edge. If τ is a proper subset of σ , we say it is a proper face and we write τ < σ.
For any vertex p ∈ σ, the face opposite to p is the face determined by the other vertices of σ, and is denoted by σ p . If τ is a j-simplex, and p is not a vertex of σ, we may construct a (j + 1)-simplex σ = p * τ , called the join of p and τ as the simplex defined by p and the vertices of τ .
The length of an edge is the distance between its vertices. The height of p in σ is D(p, σ) = d(p, aff σ p ).
A circumscribing ball for a simplex σ is any n-dimensional ball that contains the vertices of σ on its boundary. If σ admits a circumscribing ball, then it has a circumcenter, C(σ), which is the center of the unique smallest circumscribing ball for σ. The radius of this ball is the circumradius of σ, denoted by R(σ).

A.1 Complexes
Before defining Delaunay triangulations, we introduce the more general concept of simplicial complexes. Since the standard definition of a simplex as the convex hull of a set of points does not extend well to the Riemannian setting (see Dyer [15]), we approach these definitions from a more abstract point of view.
The length of an edge is the distance between its vertices. A circumscribing ball for a simplex σ is any n-dimensional ball that contains the vertices of σ on its boundary. If σ admits a circumscribing ball, then it has a circumcenter, C(σ), which is the center of the unique smallest circumscribing ball for σ. The radius of this ball is the circumradius of σ, denoted by R(σ).
). The dihedral angle between two facets is the angle between their two supporting planes. If σ is a j-simplex with j ≥ 2, then for any two vertices p, q ∈ σ, the dihedral angle between σ p and σ q defines an equality between ratios of heights (see Figure 6).

A.2 Simplicial complexes
Simplicial complexes form the underlying framework of Delaunay triangulations. An abstract simplicial complex is a set K of simplices such that if σ ∈ K, then all the faces of σ also belong to K. The union of the vertices of all the simplices of K is called the vertex set of K. The dimension of a complex is the largest dimension of any of its simplices. A subset L ⊆ K is a subcomplex of K if it is also a complex. Two simplices are adjacent if they share a face and incident if one is a face of the other. If a simplex in K is not a face of a larger simplex, we say that it is maximal. If all the maximal simplices in a complex K of dimension n have dimension n, then the simplicial complex is pure. The star of a vertex p in a complex K is the subcomplex S formed by set of simplices that are incident to p. The link of a vertex p is the union of the simplices opposite of p in S p .
A geometric simplicial complex is an abstract simplicial complex whose faces are materialized, and to which the following condition is added: the intersection of any two faces of the complex is either empty or a face of the complex.

B Geodesic distortion
The concept of distortion was originally introduced by Labelle and Shewchuk [20] to relate distances with respect to two metrics, but this result can be (locally) extended to geodesic distances.

B.1 Original distortion
We first recall their definition and then show how to extend it to metric fields. The notion of metric transformation is required to define this original distortion, and we thus recall it now.

B.1.1 Metric transformation
Given a symmetric positive definite matrix G, we denote by F any matrix such that det(F ) > 0 and F t F = G. The matrix F is called a square root of G. The square root of a matrix is not uniquely defined: the Cholesky decomposition provides, for example, an upper triangular F , but it is also possible to consider the diagonalization of G as O T DO, where O is an orthonormal matrix and D is a diagonal matrix; the square root is then taken as The latter definition is more natural than other decompositions since √ D is canonically defined and F is then also a symmetric definite positive matrix, with the same eigenvectors as G. Regardless of the method chosen to compute the square root of a metric, we assume that it is fixed and identical for all the metrics considered.
The square root F offers a linear bijective transformation between the Euclidean space and a metric space, noting that where · stands for the Euclidean norm. Thus, the metric distance between x and y in Euclidean space can be seen as the Euclidean distance between two transformed points F x and F y living in the metric space of G.

B.1.2 Distortion
The distortion between two points p and q of Ω is defined by Labelle and Shewchuk [20] as A fundamental property of the distortion is to relate the two distances d Gp and d Gq . Specifically, for any pair x, y of points, we have Indeed, and the other inequality is obtained similarly.

B.2 Geodesic distortion
The previous definition can be defined to hold (locally) for metric fields instead of metric, as we show now.
Lemma B.1. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion, in the sense of Labelle and Shewchuk where d g and d g indicate the geodesic distances with respect to g and g respectively.
Proof. Recall that for p ∈ B g (p 0 , r 0 ) and, for any pair x, y of points, we have Therefore, for any curve γ(t) in U , we have that Considering the infimum over all paths γ that begin at x and end at y, we obtain the result.
Note that this result is independent from the definition of the distortion and is entirely based on the inequality comparing distances in two metrics (Equation 3).

C Separation of Voronoi objects
Power protected point sets were introduced to create quality bounds for the simplices of Delaunay triangulations built using such point sets [2]. We will show that power protection allows to deduce additional useful results for Voronoi diagrams. In this section we show that when a Voronoi diagram is built using a power protected sample set, its non-adjacent Voronoi faces, and specifically its Voronoi vertices are separated. This result is essential to our proofs in Sections 6 and 7 where we approximate complicated Voronoi cells with simpler Voronoi cells without creating inversions in the dual, which is only possible because we know that Voronoi vertices are far from one another.
We assume that the protection parameter δ is proportional to the sampling parameter ε, thus there exists a positive ι, with ι ≤ 1, such that δ = ιε. We assume that the separation parameter µ is proportional to the sampling parameter ε and thus there exists a positive λ, with λ ≤ 1, such that µ = λε.

C.1 Separation of Voronoi vertices
The main result provides a bound on the Euclidean distance between Voronoi vertices of the Euclidean Voronoi diagram of a point set and is given by Lemma C.4. The following three lemmas are intermediary results needed to prove Lemma C.4.
The altitude of the vertex p i in the simplex σ is denoted by D(p i , σ). Proof. The spheres ∂B and ∂B +δ intersect in a (n − 2)-sphere S which is contained in a hyperplane H parallel to the hyperplane H = aff(τ ). For anyq ∈ S we have where the last equality follows from Lemma C.2 and d( H, H) denotes the distance between the two parallel hyperplanes. See Figure 7 for a sketch. Since p ∈ ∂B, p belongs to B(σ ) +δ if and only if p lies in the strip bounded by H and H, which is equivalent to The result now follows.
We can make this bound independent of the simplices considered, as shown in Lemma C.4.
Lemma C.4. Let P be a δ-power protected (ε, µ)-sample. The protection parameter ι is given by δ = ιε. For any two adjacent Voronoi vertices c and c of VD(P), we have Proof. For any simplex σ, we have D(p, σ) ≤ 2R σ for all p ∈ σ, where R σ denotes the radius of the circumsphere of σ. For any σ in the triangulation of an ε-net, we have R σ ≤ ε. Thus D(p, σ) ≤ 2ε, and Lemma C.3 yields c − c ≥ δ 2 /4ε.
Remark C.5. In this section, we have computed a lower bound on the distance between two (adjacent) Voronoi vertices. In Appendix E, we shall show that Voronoi vertices are stable under metric perturbations, meaning that when a metric field is slightly modified, the position of a Voronoi vertex does not move too much. The combination of this separation and stability will then be the basis of many proofs in this paper.

C.2 Separation of Voronoi faces (Proof of Lemma 6.5)
Similar results can be obtained on the distance between a Voronoi vertex c and faces that are not incident to c, also referred to as foreign faces. Note that we are still in the context of an Euclidean metric.
The following lemma can be found in [3, Lemma 3.3] for ordinary protection instead of power protection.  Proof. We denote by p 0 an (arbitrary) vertex of σ, and by q a vertex that is not in σ but is adjacent to p 0 in Del(P). Let x be a point in B(c, r) ∩ V E (p 0 ), with 0 < r < δ 2 /4ε. The upper bound for r is chosen with Lemma C.3 in mind: we are trying to find a condition such that x ∈ V E (p 0 ). The notations are illustrated in Figure 8. Because of the triangle inequality, we have that By power protection, we have that d(c, q) 2 ≥ d(c, p i ) 2 + δ 2 . Therefore, Without loss of generality, we can assume that q is the site closest to c and thus d( This implies that as long r < δ 2 /8ε, x lies in a Voronoi object associated to the vertices p i of σ.
Further progressing, we can show that Voronoi faces are thick, with Lemma C.7. This property is useful to construct a triangulation that satisfies the hypotheses of Sperner's lemma.

Inria
Lemma C.7. Let P be a δ-power protected (ε, µ)-net. Let V 0 be the Voronoi cell of the site p 0 ∈ P in the Euclidean Voronoi diagram VD E (P). Then for any k-face F 0 of V 0 with k ∈ [1, . . . , n], there exists x ∈ F 0 such that where ∂F 0 denotes the boundary of the face F 0 .
Proof. All the vertices of F 0 are circumcenters of VD E (P). Consider the erosion of the face F 0 by a ball of radius δ 2 16ε and denote it F − 0 . If F − 0 is empty, we contradict the separation Lemma 6.5.

D Bounds on dihedral angles
The use of nets allows us to deduce bounds on the dihedral angles of faces of the Delaunay triangulation, as well as on the dihedral angles between adjacent faces of a Voronoi diagram. Those bounds are frequently used throughout this paper, and specifically to prove stability of Voronoi vertices (see Appendix E). Since we are interested in dihedral bounds in the Euclidean setting, the point set is first assumed to be a net with respect to the Euclidean metric field. We complicate matters slightly with Lemma D.5 by assuming that the point set is a power protected net with respect to an arbitrary metric field that is not too far from the Euclidean metric field (in terms of distortion), and still manage to expose bounds with respect to the Euclidean distance.

D.1 Bounds on the dihedral angles of Euclidean Voronoi cells
Assuming that a point set is an (ε, µ)-net allows us to deduce lower and upper bounds on the dihedral angles between adjacent Voronoi faces when the metric field is Euclidean.

Lemma D.1.
Let Ω = R n and P be an (ε, µ)-net with respect to the Euclidean distance on Ω. Let p ∈ P and V(p) be the Voronoi cell of p ∈ P. Let q, r ∈ P be two sites such that V(q) and V(r) are adjacent to V(p) in the Euclidean Voronoi diagram of P. Let θ be the dihedral angle between BS(p, q) and BS(p, r). Then Proof. We consider the plane H that passes through the sites p, q and r. Notations used below are illustrated in Figure 9. Lower bound. Let m pq and m pr be the projections of the site p on respectively the bisectors BS(p, q) and BS(p, r). Since P is an (ε, µ)-net, we have that l q = |p − m pq | ≥ µ/2, l r = |p − m pr | ≥ µ/2 and L = |p − c| ≤ ε. Thus Note that since 0 < µ < 2ε, we have 0 < µ/2ε < 1.

Upper bound.
To obtain an upper bound on θ, we compute a lower bound on the angle α = qpr at p, noting that θ = π − α. Let l qr = |q − r|, and R = |c − r|. By the law of sines, we have  Since P is an (ε, µ)-net, we have l qr ≥ µ and R ≤ ε. Finally,

D.2 Bounds on the dihedral angles of Euclidean Delaunay simplices
Bounds on the dihedral angles of simplices guarantee the thickness -the smallest height of any vertex -of simplices, and thus their quality. Additionally, they can be used to show that circumcenters of adjacent simplices are far from one another, thus proving the stability of circumcenters and of Delaunay simplices.

D.2.1 Using power protection with respect to the Euclidean metric field
We first assume that the metric field g is the Euclidean metric field g E and show that the simplices of an Euclidean Delaunay triangulation constructed from a power protected net are thick.
Recall that the dihedral angle can be expressed through heights as The bound on dihedral angles is obtained by bounding the height of vertices in a simplex. An obvious upper bound on the height of a vertex p in σ is D(p, σ) < 2ε. A lower bound is already obtained in Lemma C.3: we have that D(p, σ) ≥ δ 2 /2 c − c = δ 2 /4ε. We can thus bound the dihedral angles as follows: Lemma D.2. Let P be a δ-power protected (ε, µ)-net with respect to the Euclidean metric field g 0 . Let ϕ be the dihedral angle between two facets τ 1 and τ 2 of a simplex σ of Del g0 (P). Then arcsin(s 0 ) ≤ ϕ ≤ π − arcsin(s 0 ), and A λ,ι defined as in the previous Lemma.

D.2.2 Using power protection with respect to an arbitrary metric field
When considering a Voronoi diagram built using the geodesic distance induced by an arbitrary metric field g, the assumption of a power protected net is made with respect to this geodesic distance. To prove the stability of the power protected assumption under metric perturbation, we will however need to deduce lower and upper bounds on the dihedral angles between faces of the simplices of the Riemannian Delaunay complex with respect to the Euclidean metric field . We prove here that if the point set P is a δ-power protected (ε, µ)-net with respect to an arbitrary metric field g and if the distortion between g and the Euclidean metric field g E is small, then the dihedral angles of the simplices of the Euclidean Delaunay triangulation of P can be bounded.
We first give a result on the stability of Delaunay balls which expresses that if two metric fields have low distortion, the Delaunay balls of a simplex with respect to each metric field are close. One of these metric fields is assumed to be the Euclidean metric field. A similar result can be found in the proof of Lemma 5 in the theoretical analysis of locally uniform anisotropic meshes of Boissonnat et al. [6].

Lemma D.3. Let U
⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Assume furthermore that g is the Euclidean distance (thus d g = d E ). Let B = B g (c, r) be the geodesic ball with respect to the metric field g, centered on c ∈ U and of radius r. Assume that B E (c, ψ 0 r) ⊂ U . Then B can be encompassed by two Euclidean balls B E (c, r −ψ0 ) and B E (c, r +ψ0 ) with r −ψ0 = r/ψ 0 and r +ψ0 = ψ 0 r.
Proof. This is a straight consequence from Lemma B.1. Indeed, we have for all x ∈ U that and similarly 1 Thus, On the other hand, we have giving us B ⊂ B E (c, r +ψ0 ).
Note that r −ψ0 and r +ψ0 go to r as ψ 0 goes to 1. We now use this stability result to provide Euclidean dihedral angle bounds assuming power protection with an arbitrary metric field that is close to g E . We first require the intermediary result given by Lemma D.4.

Lemma D.4 (Whitney's lemma).
Let H be a hyperplane in Euclidean n-space and τ an n − 1simplex whose vertices lie at most η from the H and whose minimum height is h min . Then the angle ξ between aff(τ ) and H is bounded from above by Proof. By definition, the barycenter of a (n − 1)-simplex has barycentric coordinates λ i = 1/d. This means that it has distance a h min /n to each of its faces. So the ball centered on the barycenter with radius h min /n is contained in the simplex. This means that for any direction in aff(τ ) there exists a line segment of length 2h min /n that lies within τ . Moreover the end points of this line segments lie at most η from H because the vertices of the τ do. This means that the angle ξ is bounded by We can now give the main result which bounds Euclidean dihedral angles, assuming power protection with respect to the arbitrary metric field. Lemma D.5. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U be a δ-power protected (ε, µ)-net over U , with respect to g. Let B = B g (c, r) and B = B g (c , r ) be the geodesic Delaunay balls of σ = τ * p and σ = τ * p , with σ, σ ∈ Del g (P U ). Assume that P U is sufficiently dense such that U contains B and B . Then the minimum height of the simplex satisfies

Anisotropic triangulations via discrete Riemannian Voronoi diagrams 29
Note that this is the height of p in σ with respect to the Euclidean metric.
We proceed in a similar fashion as the proof Lemma C.1. However, a significant difference is that we are interested here in only proving that power protection with respect to the generic metric field provides a height bound in the Euclidean setting (rather than an equivalence).
Proof. We use the following notations, illustrated in Figure 10: are the two sets of (Euclidean) enclosing balls of respectively B and B defined as in Lemma D.3.
• q is a point on S,q is a point on S and q is a point on S .
• H is the geodesic supporting plane of τ , that is {argmin( vi∈τ λ i d g (x, v i ))}.
• H E , H E and H E are the two Euclidean hyperplanes orthogonal to [cc ] passing through respectively q,q and q .
The separation between the hyperplanes H E and H E provides a lower bound on the distance d E (p, τ ), thus on the height D E (p, σ). We therefore seek to bound d E (H E , H E ). By definition of the enclosing Euclidean balls, we have that |cq| = r ψ 0 , |cq| = ψ 0 r, |c q| = ψ 0 r , and |c q| = 1 ψ 0 r 2 + δ 2 Using the law of cosines in the triangles cc q and cc q, we find where λ c r = |c − c |. Subtracting one from the other, we obtain Similarly we can calculate the distance d E (H E , H E ) to be: Lemma D.4 gives us that the angle ξ between H E and aff(τ ) is bounded by The vertices of τ lie in between H E and H E and inside B E (c , ψr ). Thus, if we restrict to thẽ qcc plane, distance between the point where aff(τ ) intersects H E and the orthogonal projection π H E (q) ofq on H E is at most ψ(r + r ). This in turn implies that the line connecting π H E (q) and q intersects aff(τ ) at most distance (r + r ) tan(ξ) from H E . This also gives us that the distance Inria Anisotropic triangulations via discrete Riemannian Voronoi diagrams 31 fromq to its orthogonal projection π aff(τ ) (q) on aff(τ ) is bounded by Here we note that h min originally referred to the minimum height of a face, but because the height of a face is always greater than the height in the simplex we may read this in a general way, that is we regard h min as a universal lower bound on the height. Because |q − π aff(τ ) (q)| bounds the height of the simplex we get the following relation: To make the expression a bit more readable, we introduce Multiplying with h min and squaring we find: We note that in the limit ψ → 1, s 1 tends to zero, we therefore expand h min is terms of s 1 . Using a computer algebra system we find We emphasize that this equation gives h min = δ 2 /4ε as ψ tends to 1. This means that for sufficiently small metric distortion the height of a protected simplex will be strictly positive.
Lemma D.6. Let g be a metric field and P be a point set defined over R n . Let ψ 0 = ψ(g, E). Assume that P is a δ-power protected (ε, µ)-net over R n . Let φ be the dihedral angle between two faces of a simplex τ ∈ Del E (P). Then with s 0 detailed in the proof.
Proof. Denote h min the lower bound on D(q, σ) obtained in the previous lemma. We also immediate have that Let ϕ be the dihedral angle between aff(σ p ) and aff(σ q ). Recall that sin(ϕ) = sin ∠(aff(σ p ), aff(σ q )) = D(p, σ) For s 0 to make sense, we want s 0 > 0, which is bound to happen as ψ 0 goes to 1, as shown in Lemma D.5: h min goes to δ 2 /4ε thus h min /4ε goes to ι 2 /4 > 0.

E Stability
The notion of stability designates the conservation of a property despite changes of other parameters. In our context, the main assumptions concern the nature of point sets: we assume that point sets are power protected nets and wish to preserve these hypotheses despite (small) metric perturbations. The stability is important both from a theoretical and a practical point of view. Indeed, if an assumption is stable under perturbation, we can simplify matters without losing information. For example, we will prove that if a point set is a net with respect to a metric field g, then it is also a net (albeit with slightly different sampling and separation parameters) for a metric field g that is close to g (in terms of distortion) In a practical context, the stability of an assumption provides robustness with respect to numerical errors (see, for example, the work of Funke et al. [17] on the stability of Delaunay simplices).

E.1 Stability of the protected net hypothesis under metric transformation
It is rather immediate that the power protected net property is preserved when the point set is transformed (see Section B.1.1) between these spaces, as shown by the following lemma.
Lemma E.1. Let P be a δ-power protected (ε, µ)-net in Ω. Let g = g 0 be a uniform Riemannian metric field and F 0 a square root of g 0 . If P is a δ-power protected (ε, µ)-net with respect to g 0 then P = {F 0 p i , p i ∈ P} is a δ-power protected (ε, µ)-net with respect to the Euclidean metric.
Proof. This results directly from the observation that

E.2 Stability of the protected net hypothesis under metric perturbation
Metric field perturbations are small modifications of a metric field in terms of distortion. Since a generic Riemannian metric field g is difficult to study, we will generally consider a uniform approximation g 0 of g within a small neighborhood, such that the distortion between both metric fields is small over that neighborhood. In that context, the perturbation of g is the act of bringing g onto g 0 . Stability of the assumption of power protection was previously investigated by Boissonnat et al. [4] in the context of manifold reconstructions.

E.2.1 Stability of the net property
The following lemma shows that the "net" property is preserved when the metric field is perturbed: a point set that is a net with g is also a net with respect to g , a metric field that is close to g.

Lemma E.2. Let U
⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U be a point set in U . Suppose that P U is an (ε, µ)-net of U with respect to g. Then P U is an (ε g , µ g )-net of U with respect to g with ε g = ψ 0 ε and µ g = µ/ψ 0 .

E.2.2 Stability of the power protection property
It is more complex to show that the assumption of power protection is preserved under metric perturbation. Previously, we have only considered two similar but arbitrary Riemannian metric fields g and g on a neighborhood U . We now restrict ourselves to the case where g is a uniform metric field. We shall always compare the metric field g in a neighborhood U with the uniform metric field g = g 0 = g(p 0 ) where p 0 ∈ U . Because g 0 and the Euclidean metric field differ only by a linear transformation, we simplify matters and assume that g 0 is the Euclidean metric field.
We now give conditions such that the point set is also protected with respect to g 0 . A few intermediary steps are needed to prove the main result: • Given two sites, we prove that the bisectors of these two sites in the Voronoi diagrams built with respect to g and with respect to g = g E are close (Lemma E.5).
• We prove that the Voronoi cell of a point p 0 with respect to g, V g (p 0 ) can be encompassed by two scaled versions (one larger and one smaller) of the Euclidean Voronoi cell V E (p 0 ) (Lemma E.12).
• We combine this encompassing with bounds on the dihedral angles of Euclidean Delaunay simplices given by Lemma D.6 to compute a stability region around Voronoi vertices where the same combinatorial Voronoi vertex of lives for both V g (p 0 ) and V E (p 0 ) in 2D (Lemma E.13). We then extend it to any dimension by induction (Lemma E.14).
The main result appears in Lemma E.16 and gives the stability of power protection under metric perturbation. We first define the scaled version of a Voronoi cell more rigorously.

Definition E.4 (Relaxed Voronoi cell).
Let ω ∈ R. The relaxed Voronoi cell of the site p 0 is The following lemma expresses that two Voronoi cells computed in similar metric fields are close.
Lemma E.5. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B g (p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a point set in U . Let V p0,g denote a Voronoi cell with respect to the Riemannian metric field g.
Proof. Let BS g (p 0 , p i ) be the bisector between p 0 and p i with respect to g. Let y ∈ BS g (p 0 , p i ) ∩ B g (p 0 , ρ), where B g (p 0 , ρ) denotes the ball centered at p 0 of radius ρ with respect to g . Now d g (y, p 0 ) = d g (y, p i ), and thus Thus d g (y, p 0 ) 2 ≤ d g (y, p i ) 2 + ω and d g (y, p 0 ) 2 ≥ d g (y, p i ) 2 − ω with ω = 2ρ 2 (ψ 2 0 − 1), which gives us the expected result.

Inria
Remark E.6. Lemma E.5 does not require g to be a uniform metric field.
We clarify in the next lemma that the bisectors of a Voronoi diagram with respect to a uniform are affine hyperplanes.
Lemma E.7. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a point set in U . Let g be a uniform metric field. We refer to g as g 0 to emphasis its constancy. Let p 0 ∈ P U . The bisectors of V ±ω0 g0 (p 0 ) are hyperplanes.
Proof. The bisectors of V ±ω0 g0 (p 0 ) are given by which is the equation of an hyperplane since g 0 is uniform.
The cells V ±ω0 g0 (p 0 ) are unfortunately impractical to manipulate as we do not have an explicit distance between the boundaries ∂ V g0 (p 0 ) and ∂ V ±ω0 g0 (p 0 ). However that distance can be bounded; this is the purpose of the following lemma.
Lemma E.8. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a point set in U . Assume furthermore that g 0 is the Euclidean metric field. Let p 0 ∈ P U . We have Proof. Let m ω0 be the intersection of the segment [p 0 , p i ] and the bisector BS −ω0 g0 (p 0 , p i ), for i = 0. Let m be the intersection of the segment [p 0 , p i ] and ∂ V g0 (p 0 ), for i = 0. We have Since (m − m ω0 ) and (p i − p 0 ) are linearly dependent, P is µ 0 -separated, which implies that Definition E.9. In the following, we use ρ = 2ε 0 , and therefore ω 0 = 8ε 2 0 (ψ 2 0 − 1). We show that this choice is reasonable in Lemma E. 15. Additionally, we define We are now ready to encompass the Riemannian Voronoi cell of p 0 with respect to an arbitrary metric field g with two scaled versions of the Euclidean Voronoi cell of p 0 . The notions of dilated and eroded Voronoi cells will serve the purpose of defining precisely these scaled cells.

Definition E.10 (Eroded Voronoi cell).
Let ω ∈ R. The eroded Voronoi cell of p 0 is the morphological erosion of V g (p 0 ) by a ball of radius ω:

Definition E.11 (Dilated Voronoi cell).
Let ω ∈ R. The dilated Voronoi cell of p 0 is: where H ω (p 0 , p i ) is the half-space containing p 0 and delimited by the bisector BS(p 0 , p i ) translated away from p 0 by ω 0 (see Figure 11).
The second important step on our path towards the stability of power protection is the encompassing of the Riemannian Voronoi cell, and is detailed below.
Lemma E.12. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a point set in U . We have These inclusions are illustrated in Figure 11.
Proof. Using the notations and the result of Lemma E.8, we have Since the bisectors BS ω0 g (p 0 , p i ) are hyperplanes, the result follows.
In Figures 11 and 12, DV +η0 g0 and EV −η0 g0 are shown in green and V g0 in yellow. The next step is to prove that we have stability of the Voronoi vertices of V g (p 0 ), meaning that a small perturbation of the metric only creates a small displacement of any Voronoi vertex of the Voronoi cell of p 0 . The following lemma shows that Voronoi vertices are close if the distortion between the metric fields g and g 0 is small. The approach is to use Lemmas E.5 and E.12: we know that each (n − 1)-face of the Riemannian Voronoi cell V g (p 0 ) and shared by a Voronoi cell V g (q) is enclosed within the bisectors of DV g0 (p 0 ) and EV g0 (p 0 ) (which are translations of the bisectors of V g0 (p 0 )) for the sites p 0 and q. These two bisectors create a "band" that contains the bisector BS g (p 0 , q). Given a Voronoi vertex c in V g (p 0 ), c can be obtained as the intersection of n + 1 Voronoi cells, but also as the intersection of n Voronoi (n − 1)-faces of V g (p 0 ). The intersection of the bands associated to those n (n − 1)−faces is a parallelotopic-shaped region which by definition contains the same (combinatorially speaking) Voronoi vertex, but for V g0 (p 0 ). Lemmas E.13 and E.14 express this reasoning, which is illustrated in Figure 12 for 2D and 13 for any dimension.
Lemma E.13. We consider here Ω = R 2 . Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a point set in U . Assume furthermore that P U is a δ 0 -power protected (ε 0 , µ 0 )-net with respect to g 0 (the Euclidean metric). Let p 0 ∈ P U . Let c and c 0 be the same Voronoi vertex in respectively V g (p 0 ) and V g0 (p 0 ). Then d g0 (c, c 0 ) ≤ χ 2 with where λ is given by Proof. We use Lemma E.12. V g (p 0 ) lies in DV +η0 g0 and contains EV −η0 g0 . The circumcenters c and c 0 lie in a parallelogrammatic region centered on c 0 , itself included in the ball centered on c 0 and with radius χ. The radius χ is given by half the length of the longest diagonal of the parallelogram (see Figure 12). By Lemma E.2, P is an (ε 0 , µ 0 )-net with respect to g 0 . Let θ be the angle of the Voronoi corner of V g0 (p 0 ) at c 0 . By Lemma D.1, that angle is bounded: Since π − θ M < θ m , χ is maximal when θ > π/2. We thus assume θ > π/2, and compute a bound on χ as follows: using sin( 1 2 arcsin(α)) = 1 Figure 12: Black lines trace Vor g0 and cyan lines trace Vor g . The cell V g0 (p 0 ) is colored in yellow and the cell V g (p 0 ) is dashed. The green regions correspond to DV +η g0 (p 0 ) and EV −η g0 (p 0 ). On the left, the configuration where θ < π/2; on the right, θ > π/2. The result obtained in Lemma E.13 can be extended to any dimension using induction and the stability of the Voronoi vertices of facets.
Lemma E.14. We consider here Ω = R n . Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Let P U = {p i } be a δ 0 -power protected (ε 0 , µ 0 )-net with respect to g 0 (the Euclidean metric). Let p 0 ∈ P U . Let c and c 0 be the same Voronoi vertex in respectively V g (p 0 ) and V g0 (p 0 ). Then where χ 2 is defined as in Lemma E.13, and ϕ 0 is the maximal dihedral angle between two faces of a simplex.
Proof. We know from Lemma E.5 that V g (p 0 ) lies in DV +η0 g0 and contains EV −η0 g0 . The circumcenters c and c 0 lie in a parallelotopic region centered on c 0 defined by the intersection of n Euclidean thickened Voronoi faces. This parallelotope and its circumscribing sphere are difficult to compute. However, it can be seen as the intersection of two parallelotopic tubes defined by the intersection of n − 1 Euclidean thickened Voronoi faces. From another point of view, this is Inria the computation of the intersection of the thickened duals of two facets τ 1 and τ 2 incident to p 0 of the simplex σ ∈ Del(P), dual of c (and c 0 ), see Figure 13 (left). Figure 13: Left, a simplex, the duals of two faces and their respective thickened duals (cylinders); the orange thickened dual is orthogonal to the purple face (and inversely). Right, the intersection of the two tubes, seen from above, illustrates the proof of Lemma E.14.
The stability radius χ is computed incrementally by increasing the dimension and proving stability of the circumcenters of the faces of the simplices. We prove the formula by induction.
The radius of each tube is given by the stability of the radius of the circumcenter in the lower dimension of the facet. The base case, R n = R 3 , is solved by Lemma E.13, and gives We now consider two facets τ 1 and τ 2 that are incident to p 0 . Denote D 1 and D 2 their respective duals. By Lemma E.5, D 1 and D 2 lie in two cylinders C 1 and C 2 of radius χ 2 . C 1 and C 2 are also orthogonal to τ 1 and τ 2 and c and c 0 lie in C 1 ∩ C 2 . The angle ϕ between C 1 and C 2 is exactly the dihedral angle between τ 1 and τ 2 . By Lemma D.6, we have arcsin(s 0 ) ≤ ϕ ≤ π − arcsin(s 0 ) with s 0 = 1 2 Let ϕ 0 = arcsin(s 0 ) = π − arcsin(s 0 ). We encompass the intersection of the cylinders, difficult to compute, with a sphere whose radius can be computed as follows (see Figure 13, right): Thus, We have assumed in different lemmas that we could pick values of ρ 0 or ω 0 that fit our need. The following lemma shows that these assumptions were reasonable. Lemma E.15. Lemmas E.12 and E.14 allow us to characterize the parameter ρ 0 more precisely. Indeed, an assumption of Lemma E.5 was that V ω0 g0 (p 0 ) is included in a ball B g0 (p 0 , ρ 0 ). If the sampling of P is sufficiently dense, such an assumption is reasonable.
Proof. By definition, the Voronoi cell V ω0 g0 (p 0 ) is included in the dilated cell DV +η0 g0 (p 0 ). Since the point set is an ε-sample, we have d g0 (p 0 , x) ≤ ε 0 for x ∈ V g0 (p 0 ). By Lemma E.14, we have for Recall from Lemma E.14 that with η 0 = ρ 2 0 (ψ 2 0 −1) µ0 The parameter ρ 0 can be chosen arbitrarily as long as it is greater than ε 0 , and we have taken ρ 0 = 2ε 0 (see Definition E.9), which imposes Recall that the parameter λ is fixed. By continuity of the metric field, lim ε→0 ψ 0 = 1, therefore the left hand side goes to 0 and Inequality (6) is eventually satisfied as the sampling is made denser.
Finally, we can now show the main result: the power protection property is preserved when the metric field is perturbed.
Lemma E. 16. Let U ⊂ Ω be open, and g and g be two Riemannian metric fields on U . Let ψ 0 ≥ 1 be a bound on the metric distortion. Suppose that U is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + , such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g (p)) ≤ ψ 0 . Assume that P U is a δ-power protected (ε, µ)-net in U with respect to g. If δ is well chosen, then P U is a δ 0 -power protected net with respect to g 0 , with Proof. By Lemma E.2, we know that P U is (ε 0 , µ 0 )-net with respect to g 0 . Let q ∈ P U , with q not a vertex in the dual of c. Let c 0 be the combinatorial equivalent of c in V g0 (P) Since P is a δ-power protected net with respect to g, we have d g (c, q) > √ r 2 + δ 2 , where r = d g (c, p). On the one hand, we have by Lemma E.13. On the other hand, for any p ∈ P U such that p is a vertex of the dual of c, we have Thus δ 0 -power protection of P U with respect to g 0 requires

This is verified if
. This gives us This condition on δ is only reasonable if the right hand side is not too large. Indeed, since P is an ε-sample, we must have d g (c, q) < 2ε. However, we have that d g (c, q) 2 > d g (c, p) 2 + δ 2 by δ-power protection of P with respect to g. Because d g (c, p) < ε, it suffices that δ < ε. We will now show that this is reasonable by examining the limit of the right hand side of Inequality (7).
We note, see Lemma E.14, that where ε 0 = ψ 0 ε and µ 0 = µ/ψ 0 = λε/ψ 0 (see Remark E.3). So that This means that the right hand side of (7) is of the form f (ψ 0 )(ψ 2 0 − 1)ε 2 , where f (ψ 0 ) is a function that tends to a constant as ψ 0 goes to 1: So the bound given in Equality (7) may be easily satisfied if the metric distortion is sufficiently small. We now provide an explicit value for δ 0 in terms of δ. Let ξ = d g (c, q) and ξ 0 = d g0 (c 0 , q 0 ). We have the following bounds on r 0 and ξ 0 : If we hadδ-power protection, we would have Therefore we can take δ 2 0 = δ 2 . Note that with this definition, δ 0 goes to δ as ψ 0 goes to 1, which proves that our value of δ 0 is legitimate.

F Embeddability of the straight Delaunay triangulation (Proofs of Section 3.2)
We first prove Lemma 3.3, recalled below, which bounds the distance between the same point on a the Karcher and the straight simplex.

Inria
Lemma. Let P be an ε-sample with respect to g on Ω. Let {p i } be a set of n + 1 vertices in M such that N = ∩ pi∈P V (p i ) = ∅. Letσ and σ be the straight and Karcher simplices that realize N . Let x be a point on the Karcher simplex σ determined by the barycentric coordinates {λ i } (see Equation 2). Let x e be the point uniquely determined by Proof. The key observation is that given a convex function f and a function f that is close, that is f − f < α with α small, then the minimum value of f is at most of min f + α. If we observe that at any point x where f (x) > min f + 2α, we also have f (x) > min f + α so x is not a minimum of f , we see that the minima of f and f can not be far apart. In particular, we have that if x f ,min is the point where f attains its minimum, then f (x f ,min ) ≤ min f + α. The precise argument goes as follows.
We again assume that (possibly after a linear transformation) the metric is close to the Euclidean one, that is: If we assume that |x − y| ≤ 4ε and ψ 0 ≤ 2, it follows that attains its minimum.
Using the bounds above, we find that We also see that Taking f to be i λ i d g ( x, p i ) 2 in the explanation above, we find that because the Euclidean barycenter x e = i λ i p i is where the function i λ i | x − p i | 2 attains its minimum. Combining these results yields which yields a distance bound of 2 · 4 3 (ψ 0 − 1)ε 2 .
Although we have formulated this metric distortion result for simplices, the same proof extends almost verbatim to continuous distributions. By this we mean that the barycenter with respect to a metric g of a continuous distribution is close to the barycenter with respect to the Euclidean metric, if g is close to the Euclidean metric. Furthermore, note that the proof does not depend on the weights being positive.
We now prove Theorem 3.4, recalled below.
Theorem. Let P be a δ-power protected (ε, µ)-net with respect to g on Ω. Proof. The lower bound on D(p, σ) given in Appendix D is proportional to ε. The proximity (upper) bound from Lemma 3.3 is proportional to (ψ 0 − 1)ε, therefore going to 0 much faster. The embeddability is thus satisfied once

G Deforming the triangulation T v
An extreme configuration can have a sphere in Figure 14 separating parts of T v from the Voronoi vertex v. In that case, we can chose different spheres to "push away" faces. Figure 14 shows the construction. We do not detail the computations.

H Equality of the Riemannian Delaunay complexes in advanced geometric settings
In this section, we explain precisely how to obtain conditions on the maximal length of a canvas edge and the quality of the sample set such that the Del d g (P) and the Del g (P) are equal.

H.1 Uniform metric field
We first investigate the setting of a subdomain of R n endowed with a uniform metric field. If, for all σ ∈ Del d g (P), there exists σ C ∈ C such that σ C witnesses σ, and ∀σ C ∈ C the simplex witnessed by σ C belongs to Del d g (P), then we say that C captures Del d g (P). A Voronoi cell is said to be captured if all its Voronoi vertices are witnessed by the canvas.
By Lemma E.1, if a point set P in R n is a δ-power protected (ε, µ)-net with respect to a uniform metric field g 0 then the point set P = {F 0 p i , p i ∈ P}, where F 0 is the square root of g 0 , is also a δ-power protected (ε, µ)-net, but with respect to the Euclidean metric. We can deduce of Theorem 6.1 for the Euclidean setting. The main result is given by the following theorem: Theorem H.1. Let P be a point set in Ω. Let g be a uniform Riemannian metric field on Ω (∀x ∈ Ω, g(x) = g 0 ). Let C be the canvas and let e C be the length of its longest edge. If then Del d g (P) = Del g (P). Proof. Consider R = Vor d g0 (P) over C. Let F 0 be a square root of G 0 . The matrix F 0 provides a stretching operator between the Euclidean and the metric spaces. Let P 0 = {F 0 p, p ∈ P} be the transformed point set and C 0 be a canvas (a dense triangulation) of the transformed space. Denote by R 0 the discrete Riemannian Voronoi diagram of P 0 with respect to g E over C 0 . Let e C,0 be the upper bound on the canvas edge length of C 0 provided by Theorem 6.1 such that D 0 is captured by C 0 .
Since P 0 is a δ-power protected net with respect to g E , we can invoke Theorem 6.1 and we must thus have for the canvas C 0 to capture Del E (P 0 ). Let C 0 be the image of C by F 0 . Note that C 0 and C 0 are two different triangulations of the same space. If any edge of C 0 is smaller (with respect to the Euclidean metric) than e C,0 , then C 0 satisfies and thus C 0 captures Del E (P 0 ). Recall that given an eigenvector v i of G 0 with corresponding eigenvalue λ i , a unit length in the direction v i in the metric space has length 1/ √ λ i in the Euclidean space. Therefore, if the bound e C of C is smaller than αe C0 , with then every edge of C 0 is smaller than e C,0 . This implies that C 0 captures Del E (P 0 ) and therefore that C captures Del g0 (P).
This settles the case of a uniform metric field.

H.2 Arbitrary metric field
We now consider an arbitrary metric field g over the domain R n . The key to proving the equality of the discrete Riemannian Delaunay complex and the Riemannian Delaunay complex in this setting is to locally approximate the arbitrary metric field with a uniform metric field, a configuration that we have dealt with in the previous section. We shall always compare the metric field g in a neighborhood U with the uniform metric field g = g(p 0 ) where p 0 ∈ U . Because g and the Euclidean metric field differ by a linear transformation, we can simplify matters and assume that g is the Euclidean metric field. The main argument of the proof will be once again that a power protected has stable and separated Voronoi vertices.
We recall the main result of this section, Theorem 5.1.

Inria
Theorem. Let g be an arbitrary metric field on Ω. Assume that P is a δ-power protected (ε, µ)−net in Ω with respect to g. Denote by C the canvas, and e C the length of its longest edge. If where e C,p is given by Lemma H.2, and if ε is sufficiently small and δ is sufficiently large (both values will be detailed in the proof), then Del d g (P) = Del g (P).
We prove Theorem 5.1 by computing for each point p ∈ P the maximal edge length of the canvas such that the Voronoi cell V g (p) is captured correctly. Conditions on ε and δ shall emerge from the intermediary results on the stability of power protected nets.
Lemma H.2. Let ψ 0 ≥ 1 be a bound on the metric distortion and g be a Riemannian metric field on U . Let U be an open neighborhood of Ω = R n that is included in a ball B g (p 0 , r 0 ), with p 0 ∈ U and r 0 ∈ R + such that ∀p ∈ B(p 0 , r 0 ), ψ(g(p), g E (p)) ≤ ψ 0 . Let P U be a point set in U and let p 0 ∈ P U . Suppose that P U is a δ-power protected (ε, µ)−net of with respect to g. Let V g (p 0 ) be the Voronoi cell of p 0 in Vor d g (P U ). If with {λ i } the eigenvalues of g 0 and 0 that is made explicit in the proof, and if ε is sufficiently small and δ is sufficiently large (both values will also be detailed in the proof), then Del d g (P) = Del g (P).

H.2.1 Approach
The many intermediary results needed to prove Theorem 5.1 are presented in Appendices C, D and E. We refer to them at appropriate times.
We use the fact that a Riemannian Voronoi cell V g (p 0 ) can be encompassed into two Euclidean Voronoi cells DV +η E (p 0 ) and EV −η E (p 0 ) that are scaled up and down versions of V E (p 0 ) (Lemma E.12). Specifically, EV −η E (p 0 ) and DV +η E (p 0 ) are defined by where H ω (p 0 , p i ) is the half-space containing p 0 and delimited by the bisector BS(p 0 , p i ) translated away from p 0 by ω 0 . The constant η is the thickness of this encompassing and depends on the bound on the distortion ψ 0 in the neighborhood, and on the sampling and separation parameters ε and µ. We have that η goes to 0 as ψ 0 goes to 1. The (local) stability of the power protected nets assumption is proved here again (Lemmas E.2 and E. 16). From this observation, we can deduce that the Voronoi vertices of the Euclidean Voronoi cell V E (p 0 ) are separated, and thus that the Voronoi vertices of EV E (p 0 ) are separated. A bound on the maximal length of a canvas edge can then be computed such that EV E (p 0 ) is captured and thus V g (p 0 ) is captured.

Difficulties
The difficulty almost entirely comes from proving the stability of the assumption of power protection under metric perturbation in any dimension, that is proving that if we assume that P is δ-power protected with respect to the arbitrary metric field g, then, in a small neighborhood around p 0 , the point set P is δ 0 -power protected with respect to g 0 = g(p). Assuming a power protected net does give us some bounds, but creates a tricky circular dependency as the coefficient δ appears in the dihedral angles (see Lemma D.2). We remedy this issue by proving that Euclidean dihedral angles are bounded assuming power protection with respect to the arbitrary metric field, with Lemmas D.5 and D.6.

H.2.2 Proof of Lemma H.2
Lemma E.5 gives us that V g (p 0 ) lies in DV +η E (p 0 ) and contains EV −η E (p 0 ). Since V g (p 0 ) contains EV −η E (p 0 ), if e C is small enough such that EV −η E (p 0 ) is captured, then V g (p 0 ) is also captured. Proving that EV −η E (p 0 ) is captured is done similarly to the Euclidean setting. While we do not explicitly have the power protected net property for the relaxed Voronoi cells (and specifically, EV −η E (p 0 )), we can still extract the critical property that the Voronoi vertices are separated, as shown by the next lemma.
Lemma H.3. Assume U , g, and ψ 0 as in Lemma H.2. Assume that the point set P U is a δ 0 -power protected (ε, µ)-net with respect to the Riemannian metric field g. Then the Voronoi vertices of EV −η E (p 0 ) are separated. Proof. By Lemmas E.2 and E.16, we have local stability of the power protection and net properties. Hence, P is δ 0 -power protected (ε 0 , µ 0 )-net with respect to g E in U .
Let L 0 = δ 2 0 /4ε 0 be the separation bound induced by the δ 0 -power protection property of P U (see Lemma C.4). Let l be the distance between any two adjacent Voronoi vertices of EV −η E (p 0 ). We know by Lemma E.14 that the parallelotopic region around a Voronoi vertex is included in a ball centered on the Voronoi vertex and of radius χ. The protection parameter ι is given by δ = ιε. We have that where ϕ represents the dihedral angle and s 0 = 1 2 . For the stability regions not to intersect, we require l to be positive. This can be ensured by enforcing that the lower bound is positive: Recall that ε 0 = ψ 0 ε, µ 0 = λε/ψ 0 and ρ 0 = 2ψ 0 ε. Using these notations, we see that l > 0 if

Inria
This condition is easy to satisfy when ψ 0 goes to 1 because the right hand side of Inequality (12) is proportional to (ψ 2 0 − 1)ε 2 . The intermediary results that we use already impose some conditions on δ and ε, and we thus would like to give the condition in Equation 12 in terms of δ, so that it may be compared with Inequality (7). In Lemma E.16, we have seen that which is equivalent to This bound is again proportional to (ψ 2 0 − 1)ε and is very similar to the bound given by Inequality (7), made explicit in Inequality (8), but Inequality (13) provides the tougher bound due to the (ε + χ) coefficient.
We can now provide an upper bound on the length of any canvas edge so that it captures EV −η E (p 0 ) and prove Lemma H.2. From Theorem 6.1, we have that if the canvas edge length is bounded as: e C < min{µ 0 /16, δ 2 0 /64ε 0 }, then V E (p 0 ) is captured as P is a δ 0 -power protected (ε 0 , µ 0 )-net with respect to the Euclidean metric field. As we want to capture EV E (p 0 ), we cannot directly use this result. We have nevertheless obtained the separation between the Voronoi vertices of the eroded Voronoi cell (Equation 11). It is then straightforward to modify the result of Theorem 6.6 by using the separation bound provided in Lemma H.3 instead of the one provided by Lemma C.4. We thus choose e 0 C,p0 = min Remark H.4. We here ignore the consequences of Appendix G as it only complicates formulas without changing the logical steps.
We should not forget that we have assumed that g 0 is the Euclidean metric field, which is generally not the case, we must in fact proceed like for the case of a uniform metric field (Theorem H.1): with {λ i } the eigenvalues of g 0 . Therefore, if the site set satisfies the previous conditions on ε and δ and the canvas is enough for all of its edges to have a length smaller than e C,i , then EV −η E (p 0 ) is captured, and thus V g (p 0 ) is captured, which proves Lemma H.2.
Denote by V (p 0 ) the Voronoi cell of p 0 with the approximate metric. We can then incorporate Lemma H.5 to obtain a result similar to Lemma E.12.
Lemma H.6. Let U , ψ 0 , g, g 0 and P U be defined as in Theorem 5.1. Then we can find η 0 and ω 0 such that These inclusions are illustrated in Figure 15.
∂V +ω g (p 0 ) The subsequent lemmas and proofs are similar to what was done in the case of exact geodesic computations and we do not explicit them.