Efficient Exact Enumeration of Single-Source Geodesics on a Non-Convex Polyhedron

Let $\mathcal{P}$ be a simplicial polyhedron in $\mathbb{R}^{3}$. When we give a source $s$ on $\mathcal{P}$, we want to compute all geodesics from $s$ to an arbitrary point $t$ on $\mathcal{P}$. Since we have infinitely many choice of $t$, we consider a query that inputs $t$ and outputs these geodesics. And yet, there are likely to be infinitely many such geodesics and we give an upperbound of their length $R$ to obtain finite result.

It is widely known that single-source shortest path problems (SSSPs) for graphs or polyhedra can be efficiently solved using a queue or a priority queue. Here, we can think of a generalization of SSSPs, namely, single-source geodesics enumeration problem on a polyhedron.

The shortest path problem on a polyhedron has been extensively researched. On a polyhedron without boundary, a shortest path is a geodesic in our sense [1], thus the shortest path problem is equivalent to the shortest geodesic problem. Shortest geodesic algorithms can be divided into exact algorithms and approximate algorithms. Since the interest of this paper is exact algorithms, approximate algorithms are not discussed in detail here. Detailed survey of this topic is given by [7, 8].

Geodesics, which are not limited to shortest, also have been researched in the context of geodesic tracing and path shortening.

Firstly, the intervals for all edges facing $s$ are generated. Figure 11 shows the cases of $s$ being inside a face (left), inside an edge (middle), and at a vertex (right). For each $I$ of these initial intervals, $I.\mbox{Parent}$ is Null, $I.\mbox{Center}$ is $s$, $I.\mbox{Depth}$ is 0, $I.\mbox{Extent}$ is the whole $I.\mbox{Edge}$, and $I.\mbox{Face}$ is the face containing both $s$ and $I.\mbox{Edge}$. For each $I$, the associated edge event and vertex event are inserted into the event queue. (Algorithm 2)

When an edge event occurs, the associated interval $I$ is projected from $I.\mbox{Center}$ into the opposite edges of $I.\mbox{Edge}$ (Figure 12). The projection result is on either one edge (left) or two edges (right). For each newly created interval $I_{i}$ (left: $I_{1}$, right: $I_{1},I_{2}$), $I_{i}.\mbox{Parent}=I$, $I_{i}.\mbox{Depth}=I.\mbox{Depth}$ and $I_{i}.\mbox{Face}$ is the triangular face in Figure 12. $I_{i}.\mbox{Center}$ is obtained by certain 3D rotation around $I.\mbox{Edge}$ to make it coplanar with $I_{i}.\mbox{Face}$. For each $I_{i}$, the associated edge event is inserted into the event queue. In the two-interval case (right), the vertex event at the intermediate vertex $v$ is associated with the interval starting at $v$ (that is $I_{2}$) and inserted into the event queue. (Algorithm 3)

When a vertex event occurs, the associated interval $I$ is recorded on $v$ as it gives the geodesic arriving at $v$. If $v$ is hyperbolic, one or more pseudo-source intervals are generated. In Figure 9, one interval in the left subfigure, two intervals in the right are generated. For each newly-created interval $I_{i}$, $I_{i}.\mbox{Parent}=I$, $I_{i}.\mbox{Center}=v$, $I_{i}.\mbox{Depth}$ is the length of the geodesic to $v$ (= the time of occurrence of this event), and the associated edge event and vertex event (if exists) are inserted into the event queue. (Algorithm 4)

After building the complete geodesic interval tree, it can accept the geodesics enumeration query which inputs a point $t$ on $\mathcal{P}$ and outputs the set $G_{st}^{R}$ of geodesics from $s$ to $t$, whose length is less than $R$. First, we must determine which intervals are responsible for the output. It can be done using the GetIntervals function:

For each obtained interval, using Algorithm 5, the geodesic is constructed from $t$ to $s$ by backtracking $I.\mbox{Parent}$. The whole procedure is given by Algorithm 6. In the pseudocode, Intersect is the function of getting the intersection, and AddFront is the operation of inserting a point at the front of the list.

For the purpose of generalizing this argument, for each hyperbolic vertex $v$, we arbitrarily choose an edge $e_{v}$ among those incident to $v$ and we fix the choice (Figure 15). It allows us to numericalize the direction of $g$ (seen from $v$) into $\alpha$ (measured along the faces). When $g$ is incoming into $v$, we call $\alpha$ the incoming angle of $g$, and when $g$ is outgoing from $v$, we call $\alpha$ the outgoing angle of $g$. Also, we use the term outgoing angle range $\iota=[\mu,\nu]$ ($\nu-\mu<\tau$) to assert that all values within this range are considered as outgoing angles, here $\tau$ is the total angle of $v$. By convention, if $\mu>\nu$ then $\iota$ is empty, and if $\nu\geq\tau$ then all values within $\iota$ are subject to the “mod $\tau$” operation.

Let us explain how it works in the example shown in Figure 16. In this figure, the five geodesics $g_{1},\dots,g_{5}$ of incoming angles $\alpha_{1},\dots,\alpha_{5}$ (respectively) arrive at the vertex in this order. Here we assume $\alpha_{1}<\alpha_{2}<\alpha_{5}<\alpha_{4}<\alpha_{3}<\alpha_{1}+\tau$. Each incoming angle $\alpha_{i}$ is mapped to the corresponding outgoing angle range $\iota_{i}=[\mu_{i},\nu_{i}]$ as follows:

1. (1)

   $\alpha_{1}$: the outgoing angle range is $\iota_{1}=[\alpha_{1}+\pi,\alpha_{1}-\pi+\tau]$ and corresponding one pseudo-source interval (not shown) is generated.

2. (2)

   $\alpha_{2}$: $\mu_{2}$ is limited by $g_{1}$ while $\nu_{2}$ is not, thus $\iota_{2}=[\alpha_{1}-\pi+\tau,\alpha_{2}-\pi+\tau]$ and two corresponding pseudo-source intervals (not shown) are created.

3. (3)

   $\alpha_{3}$: neither $g_{1}$ nor $g_{2}$ limits $\iota_{3}$, thus $\iota_{3}=[\alpha_{3}+\pi,\alpha_{3}-\pi+\tau]$.

4. (4)

   $\alpha_{4}$: $\iota_{4}$ is limited by $g_{2}$ and $g_{3}$, thus $\iota_{4}=[\alpha_{2}-\pi+\tau,\alpha_{3}+\pi]$.

5. (5)

   $\alpha_{5}$: $\iota_{5}$ is limited by $g_{2}$ and $g_{4}$ thus $\iota_{5}$ is temporarily computed as $\iota_{5}=[\alpha_{2}-\pi+\tau,\alpha_{4}+\pi]$. However, it is empty and no pseudo-source intervals are created.

Since the outline, initialization and procedure for edge events (Algorithm 1, 2, 3) are identical, we only explain the procedure for vertex events and geodesics query. All intervals generated in this way form a tree structure by the Parent pointer. we call it the reduced geodesic interval tree or the reduced GIT.

Like the naive version, when a vertex event occurs, the associated interval $I$ is recorded on $v$. If $v$ is hyperbolic, it performs the procedure of the reduction explained in the previous subsection, and, if the resulting outgoing angle range is not empty, pseudo-source intervals are generated and the associated edge events and vertex events (if exist) are inserted into the event queue. (Algorithm 7)

In this subsection, we discuss the geodesics enumeration query for the reduced geodesic interval tree, which inputs a point $t$ on $\mathcal{P}$ and outputs the set $G_{st}^{R}$ of directed geodesics from $s$ to $t$, whose length is less than $R$.

In the complete geodesic interval tree, a geodesic is given for each pair $(I,p)$ by the ConstructGeodesic function (Algorithm 5). On the other hand, in the reduced geodesic interval tree, geodesics of the same path between the last vertex $v$ and $p$ are given together for the pair $(I,p)$. The primitive geodesic between $vp$ is given by the ConstructPrimitiveGeodesic function (Algorithm 8). The geodesic is constructed from $t$ to $s$, and every time it hits a hyperbolic vertex, possible branches must be enumerated.

Let us take a look at the example illustrated in Figure 17. In this figure, each arrow indicates the direction from $s$ to $t$, although the actual query is performed backwards. Let us assume that we have just found a geodesic $g$ outgoing from $v$, and there are three geodesics $g_{1}$, $g_{2}$ and $g_{3}$ incoming into $v$, and each of $g_{1}$ and $g_{2}$ is connectable with $g$ as a geodesic while $g_{3}$ is not. Then, for each of $g_{1}$ and $g_{2}$, we check it can satisfy the length upperbound, and if so, we connect it with $g$ and perform this process recursively.

In general, we can use a simple depth-first search here. For each interval $I$, $I.\mbox{Depth}$ is the minimum of the depths of these geodesics grouped together. Since $G_{st}^{R}$ contains at least one geodesic represented by the pair $(I,t)$ if and only if the minimum of their lengths is less than $R$, the procedure of the GetIntervals function (Definition 5.2) is identical. The whole procedure is given by Algorithm 9. In the pseudocode, RemoveLast is the operation of removing the last point from the sequence, and is used for endpoints of primitive geodesics not to appear twice.

Since a geodesic is a sequence of primitive geodesics, we can reduce $G_{st}^{R}$ into a graph whose edges are primitive geodesics (Figure 18). We call it a single-pair geodesic graph, or simply a geodesic graph. It can be computed directly using the reduced geodesic interval tree.

A geodesic graph is constructed from $t$ to $s$ using Algorithm 10. We use a modified version of Dijkstra’s algorithm to determine the shortest possible length from $t$ for each primitive geodesic. We have two points to remark:

• •

  When $t$ is given, an interval can yield at most two primitive geodesics, and each primitive geodesic can be specified by the pair $(I,\mbox{\sc IsTarget})$ where $I$ is an interval and IsTarget is a Boolean value indicating whether the primitive geodesic ends with $t$ (otherwise it ends with the starting point of $I$). Note that IsTarget is set to be True only through the initialization of the algorithm since we regard $t$ to be a vertex of $\mathcal{G}_{st}^{R}$ distinct from any other vertex (Remark 6.3). Since we can assume that no duplicated intervals are supplied by the GetIntervals function, we need to check visitedness only when IsTarget is set to False (the line 11–14).

• •

  In the line 23, $l$ is the length of the primitive geodesic given by $(J,\mbox{\sc False})$. A conceptual figure is described in Figure 18. Let the red curve be the primitive geodesic and denoted as $h$. Not necessarily all geodesics in $G_{st}^{R}$ contain $h$. Let the green part and blue part be the subgraph of $G_{st}^{R}$ between $su$ and $vt$ respectively through which a geodesic containing $h$ can pass (in general the green and blue subgraph may be overlapped). Then, the distance between $su$ on the green subgraph is $J.\mbox{Depth}$, because of the construction method of the reduced geodesic interval tree. Moreover, the distance between $vt$ on the blue subgraph is $d$, because of the construction algorithm of the geodesic graph, which is derived from Dijkstra’s algorithm. Therefore the minimum length of geodesics between $st$ containing $h$ is $J.\mbox{Depth}+l+d$, and $h$ is contained in $\mathcal{G}$ as an edge if and only if it is less than $R$.

For the purpose of evaluating performance of our methods, we mainly use the Elephant mesh contained in the CGAL (Computational Geometry Algorithms Library) [15]. The proposed methods consist of two parts, i.e., the construction of a geodesic interval tree, and the geodesics query or the construction of a geodesic graph. Since the geodesics query consumes much less (often negligible) time, we only evaluate the performance of the construction of a geodesic interval tree, except Figure 34. We implemented our (single-threaded) algorithms in C++ and tested them using a machine with Intel(R) Core(TM) i9-9980XE CPU @ 3.00GHz and 128GB RAM running Linux. Except Figure 30, 30 and 34, we ran our program for 300 seconds and took statistics every second. In Figure 30 and 30, we manually recorded the amount of memory consumption measured by the OS, when it elapsed 10, 20, 30, 40, 60, 80, 100, 125, 150, 180, 210, 240, 270 and 300 seconds since each program started.

Relation of $R$ (normalized so that the mean edge length is 1) and running time is shown in Figure 25, and its log-linear plot and log-log plot are shown in Figure 25 and Figure 25. We can observe that the running time of the naive version grows exponentially to $R$, while that of improved version grows more slowly. A log-linear plot and log-log plot of the total number of generated intervals $|\mathcal{T}|$ are shown in Figure 25 and Figure 25, and they exhibit a pattern similar to that of the computation time. Figure 25 shows the slope of the log-log plot, i.e. $\Delta(\log|\mathcal{T}|)/\Delta(\log R)$, which can be considered as the exponent $\alpha$ when $|\mathcal{T}|$ is locally fit by $kR^{\alpha}$ ($k$ and $\alpha$ are constants). This value is increasing in the naive version but decreasing in the improved version, as $R$ increases. Figure 30 shows relation of $|\mathcal{T}|$ and computation time, and exhibits a pattern similar to the quasilinear growth. Memory consumption is shown in Figure 30 (to $R$) and Figure 30 (to $|\mathcal{T}|$). We can observe that the memory consumption seems to grow linearly to $|\mathcal{T}|$. Figure 30 shows the ratio of the number of propagating vertex events (vertex events that generated at least one interval) to the number of hyperbolic vertex events (vertex events that occurs at a hyperbolic vertex), among newly-processed events in one second. This value is 1 in the naive version but around 0.2 in the improved version in this instance, which means that, in the improved version, there are less vertex events that actually generate intervals compared with the naive version.

We also tested on a synthesized simple torus (Figure 31, 33, 33). We can observe that, for sufficiently large $R$, the total number of generated intervals $|\mathcal{T}|$ in the improved version is approximately $\Theta(R^{3})$. We have not established a theory, but we could guess the reason behind it as follows:

• •

  For sufficiently (very) large $R$ (and in generic cases), pseudo-source intervals are likely to “fill up” all angles around hyperbolic vertices so that no new pseudo-source intervals are generated.

• •

  In this situation, the area “swept” by the imaginary wavefront of intervals is $\Theta(R^{2})$. That is, the number of vertex events is likely to be $\Theta(R^{2})$. Since an interval in the wavefront splits into two intervals when it meets a vertex, the number of intervals in the wavefront is also likely to be $\Theta(R^{2})$. Since $|\mathcal{T}|$ can be obtained by summing the size of each generation, it is likely to be $\Theta(R^{3})$.

To evaluate relationship of smoothness of the mesh and performance, we used the recursively subdivided surfaces of Elephant using the Loop subdivision scheme [16]. In Figure 30, the computation time is shown for the original mesh (orig), and the surface subdivided once (sub1) and twice (sub2). Here $R$ is relative to the mean edge length of the original mesh. We can see that the computation time of the improved version becomes closer to that of the naive version as the mesh becomes more detailed. The reason behind it is that, pseudo-source intervals become more unlikely to overlap in the naive version, since the total angle of each vertex becomes closer to $2\pi$. The surface subdivided twice (sub2) from Elephant is used in Figure 34 to evaluate practical performance. The detail of the experiment is explained in the caption of this figure, but compared with the naive version, we can observe that the improved version took nearly half computation time and used approximately 56% memory to produce this result.