Computing the Intrinsic Delaunay Triangulation of a Closed Polyhedral Surface

Polyhedral surfaces can also be obtained from meshes, flat triangles in ${ℝ}^3$ glued along their edges. Moreover, every mesh defines a triangulation of its surface. Yet triangulations are more general than meshes: most triangulations cannot be obtained from a mesh. Some recent algorithms advantageously operate on triangulations of polyhedral surfaces without reference to a mesh [25, 31, 17]. In this context the adjective “intrinsic” is sometimes placed before the name “triangulation” to make the distinction with the particular triangulations arising from a mesh. In the mathematical community, a prominent example is that of a translation surface [21, 32, 14, 12], which arises naturally in the study of billiards in rational polygons.

Triangulations are so general that not all of them are suitable for computation, compared to meshes. Prominently, a fundamental problem on polyhedral surfaces is to compute the distance or a report a shortest path between two points. On a mesh, shortest paths can be computed in time polynomial in the number of triangles. For example, an algorithm of Mitchell, Mount, and Papadimitriou [22] propagates waves along the surface, starting from the source, in a discrete manner. See also Chen and Han [4]. On a generic triangulation however (not arising from a mesh), the number of times a shortest path visits a triangle is not bounded by any function of the number of triangles, as noted for example almost 20 years ago by Erickson [8]. Recently, Löffler, Ophelders, Staals, and Silveira [20] coined the term happiness of a triangulation, for the maximum number of times a shortest path visits a triangle. They adapted the single-source shortest paths algorithm of Mitchell, Mount, and Papadimitriou [22] from meshes to triangulations, whose time complexity now depends on the happiness of the triangulation (it is more efficient on triangulations of low happiness).

This raises the problem of replacing any given triangulation by another triangulation of the same surface whose happiness is “low”. Among the many remeshing algorithms [13, 28, 24, 29], only few have been ported to the general context of intrinsic triangulations [26], and the only solution we are aware of that compares to our main result (Theorem 1 below), by Löffler, Ophelders, Staals, and Silveira [20, Section 5], is restricted to particular inputs whose surfaces are all homeomorphic to an annulus; we will use it as a black box. Importantly, the same authors also showed that Delaunay triangulations have bounded happiness [20, Section 4.2].

Delaunay triangulations are classical objects of computational geometry [5, 11, 2], closely related to shortest paths. While mostly known in the plane, they generalize to closed polyhedral surfaces, see for example the depiction of Bobenko and Springborn [3]. To compute a Delaunay triangulation from an arbitrary intrinsic triangulation there are, to our knowledge, only two approaches, and neither compares to our main result (Theorem 1). One approach computes a Voronoi diagram with a suitably adapted multiple-source shortest path algorithm, and then derives from it a Delaunay tessellation, see for example Mount [23], Liu, Chen, and Tang [18], and Liu, Xu, Fan, and He [19]. Another approach starts from an initial triangulation and flips its edges until it reaches a Delaunay triangulation; it was proved to terminate by Indermitte, Liebling, Troyanov, and Clémençon [3, 15].

In order to state our results precisely, it now matters to make the distinction between a triangulation and the data structure representing it, and to allow for more general polygons than triangles. Following Löffler, Ophelders, Staals, and Silveira [20], we call portalgon the collection $T$ of polygons in the Euclidean plane and the partial matching of their sides. We denote by $𝒮(T)$ the associated polyhedral surface. We say that the portalgon $T$ is triangular if all polygons are triangles. The sides of the polygons, once identified, constitute a graph $T^1$ embedded on $𝒮(T)$ (the polygons themselves correspond to the faces of $T^1$): it is this graph $T^1$ that we call triangulation if $T$ is triangular, and we call $T^1$ a tessellation in general.

We consider, on a closed polyhedral surface $S$, the unique Delaunay tessellation $𝒟$ of $S$ whose vertices are exactly the singularities of $S$, with a single very special exception: if $S$ has no singularity, then $S$ is a flat torus and we let $𝒟$ be any of the Delaunay tessellations of $S$ that have exactly one vertex, for one can be mapped to the other via an orientation-preserving isometry of $S$ anyway. In any case, we say that $𝒟$ is the Delaunay tessellation of $S$, in a slight abuse.²²2Given a set $V$ of points of the surface, finite, non-empty, and containing all the singularities, our results easily extend to Delaunay triangulations whose vertex set is $V$, but this is incidental to us. It is “generically” a triangulation, but not always. If not, then triangulating the faces of $𝒟$ along any vertex-to-vertex arcs provides a Delaunay triangulation. The aspect ratio of a triangular portalgon $T$ is the maximum side length of a triangle of $T$ divided by the smallest height of a triangle of $T$ (possibly another triangle). Our main contribution is:

As already mentioned, the only two methods we are aware of for computing a Delaunay tessellation from an arbitrary triangulation are the flip algorithm and the computation of the dual Voronoi diagram. The time complexities of these algorithms are not bounded by any polynomial in $n$ and $\log (r)$.

Applications of Theorem 1 are detailed in the full version of the paper. Briefly, on the portalgon returned by Theorem 1, shortest paths can be computed in $O(n^2{\log }^{O(1)}n)$ time. And Theorem 1 enables to test whether the surfaces of two given portalgons are isometric, simply by computing and comparing the portalgons of the associated Delaunay tessellations.

We analyze our algorithms within the real RAM model of computation described by Erickson, van der Hoog, and Miltzow [10]. It is an extension of the standard integer word RAM, with an additional memory array storing reals, and with additional instructions. On such a machine, we represent each polygon of a portalgon $T$ by the list of its vertices, and each vertex is by its two coordinates, stored in the memory array dedicated to reals. So displacing (translating or rotating) the polygons in the plane provides different representations of $T$. When modifying a portalgon $T$, we actually modify our representation of $T$, using elementary operations that are easily seen to be achievable by a real RAM.

Within this model of computation, our second contribution is a lower bound that backs our main result, Theorem 1, by showing that the polynomial dependency in the logarithm of the aspect ratio is unavoidable:

Altogether Theorem 1 and Theorem 2 show that, within our model of computation, the complexity of computing the Delaunay tessellation from an arbitrary triangulation of a (closed, orientable) polyhedral surface is polynomial in the number of triangles and in the logarithm of the aspect ratio of the input triangulation.

The proof of Theorem 2 is deferred to the full version of the paper. This extended abstract focuses on the proof of Theorem 1, of which we now provide an overview.

We introduce a slight variation of happiness, more suitable to our needs, which we call segment-happiness. To prove Theorem 1, the crux of the matter is to replace the input triangular portalgon by another triangular portalgon of the same surface, whose segment-happiness is “low”. For this purpose, our approach is to first focus on portalgons $T$ whose surface $𝒮(T)$ is flat: the interior of $𝒮(T)$ has no singularity. Note that here we allow $𝒮(T)$ to have boundary, and this boundary may have singularities. The systole of $𝒮(T)$ is the smallest length of a non-contractible geodesic closed curve in $𝒮(T)$. Our key technical result is:

Sections 3–5 are devoted to the proof of Proposition 3. In Section 3 we focus on particular portalgons, whose surfaces are all homeomorphic to an annulus; the definitions and results of this section are used by the algorithm of Proposition 3. In Section 4 we describe the algorithm for Proposition 3. It is a finely tuned combination of elementary operations such as inserting and deleting edges and vertices in graphs. While the algorithm itself is relatively simple, its analysis is more involved, and is sketched in Section 5. In this section, we first provide a combinatorial analysis, and then we prepare for the geometric analysis by introducing a new parameter on the simple geodesic paths $e$ of a flat surface, enclosure, possibly of independent interest. Informally, $e$ is enclosed when a short non-contractible loop can be attached to a point of $e$ not too close to the endpoints of $e$. We then use enclosure to analyze the algorithm from a geometric point of view, proving Proposition 3.

We then extend Proposition 3 from flat surfaces to surfaces having singularities in their interior, essentially by cutting out caps around these singularities. To get a cleaner result, we also replace $2+L/s$ by the aspect ratio of $T$, and we replace segment-happiness by happiness, obtaining:

The proof of Proposition 4 is omitted from this extended abstract. We have not discussed Delaunay tessellations yet. Still, we are almost ready to prove Theorem 1. Indeed, once we have a portalgon of low happiness, we can compute shortest paths on the surface. And, as already mentioned, shortest path algorithms classically extend to construct Voronoi diagrams and then Delaunay tessellations. Formally:

The proof of Proposition 5 is omitted from this extended abstract. We insist that the proof of Proposition 5 is incidental to us, and Proposition 5 is not surprising at all. Our contribution is really the proof of Proposition 4. Theorem 1 is immediate from Proposition 4 and Proposition 5:

A portalgon $T$ is a disjoint collection of oriented polygons in the Euclidean plane, together with a partial matching of the sides of the polygons such that any two matched sides have the same length. It is triangular if all polygons are triangles. See Figure 1. Any subset of the polygons defines a sub-portalgon $T^{\prime }$ of $T$: two sides of polygons are matched in $T^{\prime }$ if and only if they are matched in $T$. In a portalgon $T$, identifying the matched sides, isometrically, and respecting the orientations of the polygons, provides the surface of $T$, denoted $𝒮(T)$; it is a 2-dimensional Riemannian manifold whose metric may have singularities. The sides of the polygons of $T$ correspond to a graph $T^1$ embedded on $𝒮(T)$, the 1-skeleton of $T$.

A polyhedral surface is any Riemannian manifold $S$ (possibly with singularities) isometric to the surface of a portalgon. And when we say that a portalgon $T$ is a portalgon of $S$, we implicitly fix an isometry between $𝒮(T)$ and $S$. A tessellation of $S$ is any 1-skeleton of a portalgon of $S$, it is a triangulation if the portalgon is triangular.

Consider a polyhedral surface $S$, a triangulation $T^1$ of $S$, a vertex $x$ of $T^1$, and the sum $a$ of the angles of faces of $T^1$ around $x$. The point $x$ is a singularity if $x$ lies in the boundary of $S$ and $a\ne \pi$, or if $x$ lies in the interior of $S$ and $a\ne 2\pi$. Every other point of $S$ is flat. This does not depend on any particular triangulation of $S$. A surface $S$ is flat if its interior has no singularity (although its boundary may have singularities). The closed flat surfaces are called flat tori.

The aspect ratio of a triangular portalgon $T$ is the maximum side length of a triangle of $T$ divided by the smallest height of a triangle of $T$ (possibly another triangle). Note that the aspect ratio is always greater than or equal to $2/\sqrt{3}>1$, because the maximum side length of a triangle is always greater than or equal to $2/\sqrt{3}$ times its smallest height.

The systole of a polyhedral surface $S$ is the smallest length of a non-contractible geodesic closed curve in $S$, except in the particular case where every closed curve in $S$ is contractible, in which case the systole is $\mathrm{\infty }$. The important thing is that for every positive real $s$ smaller than the systole of $S$, any non-contractible closed curve in $S$ is longer than $s$.

The happiness of a portalgon $T$ is the maximum number of times a shortest path in $𝒮(T)$ visits the image of a polygon of $T$, maximized over all the shortest paths of $𝒮(T)$ and all the polygons of $T$ (see [20, Section 3]). We introduce a variation, more suitable to our needs. In a polyhedral surface $S$, a segment is a simple geodesic path $e$ whose relative interior is disjoint from any singularity of $S$. The segment-happiness of $e$ in $S$, denoted $h_S(e)$, is the maximum number of intersections between $e$ and a shortest path of $S$, maximized over all the shortest paths of $S$. The segment-happiness of a portalgon $T$ is then the maximum segment-happiness $h_{𝒮(T)}(e)$, maximized over the edges $e$ of its 1-skeleton $T^1$. A priori, the segment-happiness of a portalgon $T$ differs from the happiness of $T$. Indeed a path in $𝒮(T)$ may visit many times a face of $T^1$ without intersecting any edge of $T^1$ more than once, if the face has high degree. However, if $T$ is triangular, then the happiness and the segment-happiness of $T$ do not differ by more than a constant factor.

Informally, our goal is to “improve the geometry” of a triangular portalgon $T$. We will make this precise in Section 5. Roughly, the issue is that, without any condition on $T$, the edges of $T^1$ that lie in the interior of $𝒮(T)$ can be arbitrarily long, so one of them may intersect some shortest path arbitrarily many times by wrapping around the surface, and so the segment-happiness of $T$ can be arbitrarily large. A naive way of shortening an edge is to cut the edge in two at its middle point.

InsertVertices. Given a triangular portalgon $T$, consider every edge $e$ of $T^1$ that lies in the interior of $𝒮(T)$, and insert the middle point of $e$ as a vertex in $T^1$.

Applying InsertVertices to a triangular portalgon $T$ produces a portalgon $T^{\prime }$ whose polygons are usually not triangles. We now consider transforming $T^{\prime }$ into a triangular portalgon. To do that we repeatedly cut the polygons of $T^{\prime }$. We need a definition. In the plane consider a polygon $P$, two vertices $u\ne v$ of $P$, and the rectilinear segment $a$ between $u$ and $v$. If the relative interior of $a$ is included in the interior of $P$ then $a$ is called a vertex-to-vertex arc of $P$. It is easily seen that if $P$ is not a triangle then $P$ has at least one vertex-to-vertex arc. Among the vertex-to-vertex arcs of $P$, the shortest ones are the shortcuts of $P$. We emphasize that we consider the shortest ones among all the vertex-to-vertex arcs, without fixing the endpoints, but the endpoints are chosen among the vertices of $P$. In a portalgon $T$ every polygon $P$ corresponds to a face $F$ of $T^1$, and every shortcut of $P$ corresponds to a path whose relative interior is included in $F$: we say of this path that it is a shortcut of $F$.

InsertEdges. Given a portalgon $T$, as long as there is a face of $T^1$ that is not a triangle, insert a shortcut of this face as an edge in $T^1$.

We shall apply InsertVertices followed by InsertEdges to a triangular portalgon $T$ in order to produce another triangular portalgon $T^{\prime }$, hopefully with a “nicer geometry”. The problem is now that $T^{\prime 1}$ has more vertices than $T^1$. All the other operations of the algorithm are devoted to keeping the number of vertices low.

From now on it is important that every surface considered is flat, there is no singularity in its interior. Given a triangular portalgon $T$, assuming that the surface $𝒮(T)$ is flat, we consider decreasing the number of vertices of $T^1$. To do that we naturally consider deleting some vertices. Not all vertices can be deleted. For example a vertex incident to a loop edge cannot be deleted. Also we will not delete vertices that lie on the boundary of the surface $𝒮(T)$. A vertex of $T^1$ is weak if it lies in the interior of $𝒮(T)$ and is not incident to any loop edge in $T^1$. It is strong otherwise.

DeleteVertices. Given a triangular portalgon $T$ whose surface $𝒮(T)$ is flat, construct a maximal independent set $V$ of weak vertices of $T^1$ that have degree smaller than or equal to six. For every vertex $v\in V$ delete $v$ and its incident edges from $T^1$.

After performing DeleteVertices, the polygons of $T$ are usually not triangles anymore, but this will be solved by applying InsertEdges after each application of DeleteVertices. Observe that in DeleteVertices we delete only vertices of degree smaller than or equal to six. Informally, the reason is that deleting a weak vertex of degree $d\ge 3$ creates a face of degree $d$ around it. We then insert $d-3$ edges in this face when applying InsertEdges. The problem is that only a constant number of edges can be inserted in each face without risking to destroy our improvements on the geometry of the tessellation. This is why we make sure that $d=O(1)$ beforehand. The exact bound on $d$ is not really important (although changing it would change some constants of the algorithm), but it must be at least six so that we can still remove a fraction of the excess vertices this way, at least when most of them are strong. Similar ideas can be found in the literature, see for example Kirkpatrick [16, Lemma 3.2].

The operation DeleteVertices cannot delete strong vertices, and among them the vertices that lie the interior of the surface and are incident to a loop edge. In this section we describe an operation for deleting such vertices.

In order to grasp the intuition, observe, informally, that it is possible that almost all the vertices of $T^1$ lie in the interior of $𝒮(T)$ and are incident to a loop edge. Fortunately, it turns out that in this case there must be a sub-portalgon $X$ of $T$ such that $X$ is a tube and the interior of $𝒮(X)$ contains loop edges of $X^1$. We delete such loop edges by replacing $X$ by a good biface with Proposition 7. There is one subtlety: we must choose $X$ carefully so that we replace any concatenation of tubes by a single biface when possible, in order to delete the loops in-between the tubes, instead of replacing each tube individually. That leads to:

SimplifyTubes. In a triangular portalgon $T$ whose surface $𝒮(T)$ is flat, do the following:

1. 1.

   In $T^1$ build a set $J$ of loop edges that lie in the interior of $𝒮(T)$ and are pairwise disjoint, as follows. There are two cases:

   1. (a)

      If $𝒮(T)$ is homeomorphic to a torus, do the following. Let $J$ contain two disjoint loop edges of $T^1$ if there exist two such edges, otherwise let $J=\mathrm{\varnothing }$.

   2. (b)

      Otherwise do the following. Construct a set $J^{\prime }$ of loop edges by considering every vertex $v$ of $T^1$ that lies in the interior of $𝒮(T)$ and is incident to a loop edge, and by putting one of the loop edges incident to $v$ in $J^{\prime }$. Then build a subset $J\subseteq J^{\prime }$ by removing from $J^{\prime }$ every $e\in J^{\prime }$ satisfying both of the following. First, cutting $𝒮(T)$ along the loops in $J^{\prime }$, and considering the resulting connected components, two such components are adjacent to $e$ (instead of one), say $S_0$ and $S_1$. Second, each one of the two sub-portalgons of $T$ whose surfaces are $S_0$ and $S_1$ is a tube.

2. 2.

   Cut the surface $𝒮(T)$ along the loops in $J$. Each resulting component is the surface of a sub-portalgon $X$ of $T$. If $X$ is a tube replace $X$ by a good biface $B$.

The idea behind step 1b is to remove loops from $J^{\prime }$ so that step 2 replaces a concatenation of tubes by a single good biface when possible, instead of replacing the tubes individually.

We are almost ready to give the algorithm, but there is still one important thing to describe. In step 2 of SimplifyTubes, if the good biface $B$ is thin we will not just replace $X$ by $B$, but we will also make sure to not modify $B$ ever again. In this sense $B$ becomes inactive. Doing so requires a data structure remembering which parts of the portalgon are inactive.

See Figure 3. The data structure maintains a portalgon $R$ together with a partition of the polygons of $R$. Each set $X$ of polygons in the partition defines a sub-portalgon of $R$ which we call region. One region is singularized as the active region $R_A$. The other regions are inactive. Note that the surface of the active region may be disconnected, and that the surfaces of distinct inactive regions may be adjacent.

The data structure will be initialized by setting $R_A=R$, without inactive region. Then the algorithm will apply the routines InsertVertices, InsertEdges, DeleteVertices, and SimplifyTubes to the active region $R_A$, and mark as inactive every thin biface encountered in step 2 of SimplifyTubes. The surface of $R_A$ will diminish over time as more and more regions are marked inactive. This may increase the numbers of connected components and boundary components of $𝒮(R_A)$, ruining our efforts to keep the combinatorial complexity of $R_A$ bounded. To counteract this, we introduce:

Gardening. Every connected component of $𝒮(R_A)$ is the surface of a sub-portalgon $X$ of $R_A$. If $X$ is a tube replace $X$ by a good biface $B$, and mark $B$ as inactive.

We described everything that the algorithm can do to the data structure. This immediately implies three invariants maintained by the algorithm. First, 1) Every polygon of the active region has degree at most six, and 2) Every inactive region is a good biface. For the last invariant we need a definition. Recall that in $R$ if two sides $s$ and $s^{\prime }$ of polygons are matched then $s$ and $s^{\prime }$ correspond to an edge $e$ of $R^1$. If moreover $s$ and $s^{\prime }$ belong to different polygons, and if their respective polygons belong to different regions, we say that $e$ is separating. Then $e$ is a loop, for it is a boundary edge of a biface by 2), and $e$ belongs to the interior of $𝒮(R)$. The third invariant is that 3) The separating loops are pairwise disjoint (no two of them are based at the same vertex of $R^1$).

The algorithm repeatedly applies two parts. The first part “improves the geometry” by applying InsertVertices and then InsertEdges. However this increases the number of vertices. So the second part applies SimplifyTubes, DeleteVertices, and InsertEdges, together with Gardening. The second part can only remove a fraction of the vertices at once, so it is repeated several times. It turns out that 350 repetitions suffice.

Algorithm. Given a triangular portalgon $T$ whose surface $𝒮(T)$ is flat, and $N\ge 1$, do the following. Initialize the data structure by letting $R$ be the input portalgon $T$, and by letting the active region $R_A$ be $R$ itself, without inactive region. Repeat $N$ times the following:

1. 1.

   Apply InsertVertices to $R_A$. Then apply InsertEdges to $R_A$.

2. 2.

   Repeat 350 times the following:

   1. (a)

      Apply Gardening. Then apply SimplifyTubes to $R_A$ but in step 2 of SimplifyTubes, whenever $B$ is thin, mark $B$ as inactive. Apply Gardening again.

   2. (b)

      Apply DeleteVertices to $R_A$. Then apply InsertEdges to $R_A$.

In the end return $R$.

When proving Proposition 3, we will apply Algorithm with $N=\lceil \log (2+L/s)\rceil$.

In this extended abstract, we only sketch the proof of Proposition 8.

To analyze Algorithm from a geometric point of view, we introduce, on the segments of a flat surface $S$, a parameter that we call enclosure. So consider a segment $e$ of $S$. See Figure 4.

Informally, $e$ is “enclosed” in $S$ when a short non-contractible loop can be attached to a point of $e$ not too close to the endpoints of $e$. Formally, consider a point $x$ in the relative interior of $e$. We denote by ${⟨x⟩}_e$ the minimum length of the two sub-segments of $e$ separated by $x$. Assume that there exists a loop $\gamma$ based at $x$ in $S$, such that $\gamma$ is geodesic except possibly at its basepoint. Further assume that its length satisfies . In this case $\gamma$ and $e$ are necessarily in general position: informally, they do not overlap, more formally, every sufficiently short sub-path of $\gamma$ is either disjoint from $e$ or its intersection with $e$ is a single point. There are two cases: either $\gamma$ crosses $e$ at $x$, or $\gamma$ meets $x$ on only one side of $e$. If $\gamma$ crosses $e$ at $x$, then we say that $\gamma$ encloses $e$ in $S$. Also we say that $\gamma$ encloses $e$ by a factor of ${⟨x⟩}_e/\mathrm{\ell }(\gamma )$ in $S$. The enclosure $c_S(e)\ge 1$ is the supremum of the ratios ${⟨x⟩}_e/\mathrm{\ell }(\gamma )$ over all the basepoints $x$ in the relative interior of $e$, and over all the loops $\gamma$ based at $x$ that enclose $e$ in $S$. It is conventionally set to one if there is no loop enclosing $e$ in $S$.

The segment-happiness $h_S(e)$ and the length $\mathrm{\ell }(e)$ can be bounded from above using the enclosure $c_S(e)$. Our bound depends on the surface $S$. More precisely, on the systole of $S$ and the diameter of $S$. But instead of the diameter of $S$, we consider a triangulation of $S$, and we use its number $n$ of triangles together with the maximum length $L$ of its edges. This will be more convenient to us when analyzing Algorithm. We prove:

In Proposition 9 the $O()$ notation does not depend on $S$, it involves a universal constant. In the second inequality of Proposition 9 the exact powers above $\lceil L/s\rceil$ and $n$, here $2$ and $1$, do not matter to us. We need only a polynomial in $\lceil L/s\rceil$ and $n$.

The geometric analysis ofAlgorithm consists in two properties on the enclosure and the length of the edges involved in any execution: Lemma 10 and Proposition 12 below, whose proofs we sketch in Section 5.3.1 and Section 5.3.2. Each proof relies on properties of enclosure that are independent of Algorithm, or of any portalgon, and can be seen as independent mathematical contributions of us. In this extended abstract we only explain how these properties of enclosure serve to analyse Algorithm.

In this section, we fix a portalgon $T$ of $n$ triangles, whose sides are smaller than some positive real $L$, and whose surface $𝒮(T)$ is flat. We abbreviate $S=𝒮(T)$. We apply the algorithm Algorithm to $T$, and we discuss the execution of the algorithm.

We need a last (easy) lemma: