Delaunay triangulations on orientable surfaces of low genus

Earlier work on Delaunay triangulation of point sets on the 2D ﬂat torus, which is locally isometric to the Euclidean plane, was based on lifting the point set to a locally isometric 9-sheeted covering space of the torus [28, 21, 13, 12]. Under mild conditions the Delaunay triangulation of the lifted point set, consisting of 9 copies of the input set, projects to the Delaunay triangulation of the input set. We improve and generalize this work. First we present a new construction based on an 8-sheeted covering space, which shows that eight copies suﬃce for the standard ﬂat torus. Then we generalize this construction to the context of compact orientable surfaces of higher genus, which are locally isometric to the hyperbolic plane. We investigate more thoroughly the Bolza surface, homeomorphic to a sphere with two handles, both because it is the hyperbolic surface with lowest genus, and because triangulations on the Bolza surface have applications in various ﬁelds such as neuromathematics and cosmological models [32, 4, 16]. While the general properties (existence results of appropriate covering spaces) show similarities with the results for the ﬂat case, explicit constructions and their proofs are much more complex, even in the case of the apparently simple Bolza surface. One of the main reasons is the fact that two hyperbolic translations do not commute in general. To the best of our knowledge, the results in this paper are the ﬁrst ones of this kind. The interest of our contribution lies not only in the results, but most of all in the construction of covering spaces itself and the study of their properties.


Introduction
Delaunay triangulations of the flat cubic torus [39], i.e., the flat torus for which the fundamental domain is a cube, are used in many fields of science, ranging from the infinitely small, with "nuclear pasta" [33] to the infinitely large, with cosmology [35,36,40,24] (more references can be found in [13,14,15,11]).They are needed either because the objects studied are periodic in space, or because they allow for modeling periodic boundary conditions, e.g., in fluid mechanics.Periodic Delaunay triangulations in the hyperbolic plane are also used where Γ is a discrete group of isometries acting on the Euclidean plane E2 (resp.hyperbolic plane H 2 ), with properties made precise later in the paper.
As in [15], we stick to the standard definition of a triangulation [30, Section 9.1, condition (3)] of a topological space Y as a geometric simplicial complex K such that K = σ∈K σ is homeomorphic to Y.A triangulation defined by a point set P is a triangulation whose set of vertices is P. The Delaunay triangulation of a point set P in E 2 (resp.H 2 ) is a triangulation of the convex hull (resp.hyperbolic convex hull) of P that can be defined uniquely even in degenerate cases [20].
The Delaunay triangulation of points on a 2D flat torus was usually obtained from a triangulation of the convex hull of nine periodic copies of the points in the plane E 2 , laid in a 3x3 pattern [28,22] (or 27 copies in a 3x3x3 pattern in the 3D case [21]).A more recent algorithm uses covering spaces instead, and computes the Delaunay triangulation directly in the torus whenever possible [13], thus providing all adjacency relations between simplices.It led to cgal packages [27,12], which have already been used in various fields. 1In that work, the Delaunay triangulation of a surface M is defined as follows: Definition 1.The Delaunay triangulation DT M (P) of M defined by a finite point set P in E 2 (resp.H 2 ) is the projection onto M of the Delaunay triangulation DT E (ΓP) (resp.DT H (ΓP)) of the infinite point set ΓP, if it is a triangulation.
Note that the fact that Γ is discrete guarantees that the triangulation of ΓP in E 2 (resp.H 2 ) is locally finite [15].We elaborate a bit on the projection here and check that it is well-defined.Let the projection be denoted as π : E 2 (resp.H 2 ) → M .By definition, for any given point p ∈ P, all elements of the orbit Γp = {g(p) | g ∈ Γ} project by π to the same point π(p) of M .We can assume that no two points of P lie on the same Γ-orbit: otherwise we can always choose a proper subset P ⊂ P, so that ΓP = ΓP.The Delaunay triangulation DT E (ΓP) (resp.DT H (ΓP)) is a partition of E 2 (resp.H 2 ) by vertices, edges, and triangular faces.Each g ∈ Γ is an isometry, so the empty disk circumscribing a Delaunay triangle (p, q, r), p, q, r ∈ E 2 (resp.H 2 ) is transformed by g into the disk circumscribing triangle (g(p), g(q), g(r)), which is also empty.If the set ΓP contains degeneracies, the Delaunay triangulation of a polygon formed by cocircular points P = (p 0 , p 1 , . . ., p k = p 0 ), k ≥ 4 is not uniquely defined.For any g ∈ Γ, the polygon g(P ) = (g(p 0 ), g(p 1 ), . . ., g(p k ) = g(p 0 )) is also formed by cocircular points.Once a triangulation of P is chosen, the polygon g(P ) can be triangulated by the image by g of the triangulation of P .Then, these two triangulations project under π to the same set of triangles on M .Since π is surjective, the projection π(DT E (ΓP)) (resp.π(DT H (ΓP))) forms a partition of M by vertices, edges, and triangular faces. 2 This discussion allows us to ignore degeneracies in the sequel.
The algorithm follows the classic incremental algorithm in E d [8]: for each new point p, the conflicting simplices (i.e., simplices whose circumscribing balls contain p) are deleted, and the region that they formed is then triangulated with new simplices with apex p.That algorithm relies on the fact that the union of the set of conflicting simplices is always a topological ball, which subsumes that the graph of edges of the triangulation has no cycle 20:3 of length one (1-cycles) or two (2-cycles).However, the graph of edges of the projection of the Delaunay triangulation on a manifold can have non-trivial cycles: for instance, in the extreme case when the input consists of only one point p, then the set Γp is the orbit of p under Γ, so all Delaunay vertices project to the same point, i.e., Delaunay edges project to a small set of 1-cycles on the manifold, all with the same vertex.The existence of such 1or 2-cycles can be checked using a geometric criterion.The algorithm computes in some finite-sheeted covering manifold (i.e., roughly speaking, using "copies" of the input points) satisfying the criterion: a 9-(resp 27-) sheeted covering manifold is sufficient for the 2D (resp.3D) flat torus.When sufficiently many points have been inserted so that the criterion is satisfied on the initial manifold, the computation continues in this manifold, i.e., without using more copies of points.
To the best of our knowledge, there are no known similar results on hyperbolic surfaces.Intuition is challenged there, in particular because the group Γ is non-Abelian in general.
In this paper, we first propose a construction of a family of 2 k sheeted-covering spaces, k > 0, of the flat cubic torus (Section 2).We show how this construction allows us to improve the computation of a Delaunay triangulation of that torus by reducing the number of sheets of the covering spaces.
For any closed hyperbolic surface M , we prove in Section 3 the existence of a finite-sheeted covering space M of M so that the projection of the Delaunay triangulation in H 2 onto M has no 1-or 2-cycle (Proposition 9).We also present a geometric criterion to checks that 1or 2-cycles don't exist (Proposition 8).When such a covering space M can be constructed, then the algorithm previously proposed for closed Euclidean manifolds [15] naturally extends to a set of points P on the hyperbolic surface M .
The Bolza surface M is the simplest possible surface after the 2D Euclidean torus: it is homeomorphic to a double torus.The main contribution of this paper (Section 4) is a construction of some finite-sheeted covering spaces of M, which allows us to concretely exhibit a 4-sheeted covering space in which the projection of a Delaunay triangulation never has a 1-cycle, and a 128-sheeted covering space in which it has no 2-cycle.These results mostly are of theoretical interest for the incremental algorithm above, due to the large number of points that would need to be considered for the computations; we propose a more practical algorithm in the last section.

Some finitely-sheeted covering spaces of the flat torus
Let a, b denote two orthogonal translations of the same length.We denote as G 1 the group < a, b > generated by a and b.Since G 1 is Abelian, any x ∈ G 1 can be uniquely written as x = a α b β , where α, β ∈ Z.The length λ 1 (x) of x can thus be defined as the sum |α| + |β|.
We are going to study some covering spaces of the flat torus where • denotes the Euclidean norm.Any closed square of edge length l = δ(G 1 ) whose edges are parallel to a and b is a fundamental domain of G 1 .
In the sequel, the inverse x −1 of an element x of G 1 will alternatively be denoted by x.All subgroups are normal, since G 1 is Abelian.Let us consider the subgroup G 2 consisting of elements of G 1 of even length.It is easy to see that G 2 = ab, ab .
Proof.Let ϕ : G 1 → Z 2 be the group homomorphism defined by ϕ(x) = λ 1 (x) mod 2. The subgroup G 2 is the kernel of ϕ.According to the First Isomorphism Theorem (see, e.g., [2]) The proof of the following lemma is similar to that of Lemma 1.

Lemma 2. For any
By construction, for any k ≥ 0, the generators of G 2 k are two translations whose vectors a k and b k are orthogonal and of equal length δ(G

For any
Delaunay triangulations via 2 k -sheeted covering spaces, k > 0 As shown in [13] the Delaunay triangulation DT T (P) as defined in introduction (Def. 1) does not exist for each set of points P ∈ E 2 , because π(DT (GP)) can have 1-or 2-cycles.However, there are always some covering spaces of the torus in which the Delaunay triangulation exists and can be computed.Let ∆(S) denote the diameter of the largest disk in E 2 that does not contain any point of a set S in its interior.Note that for any 2 , then for any finite for any finite P ⊇ P, P ∈ E 2 .
), for any p ∈ P. Therefore, by Proposition 5, π 2 k (DT (G 1 P)) is a simplicial complex.Let us now consider the case k = 3. Corollary 6.For any finite set P ⊂ D 1 of n > 1 points, π 8 (DT (G 1 P)) is a simplicial complex.
Proof.The maximum empty disk diameter ∆(G 1 p), for any p ∈ P, is √ 2l, that is equal to 2).We are going to show that ∆(G 1 P) is strictly less than √ 2l.Let us suppose that B is a disk with diameter √ 2l centered at some point c ∈ E 2 and containing no point of G 1 P in its interior (See Figure 2).B contains D 1 (c), the original domain centered at c. Since D 1 (c) is a half-open square, there is only one point of D 1 (c) on the boundary of B. For p and q = p, q ∈ P, there is a representative of G 1 p and a representative of G 1 q in D 1 (c).At least one of these representatives lies in the interior of B. This contradicts that the interior of B is empty.If we add more points, the diameter of the largest empty disk cannot become larger.By Proposition 5, π 8 (DT (G 1 P)) is a simplicial complex.
Algorithm.For a finite set of points P ⊂ D 1 , the algorithm of [13], inspired by the standard incremental algorithm [8], computes DT T (G 1 P) via a 9-sheeted covering space.As soon as the condition of Proposition 4 is fulfilled upon insertion of a point, it switches to computing in T 2 .Assuming that P contains at least two points, we modify this algorithm and use the covering spaces T 2 8 , T 2 4 , and T 2 2 using Proposition 5 and Corollary 6. • The algorithm starts by precomputing π 8 (DT (G 1 {p, q})), for any {p, q} ⊂ P. Then it adds points one by one in the Delaunay triangulation of the 8-sheeted covering space T 2 8 .
• If after insertion of some set of points S ⊆ P, the maximum empty disk diameter becomes less than l, then π 4 (DT (G 1 S )) is guaranteed to be a triangulation for any finite S ⊇ S. So, we can discard all redundant periodic copies of simplices of π 8 (DT (G 1 S)) and switch to incrementally computing π 4 (DT (G 1 S )) in the 4-sheeted covering space T 2 4 for S ⊆ S ⊆ P.
• Similarly, if after some more insertions the maximum empty disk diameter becomes less than √ 2 2 l, then we switch to incrementally computing triangulations in the 2-sheeted covering space T 2  2 .• If it becomes less than 1  2 l, then we switch to computing π(DT (G 1 P)) in T 2 .For any S ⊂ P, π 2 k (DT (G 1 S)) contains 2 k periodic copies of each point of S. Hence, using covering spaces whose number of sheets is as small as possible improves the efficiency of the algorithm, even though the asymptotic complexity is unchanged: it stays randomized worst-case time and space optimal [13].

3
Delaunay triangulations on hyperbolic closed surfaces

Background
As the Poincaré disk is a conformal model of the hyperbolic plane (i.e., it preserves angles), it is widely used in applications [25,32,31,4,16].In this model [10]   is denoted as Isom(H 2 ), and the group of orientation-preserving isometries is denoted as There are three types of orientation-preserving isometries of H 2 ; in this paper we only consider hyperbolic isometries, also called translations.A translation g fixes no point of B and fixes two points on H ∞ .The axis X g of g is the hyperbolic line ending at the two fixed points (Fig. 3(right)).Any point is translated by g along a curve at constant distance to the axis.All such curves have the fixed points of g as infinite points.The distance by which g moves all points of X g is the same, it is called the translation length of g, l(g) = 2 • acosh 1  2 trace(m g ) .If p ∈ X g , then the distance by which g moves p is strictly greater than l(g).It is worth noticing that -two translations do not commute in general.-A Fuchsian group is a discrete subgroup of Isom + (H 2 ) [26].A hyperbolic surface is a connected 2-dimensional manifold, such that every point has a neighborhood isometric to a disk of H 2 .In this paper all hyperbolic surfaces are closed surfaces, i.e., they are compact.
-A closed orientable hyperbolic surface M is of the form H 2 /Γ, --where Γ is a Fuchsian group only containing hyperbolic translations (and the identity) -(See [30,37,38,41].)The projection map π : H 2 → M = H 2 /Γ is a local isometry and a covering projection.The Dirichlet region D p (Γ) for Γ centered at p is defined as the closure Cl(V p (Γp)) of the open cell V p (Γp) of p in the Voronoi diagram VD H (Γp) of Γp in H 2 .From compactness of M follows that D p (Γ) is a compact hyperbolic convex polygon with finitely many sides.The area of M is finite: area(M ) = 4π(genus(M ) − 1).Each Dirichlet region D p (Γ) is a fundamental domain for the action of Γ on H 2 , i.e., (i) g∈Γ g(D p (Γ)) = H 2 , and The systole sys(M ) is the least length of a non-contractible loop on M .In the sequel, it will play a role similar to the role played by the edge length of the fundamental domain above.There are upper bounds depending on the area and the genus of M [23,9,17].The shortest closed curves on M are simple, they are closed geodesics and they are in one-to-one correspondence with the conjugacy classes of elements of Γ.The length of the closed geodesic γ corresponding to a conjugacy class [h], h ∈ Γ, is the translation length l(h) [30].The axis X h projects to γ by π; it can be thought of as "winding" X h on γ infinitely many times.Any gX h , where g ∈ Γ, projects onto the same closed geodesic γ by π; the hyperbolic line gX h is the axis of the translation ghg −1 .Note that the conjugacy classes [h] and [h 2 ] are different; the closed geodesic corresponding to [h 2 ] traverses twice the closed geodesic corresponding to [h].Also, the closed geodesics corresponding to [h] and [h −1 ] are different.

Projecting DT H (ΓP) on a closed hyperbolic surface
A triangulation defined by a point set P in H 2 is a Delaunay triangulation iff the circumscribing disk of each triangle does not intersect H ∞ and is empty, i.e., it does not contain any point of P in its interior.
Let M = H 2 /Γ be a closed hyperbolic surface, with projection map π.The images by π of two vertices of a triangle σ of DT H (ΓP) may coincide; in this case at least one of the edges of π(σ) is a 1-cycle.Any 1-cycle in the graph of edges of π(DT H (ΓP)) is freely homotopic to a simple closed geodesic.A 1-cycle e with a vertex v is a geodesic loop.Note that e is not necessarily a closed geodesic: it may form a non-flat angle at vertex v.We refer the reader to the preprint [7] for the more technical characterization of 2-cycles.
Let us consider disks that do not contain any element of the Γ-orbit of a given point p.The Voronoi diagram of Γp is the partition of H 2 by the images by all elements of Γ of Dirichlet regions D p (Γ). Largest empty disks are centered at images of the vertices of D p (Γ) that are furthest to p; δ p denotes their diameter.We define where O denotes the origin.δ M is completely determined by the group Γ.The following inequality holds: (proof in preprint [7]) For a set of points P, we denote as δ P the largest diameter of disks that do not contain any point of P. Clearly, for p ∈ P, δ P ≤ δ p ≤ δ M .Proposition 8.If sys(M ) > δ P (resp.sys(M ) > 2δ P ) then the graph of edges of π(DT H (ΓP)) has no 1-cycle (resp.no 2-cycle).
(proof in [7]) Let us now focus on the specific case of the Bolza surface.

The Bolza surface
Consider a regular octagon in H 2 centered at the origin O with angle π/4 at each vertex, and the four hyperbolic translations a, b, c, d that pairwise identify opposite sides of the octagon (Fig. 4).They generate a Fuchsian group G for which the octagon is a fundamental domain.The inverse of an element g ∈ G is denoted as g.The group G can be written as a finitely presented group  The Bolza surface is M = H 2 /G.It is a closed orientable hyperbolic surface homeomorphic to a torus having two handles, i.e., a double torus.

S o C G 2 0 1 6 20:8 Delaunay triangulations on low genus surfaces
The group G shows interesting properties in relation to some triangle groups.For integers p, q, r ≥ 2, 1 p + 1 q + 1 r < 1, the triangle group T (p, q, r) is defined as the group generated by reflections in the edges of a triangle ∆(p, q, r) whose angles are π p , π q , π r [30].The set of triangles {g(∆(p, q, r)), g ∈ T (p, q, r)} partitions the plane H 2 .G is a normal subgroup of index 16 without fixed points of T (2,8,8).G is a normal subgroup of index 96 without fixed points of T (2,3,8).These properties of G are illustrated geometrically on Fig. Proposition 10.If two points p and q of H 2 lie in the same orbit under T (2, 3, 8), then the Dirichlet regions D p (G) and D q (G) are isometric.
As will be seen in the sequel, this is not true otherwise, which is one of the striking differences with the Euclidean case.
Proof.q = r(p) for r ∈ T (2, 3, 8).Since G is normal in T (2, 3, 8), rGr −1 = G.Then the orbit Gq = rGr −1 q = rGp is the image of the orbit Gp under r.VD H (Gq) is the image of VD H (Gp) by r.The result follows since D q (G) = Cl(V q (Gq)).
Notation.For simplicity, in the sequel we denote as g the image g(O) of O by g ∈ G.We also label each side of D p (G) by the translation of G that maps D p (G) to the Dirichlet region adjacent through this side.Let ρ ∈ T (2, 3, 8) be the rotation around the origin O by π/4, and let A = (a, b, c, d, a, b, c, d) be the alphabet formed by the generators of G and their inverses, considered as an ordered sequence (Fig. 4).Conjugation by ρ acts as a permutation on A.
The Dirichlet regions for G have been studied in [29].The rest of this section presents a more intuitive description.
For x ∈ A, we denote as v x the vertex of D O (G) that is equidistant to O, x, and ρ −1 xρ.The eight elements of A are the elements of the orbit GO that are closest to O.In the same way, the eight elements xz, z ∈ A are closest to x.Any z ∈ A can be written as ρ k xρ −k , where 0 ≤ k < 8, so the above xz, z ∈ A can be written as xz = x • ρ k xρ −k (O).A careful but straightforward observation of the Dirichlet regions of O, x, x • ρxρ −1 , etc., shows: Observation 11.Each vertex v x is incident to eight isometric Dirichlet regions.See Fig. 6, where Without loss of generality, we now examine more closely the vertex v of D O (G) that is incident to D a (G) and D d (G) (see Fig. 4).The eight Dirichlet regions incident to v are centered at the images of O by the following translations, in clockwise order: F = ( , a, ab, abc, abcd = dcba, dcb, dc, d) (Fig. 7).The eight Dirichlet regions For a point p ∈ H 2 , we define Notation C can be remembered as standing for "crown".For conveniency, we will sometimes use the same notation C(p) for the set of points {g(p) | g ∈ C(p)} (the context will disambiguate).It is easy to check that for q = t(p), t ∈ T (2, 3, 8) we have Let p be a point in triangle SOQ shown in Fig. 8. Let us consider the Dirichlet regions that are adjacent to D p (G).  9 shows the Dirichlet regions of three different points q, obtained by computing the Delaunay triangulation in H 2 of sufficiently many points of Gq, with the software described in [6].Proposition 13.The largest empty disk diameter satisfies 4.89 < δ M < 6.62.(proof in [7]) There exist closed geodesics that are < δ M [3] (see preprint [7]).

Covering spaces of the Bolza surface
The function LowIndexSubgroupsFpGroup( Gamma, H, k, excluded ) of gap [1,34] provides in principle a way to compute the two covering spaces of Prop.9.However, its huge memory consumption does not allow us to use it for the Bolza surface.We propose a construction of coverings spaces of M in which, for any set P, the projection π(DT H (ΓP)) does not contain any 1-or 2-cycle.In a similar spirit as for the flat torus, here we look for subgroups G of G not containing elements corresponding to short loops, and whose associated Dirichlet regions keep some symmetry.
G is a normal subgroup of the triangle groups T (2, 3, 8), T (2, 4, 8), T (2,8,8).They contain the rotation ρ defined earlier, and the reflection τ with respect to the axis X a of translation a.Obviously, if G is a normal subgroup of one of these three triangle groups, then G = ρ Gρ −1 and G = τ Gτ −1 (which holds for G = G).
Let us introduce the subset N O ( G) ⊂ G of elements g such that the Dirichlet regions D O ( G) and D g ( G) are adjacent (i.e., they share a common edge).
, and since G = ρGρ −1 , so ρhρ −1 is in G. Then the translation ρhρ −1 • g is in G, and it fixes O. Since G has no fixed point, it must be the identity, and ρhρ −1 = g.This is true for each The same proof works for τ : the only property of ρ that we used is that ρ(O) = O, which is also true for τ .
It follows as a corollary that if G is normal in a triangle group, then D O ( G) is invariant under ρ and τ .The three triangle groups T (2, m, 8) mentioned above can be generated by ρ, τ and some other rotation m .Before explicitly constructing subgroups of G, let us derive a condition that will be used later (Prop.so that the geodesic segment with endpoints O and g is an edge of DT H (GO).
In the sequel, we use the standard notation [g, h] for the commutator of two isometries g and h, i.e., [g, h] = ghg −1 h −1 .
We first define the subgroup G 2 of G generated by words of even size on A (Fig. 11).G 2 is a subgroup of index 2 in G, hence it is normal in G.
Proof.It is sufficient to prove that any word of size 2 on A can be generated by this subset.Consider yz, y, z ∈ A. As mentioned in Section 4, ρ acts as a permutation on A, i.e., y = ρ k1 xρ −k1 and z = ρ k2 xρ −k2 for some The subgroup G 2 can be written in a finitely presented form: (proof in preprint [7]) The generators of G 2 all have the same translation length l(a) = l(b) = l(c) = l(d).As a result, sys(M 2 ) = sys(M).However, there are fewer shortest loops on M than on M, since the elements of A are not in G 2 .
As for the flat case, we can continue and define inductively the group G 2 k consisting of elements of even size in the alphabet of generators of 6).

S o C G
Proof.Let Π be the polygon with vertices ( Using the fact that N O (G 4 ) generates G 4 , we get: In other words, the generators of G 4 are ac, abcd, db, and their conjugates by ρ k , 1 < k < 8 (Fig. 12).In a similar way as done for G 2 , we can show that G 4 is the only subgroup of index 2 of G 2 for which the Dirichlet region centered at O is invariant by ρ.Thus sys(H 2 /G 4 ) > l(a).The second shortest length of a non-contractible loop corresponds to geodesics of length l(abcd) [3] (see preprint [7]).Since abcd ∈ G 4 (violet point in Fig. 12), the result follows.
The systole of M 4 is greater than sys(M), but still smaller than δ M .However, as a consequence of Proposition 15: Proposition 22.For any P, the projection of DT H (GP) on M 4 has no 1-cycle.
Due to lack of space, we only mention the following result (see preprint [7]): Proposition 23.The index of a subgroup G for which sys(H 2 / G) > 2δ M is greater than 32.There is a unique subgroup G 128 of G of index 128 that is normal in T (2,3,8).For any set of points P, the graph of edges of the projection of DT H (GP) to H 2 /G 128 has no 2-cycle.

Triangulating the Bolza surface in practice
As mentioned in introduction, the incremental algorithm for computing a Delaunay triangulation can be extended to M, using a covering space of M in which the graph of its edges has no 2-cycles, similarly to what was done for Euclidean manifolds [15].However, computing in H 2 /G 128 (Proposition 23) would not be efficient, since we would have to consider 128 images of P.
In this section, we propose a heuristic method instead, inspired from [11,12].It consists in first initializing the triangulation with a well chosen set Q of 14 points for which sys(M) > 2δ Q .Inserting new points in the triangulation cannot increase the size of empty disks, i.e., δ P ∪Q ≤ δ Q for any set P , so by Proposition 8, π(DT H (G(P ∪ Q))) has no 2-cycle.Thus, the incremental algorithm can be used to compute the Delaunay triangulation of P ∪ Q directly in M. In a lot of applications, keeping these 14 extra vertices in the triangulation would not be an issue.However, there may be cases when this should be avoided.If P contains sufficiently many points, and if they are reasonably well distributed, which is likely to be the case in applications, 3 then the final triangulation π(DT H (GP)) would have no 1-or 2-cycles.In such cases, all points of Q can be removed from π(DT H (G(P ∪ Q))) (basically, as in a 2D Euclidean triangulation [19] -practical details are omitted).However, for badly distributed point sets, such short cycles will appear during the process, which can be checked after each removal using the geometric criterion; in such cases, points of Q cannot be removed further.
The rest of this section is devoted to presenting the construction of Q.
The group G is a normal subgroup without fixed point of index 48 in the triangle group T (2,4,8).The octagon D O (G) can be partitioned into 48 triangles that are fundamental domains for T (2, 4, 8) (Fig. 13  Let δ Q denote the largest empty disk diameter for GQ, and d 1 (resp.d 2 ) denote the diameters of the disks circumscribing ∆ 1 (resp.∆ 2 ).There is a triangle in DT H (GQ) whose circumscribing disk has diameter δ Q , and since any triangle in DT H (GQ) is the image of ∆ 1 or ∆ 2 under some isometry of T + (2, 4, 8), either d 1 or d 2 is equal to δ Q .Both d 1 and d 2 are less than 1  2 sys(M), which actually is the radius of the circle inscribed in D O (G).So δ Q < 1  2 sys(M).Fig. 15 shows unique representatives in H 2 of the triangles of π(DT H (GQ)).The point set Q contains 9 points in the interior of D O (origin and red points), 4 points on the boundary of D O (green points) and 1 point at its vertex (blue point).In total Q consists of 14 points.

Open problems
Further work directions are twofold.On a theoretical side, we are planning to reduce the gap between the lower bound 32 and the upper bound 128 in Proposition 23; we conjecture that 32 sheets are always sufficient to avoid 2-cycles.We will also investigate generalizations of the results given here to closed hyperbolic surfaces of higher genus.On a practical side, 3 this was observed and analyzed in the 3D Euclidean case [11,13] S o C G 2 0 1 6 low genus surfaces ker ϕ is a normal subgroup of G 1 , and

Figure 3 (
Figure 3 (left) Lines and circles.(right) Translation g of axis Xg.

Section 4 .
1 shows a covering of the Bolza surface in which the graph π(DT H (ΓP)) has a 2-cycle formed by edges of length δ P , i.e., sys(M ) = 2δ P , so the criterion is sharp.Let now Γ be a subgroup of finite index k in Γ.The quotient space M = H 2 / Γ is a closed hyperbolic surface with projection map π : H 2 → M .The surface M is a k-sheeted covering space of M with covering map and local isometry f : M → M such that π = f • π.Note that Area( M ) = k • Area(M ), genus( M ) = k(genus(M ) − 1) + 1, sys( M ) ≥ sys(M ), and δ ‹ M = δ M .Proposition 9.

Figure 4
Figure 4 Translations a, b, c, d and their inverses.
i.e., the quotient of the group a, b, c, d generated by a, b, c, d, by the normal closure (i.e., the smallest normal subgroup) of the element abcdabcd in a, b, c, d [18, Chapter 5.5].

Figure 8
Figure 8 SOQ.-If p is in the interior of SOQ, the Dirichlet region D p (G) is a polygon with 18 sides.Each vertex of D p (G) has degree 3. The sides of D p (G) have the following labels in counterclockwise order: a, ad, bcd, b, ba, c, d, a, b, cd, c, cba, da, d, dc, dcb, abc, ab.-If p lies in the interior of edge OQ, the sides dc and abc degenerate to a vertex of D p (G).On SO, the sides ab and dcb degenerate to a vertex, and on QS b and c degenerate to a vertex.Such a new vertex has degree 4. The Dirichlet region D p (G) is a polygon with 14 sides.-If p ∈ {O, Q, S}, D p (G) degenerates to an octagon (Proposition 10).Its adjacent octagons are centered at some points in C(p).In particular: D p (G) shares a common vertex with D abcd(p) (G) iff p = O.Fig.9shows the Dirichlet regions of three different points q, obtained by computing the Delaunay triangulation in H 2 of sufficiently many points of Gq, with the software described in[6].

Figure 9
Figure 9Dirichlet region for G centered at different points q (rightmost picture: q = O).

Observation 14 .
The Dirichlet region D O ( G) for G centered at O is invariant by ρ (resp.τ ) if and only if G = ρ Gρ −1 (resp.G = τ Gτ −1 ).Proof.If G = ρ Gρ −1 , then by definition of ρ, ρ GO = ρ Gρ −1 O = GO, and the Dirichlet region D O ( G), as a cell in VD H ( GO), is invariant by ρ (and similarly for τ ).Now assume that D O ( G) is invariant by ρ.Then in particular N O ( G)(O) = ρN O ( G)(O), and we can rewrite

Figure 11
Figure 11 DO(G2).andτ .Indeed, conjugation by ρ or τ leaves A invariant, so it preserves the size of a word in G.By Observation 14, G 2 = ρG 2 ρ −1 = τ G 2 τ −1 .We can actually show: Proposition 18. G 2 is the only subgroup of index 2 of G for which the Dirichlet region centered at O is invariant by ρ.(proof in preprint[7])

Figure 15
Figure 15The set Q.
We call it the original domain of G 2 k and denote it as D 2 k .See Figure 1; for simplicity, x denotes both an element x of G 1 and the image xO of the origin O by x.The closure of the original domain D 2 k is a fundamental domain of G 2 k ; it is also the cell of O in the Voronoi diagram of its orbit G 2 k O.

0 1 6 20:6 Delaunay triangulations on low genus surfaces pencil
[5, Chapter 19]the hyperbolic plane H 2 is represented as the open unit disk B of the Euclidean plane E 2 .The points on the boundary of B are the points at infinity; their set is denoted as H of circles that it generates with H ∞ (Fig.3(left)).The group of isometries of H 2 ∞ = {p ∈ E 2 : p = 1}.Hyperbolic lines, or geodesics, are represented either as arcs of Euclidean circles orthogonal to H ∞ or as diameters of B; a hyperbolic circle is represented as a Euclidean circle contained in B, and its hyperbolic center is the limit point of theS o C G 2