$𝜺$-Net Algorithm Implementation on Hyperbolic Surfaces

We refer the reader to textbooks [9, 2, 7] for more details on topology and hyperbolic surfaces. In this paper, we use the term hyperbolic surface for a closed (connected, compact, and without boundary) oriented hyperbolic surface. Such a surface $S$ can be seen as the quotient of the hyperbolic plane ${ℍ}^2$ under the action of a discrete subgroup $\mathrm{\Gamma }$ of the group of orientation-preserving isometries of ${ℍ}^2$. The surface $S$ is locally isometric to ${ℍ}^2$ and thus has constant curvature -1. Note that Hilbert’s theorem prevents to isometrically embed ${ℍ}^2$ or any hyperbolic surface in ${ℝ}^3$, thus all our illustrations will be drawn in the Poincaré disk model. The computation of distances in ${ℍ}^2$ will be detailed in Section 6.2. We denote as $g$ the genus of $S$, that is its number of handles, and $\sigma$ its systole which is the length of the shortest non-contractible closed geodesic.

The action of $\mathrm{\Gamma }$ on ${ℍ}^2$ induces a projection $\rho :{ℍ}^2\mapsto {ℍ}^2/\mathrm{\Gamma }=S$. For $x\in S$, an element $\tilde{x}\in \mathrm{\Gamma }x:={\rho }^{-1}(x)\subset {ℍ}^2$ is called a lift of $x$. Throughout this paper, objects written with a tilde $\tilde{\cdot }$ are in ${ℍ}^2$ while objects on $S$ are written without.

A fundamental domain for $S$ is a connected subset of ${ℍ}^2$ containing exactly one lift of every point of $S$, except on its boundary. The Dirichlet domain ${𝒟}_{\tilde{b}}$ of a point $\tilde{b}\in {ℍ}^2$ is defined as the closed Voronoi cell of $\tilde{b}$ in the Voronoi diagram of the point set $\mathrm{\Gamma }\tilde{b}$. It is a fundamental polygon for $S$: the interiors of two copies of the domain are disjoint, $\mathrm{\Gamma }{𝒟}_{\tilde{b}}={ℍ}^2$, and its boundary is made of geodesic segments called sides. To each side $s$ of ${𝒟}_{\tilde{b}}$, there is a unique element ${\gamma }_s\in \mathrm{\Gamma }$, called side pairing, such that ${\gamma }_s(s)$ is also a side. The side pairings generate $\mathrm{\Gamma }$, so that ${𝒟}_{\tilde{b}}$ and its side pairings determine the metric of the surface.

We work with the Poincaré disk model in which ${ℍ}^2$ is represented by the open unit disk of $ℂ$ (Figure 1). The geodesics are either diameters of the disk, or circular arcs orthogonal to the unit circle. Hyperbolic circles are Euclidean circles but their centers do not coincide in general.

For a triangulation $T$ of $S$, we note $\tilde{T}$ the infinite periodic triangulation of ${ℍ}^2$ whose vertices, edges and faces are all lifts of those of $T$. Conversely, any vertex, edge or face of $\tilde{T}$ projects on $S$ as a vertex, edge or face of $T$. The triangulation $T$ is a Delaunay triangulation if $\tilde{T}$ is a Delaunay triangulation of ${ℍ}^2$, that is, no circumcircle of a triangle of $\tilde{T}$ contains a vertex of $\tilde{T}$ in its interior.

Let us recall definitions [8]. Let $(X,d)$ be a metric space and $\epsilon >0$. A subset $P\subset X$ is an $\epsilon$-covering if $d(x,P)⩽\epsilon$ for all $x\in X$, i.e., the closed balls of radius $\epsilon$ centered at points of $P$ cover $X$. It is an $\epsilon$-packing if $d(p,q)⩾\epsilon$ for all $p\ne q\in P$, that is, the open balls of radius $\epsilon /2$ centered at points of $P$ are pairwise disjoint. If $P$ is both an $\epsilon$-covering and an $\epsilon$-packing, then it is an $\epsilon$-net, also known as an $(\epsilon ,\epsilon )$-Delone set.

When , we say that $S$ is $\epsilon$-thick. The number of points in an $\epsilon$-packing of a hyperbolic surface $S$ is upper-bounded by $16(g-1)\left(1/{\epsilon }^2+1/{\sigma }^2\right)$, or $16(g-1)/{\epsilon }^2$ for an $\epsilon$-thick surface [17].

Let us summarize the algorithm [17], from which our implementation derives. The algorithm is inspired by Shewchuk’s Delaunay refinement [35]. It starts with a Delaunay triangulation of $S$ with a single vertex $b$. In a nutshell, as long as there is a large triangle in the triangulation, its circumcenter is inserted and the triangulation is updated (see Figure 2): the triangle in which the circumcenter lies is first split into three new triangles, then the Delaunay property is retrieved using edge flips [30, 18]. The set of vertices of the final Delaunay triangulation is an $\epsilon$-net of $S$, and all the intermediate sets of vertices are $\epsilon$-packings.

The algorithm also stores the Dirichlet domain ${𝒟}_{\tilde{b}}$ of a lift $\tilde{b}$ of $b$ together with its side pairings, which can be computed from any fundamental domain for $S$ [14]. The data structure stores the lift in ${𝒟}_{\tilde{b}}$ of each vertex of the triangulation (or one of its lifts if it lies on the boundary of ${𝒟}_{\tilde{b}}$), and for each triangle, the information needed to retrieve its (at most three) lifts having at least one vertex in ${𝒟}_{\tilde{b}}$.

The most crucial step of the algorithm is the location of the circumcenter $c_{\mathrm{\Delta }}$ of a triangle $\mathrm{\Delta }$ in a triangulation of $S$. To this aim, a lift of $c_{\mathrm{\Delta }}$ is located in the lifted triangulation in ${ℍ}^2$, as follows. First, the circumcenter ${\tilde{c}}_{\mathrm{\Delta }}$ of a lift $\tilde{\mathrm{\Delta }}$ of $\mathrm{\Delta }$ with at least one vertex in ${𝒟}_{\tilde{b}}$ is computed. The point ${\tilde{c}}_{\mathrm{\Delta }}$ is a lift of $c_{\mathrm{\Delta }}$ that does not necessarily lie in ${𝒟}_{\tilde{b}}$. To locate ${\tilde{c}}_{\mathrm{\Delta }}$, the algorithm first walks in the tiling $\mathrm{\Gamma }{𝒟}_{\tilde{b}}$ to find the element ${\gamma }_c\in \mathrm{\Gamma }$ such that ${\tilde{c}}_{\mathrm{\Delta }}\in {\gamma }_c{𝒟}_{\tilde{b}}$. The lift $\widetilde{c_b}={\gamma }_c^{-1}{\tilde{c}}_{\mathrm{\Delta }}$ of $c_{\mathrm{\Delta }}$ lying in ${𝒟}_{\tilde{b}}$ can then be computed. Second, the triangle in ${𝒟}_{\tilde{b}}$ in which $\widetilde{c_b}$ lies is identified by checking, for each triangle, its (at most three) lifts having at least one vertex in ${𝒟}_{\tilde{b}}$. It can be shown that the walk in the tiling $\mathrm{\Gamma }{𝒟}_{\tilde{b}}$ traverses a bounded number of Dirichlet domains. The complexity of the algorithm is bounded by $O\left(N^2\right)$, where $N$ is the number of points in the $\epsilon$-net.

At each step of the algorithm, a large triangle is considered. In the original algorithm, all triangles of the current triangulation are checked until a large one is found. Consequently, all triangles that are not affected by an insertion will be checked again for the next insertion. This choice has no effect on the theoretical complexity of the algorithm, but it is time-consuming in practice. Instead, we maintain a list of triangles to be processed by the algorithm. We ran experiments to guide our choice towards an effective way of maintaining this list. The details of these experiments can be found in [16, App A]. The results show that any attempt to maintain the list using criteria involving circumradii, like storing only large triangles or ordering them according to their circumradius, is highly inefficient.

More precisely, we maintain a list $ℒ$ of darts representing the triangles to be processed. The only operations on the list are: pop the front dart, and push new darts to the back. Each dart has a mark represented by a Boolean and we maintain the property that each triangle represented in $ℒ$ has exactly one marked dart; A triangle can have several darts in $ℒ$, but only one is marked. The computation of a circumradius is only performed when a marked dart is popped from the front. The list is initialized with one marked dart for each triangle of the input triangulation. Remark that for a reasonable choice of $\epsilon$, all these triangles are large.

When a dart is popped out from $ℒ$, if it is marked, we unmark it and the circumradius of the triangle it represents is computed. If the triangle is large, its circumcenter splits the triangle containing it as detailed in Section 4, and one dart is marked and pushed to the back of $ℒ$ for each of the three new triangles. Darts are added to $ℒ$ without checking the size of the circumradius of the corresponding triangle and without checking if it is a Delaunay triangle. The Delaunay property is restored using the flip algorithm described in [12] and implemented in [13]. After a flip, in each of the two new triangles, we maintain the property that only one dart is marked. Indeed, such a new triangle may have 0, 1 or 2 marked darts (Figure 6). If it has none, then we mark one and push it to the back of $ℒ$; if it has two, then we leave them in $ℒ$ and unmark one; otherwise there is nothing to do.

When a non-marked dart is popped out from $ℒ$, then nothing is done: it means that the triangle it belongs to is represented by another marked dart in $ℒ$.

The algorithm proceeds with the new front dart until $ℒ$ is empty. Note that in practice, since darts representing new triangles are always pushed at the end of $ℒ$, large triangles tend to be close to the front and are processed before triangles with smaller circumradii.

Remark that no element of $ℒ$ is ever removed from the list, except at the front. As explained in Section 4, when a triangle disappears from the triangulation, its darts stay in the data structure, but their ${\beta }_1$ pointers are modified. So, it can be the case that a dart represented a large triangle when it was inserted in $ℒ$, but it represents a triangle whose circumradius is not greater than $\epsilon$ when it comes to the front of $ℒ$.

When a marked dart is popped out from $ℒ$, the circumcenter of a lift of the triangle represented by the dart is first constructed in order to compute its circumradius. The lifted triangle is the one stored in the anchor associated to the dart. The coordinates of a circumcenter of three vertices with rational coordinates are algebraic numbers of degree two [4]. When it is inserted in the triangulation, a circumcenter becomes a vertex. Handling exact circumcenters would thus lead to cascading the algebraic degrees of coordinates, which would result in computations that are known to be impossible to handle in practice.

The coordinates of the circumcenter ${\tilde{c}}_{\mathrm{\Delta }}$ of a triangle $\tilde{\mathrm{\Delta }}$ are represented by the type CGAL::Sqrt_extension, which handles algebraic numbers of degree two. We do not insert the point ${\tilde{c}}_{\mathrm{\Delta }}$ but a rational approximation denoted $\tilde{c}$. We round each coordinate of ${\tilde{c}}_{\mathrm{\Delta }}$ to a double precision number using the to_double() method of cgal and convert it to an exact rational type. More specifically, coordinates of all vertices of our triangulation are represented by the CGAL::Exact_rational number type, which is a wrapper for the arbitrary-precision rational type mpq_t provided by gmp [21]. All the computations on our data structure are performed using this CGAL::Exact_rational number type for points and cross-ratios. Note that the bitsize of the rationals encoding an inserted point is bounded whereas the bitsize of the rationals encoding a cross-ratio grows when it is updated during a flip.

Let us detail the computation of the largeness test for triangles that we use to avoid as much as possible the use of hyperbolic trigonometric functions. In the Poincaré disk, the distance between two points $\tilde{u}$ and $\tilde{v}$ is $d_{{ℍ}^2}(\tilde{u},\tilde{v})=\mathrm{arcosh}(1+\lambda (\tilde{u},\tilde{v}))$, where $\lambda (\tilde{u},\tilde{v})={\displaystyle \frac{2{\Vert \tilde{u}-\tilde{v}\Vert }^2}{(1-{\Vert \tilde{u}\Vert }^2)(1-{\Vert \tilde{v}\Vert }^2)}}$ and $||\cdot ||$ is the Euclidean distance. We compute the minimum $\lambda$ between the approximate circumcenter $\tilde{c}$ and the three vertices of the anchor. If it is greater than a certified upper bound of $\cosh (\epsilon )-1$, then we consider the triangle large and insert $\tilde{c}$ in the triangulation. The certified upper bound on $\cosh (\epsilon )$ is computed with the Boost interval library [19]. Due to approximations for the computation of $\tilde{c}$, this process may miss a large triangle and thus does not enforce the $\epsilon$-covering property. Similarly, even if the point $\tilde{c}$ that is inserted into the triangulation is at distance at most $\epsilon$ from the vertices of its triangle, it may be at distance smaller than $\epsilon$ from another vertex of the triangulation. We will check these properties in the final step of the algorithm (Section 6.4).

Our implementation maintains a Delaunay triangulation of $S$ using the data structure presented in Section 4. Each point to be inserted is located in this triangulation by locating a lift in the lifted triangulation. We do not store a Dirichlet domain as in the original algorithm (Section 2.3), but instead directly walk in the current triangulation itself.

We tested two walk algorithms: the straight walk and the visibility walk (Figure 7) mentioned in Section 3. In our case, the query point $\tilde{c}$ is the approximate circumcenter of the lifted triangle $\tilde{\mathrm{\Delta }}$ given by the anchor of the triangle $\mathrm{\Delta }$ being processed by the algorithm. The walk starts from a vertex $\tilde{p}$ of the triangle $\tilde{\mathrm{\Delta }}$. Lifts of triangles must be constructed along the walk to find the lifted triangle containing $\tilde{c}$.

Both walks are based on orientation tests in ${ℍ}^2$ with vertices of lifted triangles. The combinatorial map encoding the triangulation provides a direct access to the neighbor ${\mathrm{\Delta }}^{\prime }$ of the triangle $\mathrm{\Delta }$ adjacent through one of its edges. However, the lift ${\tilde{\mathrm{\Delta }}}^{\prime }$ of this triangle given by its anchor may not be adjacent to the lift $\tilde{\mathrm{\Delta }}$. The lift of ${\mathrm{\Delta }}^{\prime }$ that is adjacent to $\tilde{\mathrm{\Delta }}$ is computed from the vertices of $\tilde{\mathrm{\Delta }}$ and the cross-ratio of the common edge as explained in Section 4.

Running the epsilon_net method with both walks shows that they perform equivalently (see [16, App B]): the running time of the method remains the same, and they compute a similar number of lifts. The visibility walk does not have to handle the degenerate cases of the straight walk, which may go through a vertex or along an edge; we therefore choose the visibility walk. The bound on the length of the straight walk for the original algorithm, using Dirichlet domains, does not apply to the visibility walk. However we observe in Section 7.2 that the walk has constant length in practice.

As explained in Section 6.2, our algorithm is robust in the sense that the output triangulation is the Delaunay triangulation of a set of points with rational coordinates. On the other hand, the vertices of this triangulation may not be an $\epsilon$-net, due to approximations of circumcenters. To check the $\epsilon$-covering property, it is sufficient to check that there is no large triangle. For every triangle $\mathrm{\Delta }$, we compute the exact circumcenter ${\tilde{c}}_{\mathrm{\Delta }}$ of its anchor, which has coordinates in a degree two algebraic extension, using the number type CGAL::Sqrt_extension. We then compute $\lambda ({\tilde{c}}_{\mathrm{\Delta }},\tilde{v})$ with the same type (see Section 6.2) for a vertex $\tilde{v}$ of the triangle, and check that it is less than a certified lower bound on $\cosh (\epsilon )-1$. To check the $\epsilon$-packing property, it is sufficient to check that all Delaunay edges are longer than $\epsilon$. Since all vertices have rational coordinates, the value $\lambda ({\tilde{v}}_i,{\tilde{v}}_j)$ for two vertices of an edge is rational and it is checked to be larger than the upper bound on $\cosh (\epsilon )-1$ computed with the Boost interval library [19]. When these tests succeed, the vertices of the output triangulation are then certified to be an $\epsilon$-net of the surface. Otherwise, the failure of these tests means that the double precision used for the approximation of circumcenters was not enough.

Our experiments presented in Section 7 show that for surfaces with a long systole, which is the most common case (see Section 7.1), our algorithm is successful. Using a double precision to set the rational coordinates of approximate circumcenters and to compute the bounds on $\cosh (\epsilon )-1$ actually constructs valid $\epsilon$-nets. On the other hand, we also observe in the full version of this article [16, Sec 7.3] that higher precision would be needed to handle a surface with a very small systole.

Our future work will focus on certifying the largeness test for triangles and check the $\epsilon$-packing property after each insertion. If the $\epsilon$-packing check fails, this means that the inserted point was not a good enough approximation of the circumcenter and iteratively refining it (which is possible via the cgal Algebraic Kernel [3]) would eventually recover the $\epsilon$-packing property. Note that iterative refinement of the rational bounds on $\cosh (\epsilon )-1$ must also be computed (which is possible via the Boost library).

To experiment with an algorithm in practice, it is important to sample the space of possible input with some distribution. We discuss this issue below and explain why our random generator of surfaces is likely to produce surfaces with long systole.

For each of the 180 random surfaces, an input for the $\epsilon$-net algorithm is computed as a single vertex Delaunay triangulation of the surface, as explained in Section 5. The epsilon_net method is run with 50 values of $\epsilon$ on every input, decreasing from 0.5 to 0.01 with a step of 0.01. Table 1 presents an overview of the results.

We first observe that the number of vertices in the computed $\epsilon$-nets is close to the theoretical upper bound for an $\epsilon$-thick surface, which is $16(g-1)/{\epsilon }^2=16/{\epsilon }^2$ ($g=2$ for all surfaces of the experiments). In proportion, the number of vertices is 54% of the upper bound on average, with a standard deviation of 2%. The minimum is 47% and the maximum is 63% over all tested surfaces and values of $\epsilon$.

Even though, as in the Euclidean case [23], we have no complexity analysis of the visibility walk for the point location, we observe that the walk traverses a (small) constant number of triangles. The walk locates the approximate circumcenter of a triangle in ${ℍ}^2$, starting from this triangle. The average number of computed lifted triangles during each walk tends to decrease when $\epsilon$ becomes smaller, while the average proportion of approximate circumcenters located in their own triangle tends to increase, as shown in Figure 8. On average, 68% of the approximate circumcenters lie in their triangle. Over all surfaces and all tested values of $\epsilon$, the farthest located points were 4 triangles away from the starting triangle of the walk.

We counted the number of flips done to restore the Delaunay property of the triangulation after each split. On average, it decreases when $\epsilon$ becomes smaller, going from 3.41 for $\epsilon =0.5$ to 2.41 for $\epsilon =0.01$, as shown in Figure 9. This shows that the number of flips at each insertion is a small constant in practice. The total number of flips is thus, in practice, linear in the number of points, which is in contrast with the theoretical quadratic bound (see Section 2.3).

In the full version of this article [16, App C], we show an order of magnitude for the running times to give an idea of the practical complexity of our implementation. It appears to be slightly faster than quadratic in $1/\epsilon$.

To visualize an $\epsilon$-net, we draw a fundamental domain of the surface in the Poincaré disk using its Delaunay triangulation. To do so, it suffices to draw a lift of each triangle in a connected way. We cannot just use the anchors since they would generally not lead to a connected domain. So, we access the anchor of a chosen triangle and build lifts of the other triangles starting from this initial anchor as explained in Section 4. As objects appear smaller near the boundary of the Poincaré disk, we center the drawing at the origin to be able to observe the details. To achieve this, we translate the vertices of the initial anchor such that one of them is 0 before computing the rest of the drawing. If we want similar drawings when running the $\epsilon$-net algorithm with different values of $\epsilon$ on a given surface, we always use an anchor having a fixed lift $\tilde{b}$ of $b$, the unique vertex of the input triangulation (Section 5), as the initial anchor; we then apply the translation that translates $\tilde{b}$ to 0.

In Figures 10 and 11, the drawings are computed by the lift method of the cgal package 2D triangulations on hyperbolic surfaces [13], which uses a weight on edges to order the lifts of triangles of the triangulation $T$ of $S$. The weight of an edge $(\tilde{x},\tilde{y})$ is defined as ${|\tilde{x}|}^2+{|\tilde{y}|}^2$ ($|\cdot |$ being the complex modulus). The drawing is then iteratively computed: given $T^{\prime }$ the set of triangles that have been lifted, the next triangle to be lifted is the one in $T\setminus T^{\prime }$ that is incident to the edge of least weight in $T^{\prime }$.

Figure 10 shows the obtained Delaunay triangulation of the surface with a small systole studied in the full version of this article [16, Sec 7.3]. Since the systole is very small, there is a very long collar around the shortest non-contractible closed curve.

The drawings of triangulations shown in Figure 11 naturally look like Dirichlet domains since the weights on the edges ensure that the lifts of triangles entirely contained in ${𝒟}_{\tilde{b}}$, the Dirichlet domain of $\tilde{b}$, are drawn. This becomes clear when $\tilde{b}$ is translated to the origin (Figure 12). Since the input of the $\epsilon$-net algorithm is a Delaunay triangulation with the single vertex $b$, we obtain ${𝒟}_{\tilde{b}}$ by computing the circumcenter of all the lifted triangles incident to $\tilde{b}$ in the input triangulation, before running the $\epsilon$-net algorithm.

An alternative drawing is obtained by ordering the lifts of the triangles around the initial anchor following a Breadth First Search (BFS) algorithm on the adjacency graph. Such a drawing represents a combinatorial Dirichlet domain (see Figure 13). We observe that, when $\epsilon$ decreases, the combinatorial Dirichlet domain fits the Dirichlet domain ${𝒟}_{\tilde{b}}$ better. Formalizing this convergence is an interesting open question.