Scalable Ultrafast Almost-optimal Euclidean Shortest Paths

1 Introduction

The Euclidean shortest path problem (ESPP) asks for the computation of a shortest path in the Euclidean plane avoiding a set of polygonal obstacles, see Figure 1. Here, three paths from s to t are depicted with the dashed/blue one shorter than the red/dotted one and the green one. Algorithms for this problem are usually analyzed and evaluated with respect to the complexity of the problem instance, which is characterized by the number n of obstacle vertices (in this example n = 26). The importance of the problem is quite evident; for example, when computing optimal walking paths in a pedestrian zone with many open space areas, we are not restricted to an underlying network of walkable paths but can essentially roam freely in the Euclidean plane, of course avoiding obstacles like buildings. In highly automated plants, we might allow mobile robots to move not only on predefined trajectories but cross a factory hall on a beeline path, see [Mac et al., 2016] for a survey of approaches for the robotics domain. Lastly, ignoring the effects of wind and currents, routing ships on the planet's oceans could be viewed as an instance of ESPP on the sphere.

In this paper we focus on the offline variant of ESPP, i.e., we are allowed to aggregate auxiliary information about the problem instance, such that subsequent queries between vertices can be answered as fast as possible. We put special emphasis on a scalable solution which can deal with massive problem instances with several million obstacle vertices.

1.1 RelatedWork

The efficient planning of routes in road networks has enjoyed tremendous attention in the research community over the last few decades. As a result, techniques have been developed that can compute provably optimal routes within microseconds even on continent-sized networks with hundreds of millions of vertices, for which a good implementation of Dijkstra's algorithm [Dijkstra, 1959] takes seconds. See [Bast et al., 2016] for a survey. The continuous variant of the problem - ESPP - has enjoyed far less attention.

A natural way to solve ESPP is via the construction of a visibility graph, which can be computed in O(n log n + K) time where K denotes the number of edges of the resulting visibility graph, see [Ghosh and Mount, 1991]. Unfortunately, K even in practice often gets very large (in the worst case K ∈ Θ(n2 ), so visibility-graph-based solutions tend not to scale on instances with millions of obstacle vertices. Another group of algorithms based on the idea of a continuous Dijkstra is able to beat the inherent Ω(n2 ) lower bound of visibility-graph-based approaches. [Mitchell, 1996] achieved a running time of O(n3/2+δ ), for any δ > 0, which later was improved to an optimal O(n logn) in [Hershberger and Suri, 1999]. Unfortunately, these algorithms are extremely complex and we do not know of any robust implementation of them. A conceptually simple fully polynomial time approximation scheme (FPTAS) which guarantees a result of cost at most (1 + O(ϵ)) times the optimum in running time O( nϵ log 1ϵ ( 1 √ϵ + logn)) was presented in [Aleksandrov et al., 2000]. It essentially computes a very fine discretization to meet the quality guarantee and then runs Dijkstra on that discretization. More on the practical side, Polyanya [Cui et al., 2017] carefully instantiates A* on a navigation mesh and hence does not require any precomputation (apart from an underlying navigation mesh). It computes the optimal shortest path and has been shown to outperform other state-of-theart algorithms. [Shen et al., 2022] combines Polyanya with so-called compressed path databases (CPD) [Strasser et al., 2014] to obtain even faster (and still optimal) queries. Their approach requires the precomputation and storage of allpairs-information for 'interesting' (convex) vertices, though, so this approach is not applicable for instances with millions of nodes. For dynamic settings, where even the precomputation of a navigation mesh is not feasible, [Hechenberger et al., 2020] allows for efficient queries, though considerably slower than e.g. [Cui et al., 2017].

1.2 Our Contribution

In this paper we show how to adapt speed-up techniques for route planning from the network-constrained domain to the Euclidean plane. Based on a (refined) constrained Delaunay triangulation as a navigation mesh, the use of Contraction Hierarchies with local path optimization strategies allows computing paths which are less than 2% longer than the optimum within a few milliseconds even on problem instances with millions of obstacle vertices. We draw upon the OpenStreetMap project to be able to generate problem instances of almost arbitrary size and develop a lower bounding technique for very large problem instances for which optimal solutions cannot be easily computed.

2 Preliminaries

2.1 Visibility-Graph-based Solution

It is quite easy to realize that an optimal path in an ESPP instance consists of consecutive straight line segments and only changes direction at convex obstacle vertices, e.g., the green path in Figure 1 cannot be optimal. This simple observation immediately leads to an optimal (in terms of output quality) algorithm for ESPP if source and target are obstacle vertices (the case which we focus on): in a preprocessing phase we construct a so-called visibility graph which contains as vertices all obstacle vertices. For each pair v,w of mutually visible obstacle vertices, there is an edge {v,w} in the visibility graph of weight |vw|. A query between vertices s and t can then be answered by running Dijkstra's algorithm on the visibility graph. While this yields optimal paths, the main problem with this approach is its space consumption. The visibility graph can easily have Θ(n2 ) edges which makes this approach impractical for non-miniscule problem instances.

2.2 Contraction Hierarchies (CH)

Contraction Hierarchies (CH) [Geisberger et al., 2012] are a technique to accelerate shortest path planning in road networks. In a preprocessing phase, the nodes of the weighted input graph G(V, E, c) are contracted one after another. The contraction process for node v consists of removing v and all of its incident edges from G and inserting so-called shortcut edges between former neighbors u and w of v in case the only shortest path from u to w was uvw. We then call v the skip node of edge uw and the shortcut receives the weight c(uw) = c(uv) + c(vw). This ensures that at any stage of the preprocessing, the shortest path costs between all so far uncontracted nodes are still the same as in the original graph.

After all nodes are contracted, a new graph G′ (V,E ∪ E+, c) is constructed which consists of the original graph plus all shortcut edges E+ that were created during the contractions. Furthermore, the order in which the nodes were contracted is stored by assigning each node v its rank r(v) in that order. To answer shortest s-t path queries, a bidirectional Dijkstra run is conducted in the CH-graph G′ , where in the forward run from s only upward edges are considered, that is edges ab with r(a) < r(b), and in the backward run only downwards edges with r(a) > r(b). It is easy to see that this Dijkstra run settles the node with highest rank on the shortest s-t-path in both the forward and the backward run with the correct shortest path cost. If not only the shortest distance is of interest, one can store with each shortcut its two skipped edges during preprocessing and 'unfold' the shortest path after the bidirectional Dijkstra.

CH works correctly with any contraction order of the nodes. However, the speed-up over Dijkstra depends on the contraction order. There are very fast ordering heuristics which achieve excellent results in practice, we use the socalled edge difference from [Geisberger et al., 2012] which reduces query times by four orders of magnitude compared to plain Dijkstra on continent-sized road networks after a few minutes of preprocessing.

There are other speed-up techniques for road networks, yet CH is still one of the most popular ones in practice due to its simplicity. Several other (even faster) speed-up schemes use CH as their basis. See [Bast et al., 2016] for a survey. Note that CH does not work equally well for all graph instances, already more complex edge costs decrease the efficiency of CH by at least one order of magnitude in practice, e.g., see [Funke et al., 2017].

2.3 Planar Subdivisions/NavigationMesh

Having algorithms operate directly on a set of polygons in the plane seems somewhat challenging, as the set as such provides very little structure in terms of combinatorial proximity structure (e.g., determining nearby other polygon is nontrivial). Hence most approaches like [Aleksandrov et al., 2000; Cui et al., 2017] for solving the ESPP operate on a planar subdivision of the free space 'between' the obstacle polygons, sometimes also called navigation mesh. An example for such a planar subdivision/navigation map is shown in Figure 2 where a constrained triangulation of the blue polygons is depicted. A constrained triangulation is a triangulation of the polygon vertices where all polygon edges are also part of the triangulation.

3 Discretization and Acceleration

The overall idea of our approach is to first compute a 'good' navigation mesh, such that computing shortest paths on that yields reasonably short paths. As shortest path computation on the mesh (simply treating the mesh as a graph with vertices and edges) resembles the network-constrained scenario, speed-up techniques like CH might be applicable.

3.1 CDT-based Discretization

A natural choice for a navigation mesh is the Constrained Delaunay Triangulation (CDT) [Chew, 1989] which is a triangulation of the vertex set preserving the obstacle edges, see Figure 3, left. While a CDT tries to avoid triangles with very small angles, this often cannot be avoided due to the constraining edges, which in turn leads to very thin triangles as can be seen in Figure 3, left. As we want to use the triangulation graph as a basis for our path computation, this would lead to paths which are rather different from the actual shortest path which we are aiming for. One way to ameliorate this problem is to allow the insertion of additional points into the triangulation to obtain 'well-shaped', non-skinny triangles. We employ Ruppert's algorithm [Ruppert, 1995] to obtain such a triangulation. It repeatedly inserts the circumcenters of skinny triangles as Steiner points until all triangles are well-shaped, see Figure 3, right. Due to the good-natured shape of the triangles, it seems intuitive that an optimal path can be transformed into a path in this triangulation with only slightly higher cost. The well-shaped triangles also facilitate the search for the shortest path between points that are not obstacle vertices (not the focus of this work).

3.2 CH Precomputation

We first consider the shortest path problem on the triangulation only, that is, we associate with every triangulation edge which is not in the interior of an obstacle its Euclidean distance. We then compute shortest paths between vertices using, for example, Dijkstra's algorithm. To speed-up such shortest path calculations we employ Contraction Hierarchies, see Section 2.2. It is not clear a priori, though, how well CH work for these instances, since there is not a natural hierarchy as in road networks (interstates vs. motorways vs. small highways etc.). Yet, our experiments will show that a considerable speedup of orders of magnitude compared to plain Dijkstra's algorithm can be achieved, though smaller than in the road network scenario.

So in a preprocessing step we perform the CH precomputation on the refined CDT of our problem instance. This allows us to answer subsequent queries on the triangulation graph quickly and exactly. It remains to investigate howwell the resulting paths compare to the actual Euclidean shortest paths and how they could potentially be improved.

4 Path Optimization

The paths computed by the CH are quite certainly not optimal due to the loss in precision by the discretization, so in a postprocessing step we try to improve the paths computed on the navigation mesh.

The basis of these postprocessing strategies is a routine to decide whether two vertices are mutually visible. For example, in Figure 4, when trying to decide whether P1 and P2 are mutually visible, one starts traversing the triangulation along the segment P1P2, in the example visiting triangles t1, t4, and t5 before being blocked by an obstacle. P1 and P3 are seen to be mutually visible by traversal of triangles t2 and t3. The cost of such a check is essentially linear in the number of traversed triangles. We also use this routine to compute the visibility graph (which provides us with the ground truth of shortest path distances).

Given a path P1P2 . . . Pk which is optimal in the discretization, we aim at decreasing its length.

4.1 Full-Optimization

This improvement strategy computes the visibility graph restricted to the node set {P1, . . . Pk} with the visibility edges bearing their Euclidean distances. In this graph we run Dijkstra's algorithm. This clearly yields the shortest path connecting P1 to Pk amongst all paths that can only turn at {P1, . . . Pk}. Unfortunately, this does not come for free: Θ(k2 ) visibility checks have to be performed (at query time) which might be too costly, in particular for large k.

4.2 Greedy-Optimization

The greedy strategy tries to avoid excessive calls to the visibility check at the cost of a worse solution quality. Starting with P1 we determine the minimum j1 > 1 such that P1 and Pj1 are not mutually visible; clearly j1 > 2. After that we repeat on the path Pj1-1 . . . Pk etc. This clearly requires only O(k) visibility checks, each of which could still be quite expensive, though. In the experiments we see that the loss in quality is not that bad, though. For an example of the greedy algorithm exhibiting worse quality than full optimization, see Figure 5.

4.3 Funnel

In [Lee and Preparata, 1984] the authors developed an algorithm (later called funnel algorithm) to compute the shortest path between two points within a given simple n-vertex polygon in time O(n log n), which later was adapted to compute the shortest path of a given homotopy type in a polygonal domain with obstacles. In our case we can simply interpret the union of non-obstacle triangles touched by the edges of our path as the polygon and use their algorithm to compute the optimum path from P1 to Pk in time O(k log k). Note that in contrast to the previous two strategies, the resulting path might turn at vertices other than vertices of the path itself (in that sense, it is more flexible). On the other hand, the funnel algorithm can only return a path of the same homotopy type as the original path, which is not a restriction for the two visibility-based approaches. So funnel does not dominate the others, yet, the funnel algorithm has the advantage of a guaranteed running time of O(k log k) when improving a path consisting of k segments.

5 Shortcut-Optimized-CH (SO-CH) and Lower Bounds

While path optimization is quite effective in terms of quality of the solution, in particular its Full-Opt variant is quite expensive to perform at query execution, where it then dominates the overall query time as the CH-supported query is typically very fast. Hence it is natural to investigate whether the path optimization could (at least partly) be performed as a preprocessing step.

Furthermore, while we will be able to compare to the actual optimum path for smaller problem instances via the visibilitygraph-based baseline, this is becomes difficult as soon as the problem instances get larger. Yet, for those instances we would also like to have at least a provable upper bound on the deviation from the optimum.

5.1 Shortcut-Optimized CH (SO-CH)

In the course of the CH construction shortcuts are created which represent shortest paths in our CDT. If we are willing to invest some more preprocessing time, we can tune the constructed CH further as follows: Consider a shortcut (u,w). If u and w are mutually visible, we simply update the cost of the shortcut (u,w) to its Euclidean distance. To make full use of the updated distances, though, we also better update shortcuts which do not have mutually visible endpoints but contain as subpath an updated shortcut. This can be achieved by considering all shortcuts ordered increasingly according to smaller level of their vertices. When processing shortcut (u,w) we first check for mutual visibility; if so, set its cost accordingly. Otherwise set the cost of (u,w) to the sum of the costs of of (u, v) and (v, w), where (u, v) and (v, w) are the edges bridged by shortcut (u, w).

While this scheme - which we call SO-CH - increases preprocessing times, our experimental results show that the resulting path quality (even without any optimization during query execution) improves considerably.

5.2 Lower Bound Computation

In order to lower bound the length of shortest paths between two obstacle vertices, we compute the shortest path between them in a heavily simplified problem instance, as seen in Figure 6. Simplification is achieved separately for each obstacle polygon by repeatedly removing polygon vertices (and replacing the two incident edges with only one edge) until some specified removal threshold is met: 99.5% of vertices for polygons with more than 2000 vertices, 95% otherwise.

We only cut off convex corners to guarantee that any path which doesn't intersect obstacles in the original instance also does not intersect obstacles in the simplified instance. Candidate corners to be cut off are maintained in a priority queue ordered according to the area which would be removed when cutting off the respective vertex. To cut off corners in a balanced way we add the priority with which a corner was removed to its neighboring two corners. Note that in order for a corner v to be considered convex it must not only exhibit an interior angle of at most π, but the interior of the triangle defined by the vertex and its two neighboring vertices u and w must not be intersected by any of the remaining polygon edges. The latter condition is verified by traversing a constrained triangulation of the polygon from u to w and checking for intersected constraint edges.

Finally we discard all polygons with less than 4 remaining vertices and we compute a visibility graph for the simplified instance which has reasonable space consumption due to the much lower vertex count. Now a simple Dijkstra query on the visibility graph provides lower bounds.

6 Experiments

We implemented all algorithms in C++ and evaluated them on a Ryzen 9 5900x 12-core CPU/128GB RAM running Ubuntu Linux 22.04. We use the CGAL library [The CGAL Project, 2023], in particular its geometry kernel, the exact geometric predicates, as well as the constrained Delaunay triangulation code. Unless stated otherwise, averages and maxima were calculated over 1000 trials. Source code and data sets are available on a companion page [Funke, 2024].

To be able to evaluate our algorithms on instances of (almost) arbitrary size, we use a Mercator projection of the coastline polygons of the OpenStreetMap project [The OpenStreetMap Project, 2024]. The complete Planet data set contains 645,628 obstacle polygons (corresponding to continents and islands) with around 15M vertices in total. For this large graph, computing the visibility graph is infeasible, so we extracted three smaller problem instances Aegaeis, Medi(terranean), and Pata(gonia). See Table 1 for their characteristics including the visibility graph sizes (and its parallel construction time). The visibility graph for the Pata instance is relatively small due to the high density of obstacles/islands (also see Figure 7). Quite obviously it is infeasible to compute the visibility graph for the Planet instance even with considerably more powerful hardware than what was available.

6.1 Discretization

The first step of our approach is the construction of a navigation mesh. We evaluated two strategies: an 'ordinary' CDT as well as a variant with additional Delaunay refinement as described in Section 3.1. Construction of the CDT for the Planet data set took 178s, its refinement 409s (based on the CGAL code for constructing CDT in the plane). In Table 1 we also state the number of triangles before and after refinement. For example in the Pata instance, the number of triangles increased from 1.1 million to around 3.5 million after refinement.

It remains to see whether refinement - and hence increased graph size - is actually worth the effort. To analyze the effect of refinement we compared shortest paths of 1000 random source target pairs with the ground truth visibility graph, see Table 2. The maximum deviation from the optimum is considerably worse for the unrefined graph. On average, the refined graph fares considerably better for Aegaeis and Medi and slightly worse on Pata. In anticipation of the results for improvement steps (Section 6.4), choosing the refined graph is even more worthwhile. Here the refined version fares consistently better. Hence from now on, all our experiments will only use the refined triangulation.

6.2 CHAcceleration

To accelerate queries on the navigation mesh we precompute a CH on the triangulation graph (refined). Table 3 shows that the number of shortcuts created is about the same as the number of original edges in the graph. This is very similar to CH constructions on road networks, the preprocessing time is considerably higher, though; for the planet graph it takes more than 30 minutes (multithreaded); on a road network of comparable size, it would take around 5 minutes. The speedup compared to Dijkstra on the largest graph is about 3 orders of magnitude (almost 2 seconds vs 0.89 milliseconds), which is about one order of magnitude less than typically experienced in the road network context. Overall we were quite pleased with the achievable speedup, even though there is no obvious hierarchy in our setting. Yet, some vertices probably appear in quite a lot of shortest paths, e.g., the Cape of Good Hope or the Strait of Gibraltar; this might suffice to make CH

acceleration effective.

6.3 Shortcut Optimized CH (SO-CH)

Let us now evaluate the shortcut improvement scheme. Table 4 shows construction time, query time and improvement over the standard CH. First, we can see that the construction (for a given CH) is very fast when compared to the rest of the preprocessing, namely building the CH. In terms of query times, SO-CH even improves upon standard CH queries. This might be due to the fact that fewer edges have to be relaxed. Most importantly, SO-CH results in significantly shorter paths. In all cases path lengths are reduced by at least 3% on average (no ground truth for the Planet data set).

6.4 Path Optimization

A key ingredient of our technique is the improvement of the path obtained from the triangulation graph. Three schemes were presented in Section 4: Full, Greedy, and Funnel.

The quality of the results is shown in Table 5. We compare the length of the resulting (improved SO-CH) paths to the paths returned by a standard CH as well as by SO-CH, both without improvements. For example, on the Pata data set, the standard CH returned a path about 4.1% longer than the optimum on average. SO-CH without improvement lowered that already to less than 1%, Full, Greedy and Funnel reduced that even further to 0.3%, 0.3% and 0.2% respectively. Overall, in terms of average and maximum error, Funnel seems to fare best with Full getting really close. There are big differences with respect to the running times, though, as can be seen in Table 6. Even on average, Full is considerably worse than the other two approaches. This comes as no surprise as there are potentially Θ(k2 ) visibility checks required, each of which might take quite long itself. In contrast to that, Funnel guarantees a near-linear running time below 2 milliseconds even in the maximum. So overall, the Funnel improvement fares clearly best amongst the improvement strategies.

6.5 Lower Bounds

As described in Section 5.2 we construct a 'shrunken' version of our problem instance to obtain lower bounds; this takes 19 seconds for the whole planet data set, which reduced the number of vertices to around 188k - a size amenable for the construction of the visibility graph which could be done in 9 seconds (fully using 12 CPU cores). In Table 7 we see that the resulting visibility graph has around 173 million edges. The respective smaller subgraphs have accordingly fewer nodes and edges. To assess the quality of the resulting lower bounds, we sampled 1000 source-target queries and compared their lengths in the shrunken and the original visibility graph (only possible for the smaller subgraphs). While on average the lower bound is more than 90% the length of the actual shortest path, for some queries the lower bound is as bad as only 30.7% of the actual length; this mostly stems from short range queries where the relative difference becomes large as soon as one of the involved obstacle polygons is shrunken or even dropped. For queries inside a deeply nested area like Patagonia, our lower bounding technique also reaches its limits as can be seen in Figure 7. Here, in the shrunken instance, paths get considerably shorter due to shrunken or even pruned (due to small size) obstacle polygons. We want to emphasize that the shrunken planet data set was used to obtain lower bounds also for the small data sets, even though higher quality lower bounds could easily be achieved with still very moderate memory consumption by shrinking less.

6.6 Putting It All Together and ComparisonWith Existing Approaches

In Table 8 we summarize the performance of our approach (using SO-CH and the Funnel improvement) for all our data sets. Except for the largest planet data set where the construction of the visibility graph is infeasible, we compare to the actual optimum path length. For the smaller data sets the average deviation from the optimum is around 0.1% only. CH query and improvement times are together below one millisecond. For the large planet data set, we compare to the lower bound computed via the visibility graph of the shrunken polygons. Here we have an average deviation of 1.6% from the lower bound but conjecture that in reality we are much closer to the optimum. Note that when comparing to the lower bound only, the planet queries appear better than the ones on the smaller data sets. This is due to random queries typically being long range with large parts of the (near) optimal path traversing open space which can be easily be optimized by our improvement step. Even on this largest data set, a query can be answered in less than 2ms.

We ran Polyanya, a state-of-the-art solver for solving ESPP exactly, on our data sets. While running the experiments on the smaller data sets was quite hassle-free, on the largest Planet data set, we ran out of memory for many of the queries, which is why for this Planet data set, the running times were taken over only 100 queries that could be completed. For comparison: the maximum memory allocation for all steps

(precomputation and query) of our pipeline is less than 64GB of RAM. In any case, running times of Polyanya were orders of magnitude slower than our approach, which is not surprising, though, as no precomputated information (apart from the navigation mesh) can be used during the query answering phase.

In Table 10 we examined the FPTAS [Aleksandrov et al., 2000] approach on our data sets. As to be expected, the FPTAS graph gets extremely large, even for a moderate choice of ϵ = 0.5. There is a very bad dependence on small angles in the navigation mesh; they lead to a very fast increase of vertices and edges in the discretization. Depending on the instance, such small angles are unavoidable in spite of techniques like Delaunay refinement. Even for Aegaeis, the FPTAS graph contains an extremely large number of edges. The query times in Table 10 could only be achieved via an implicit representation of the graph (and applying further pruning techniques as described in [Aleksandrov et al., 2000]). Clearly, this approach is not suitable for larger graphs even though the quality of the resulting paths is reasonably good.

7 Conclusions

We presented a preprocessing scheme that allows for the ultrafast answering of approximate Euclidean shortest paths queries even for large problem instances. Our approach combines Delaunay mesh refinement with speed-up techniques from the network-constrained domain, local improvement techniques as well as a lower bounding technique, which allows us to give instance-based quality guarantees even when the optimum solution is not easily computed. Compared to state-of-the-art algorithms that can compute the optimum solution for such problem sizes, we reduce the query times from minutes to milliseconds with moderate preprocessing effort. Experiments show that our results are less than 2% above the optimum solution. For smaller problem instances we can even measure a deviation of roughly 0.1% from the optimum on average. We conjecture this to hold true for the larger problem instances as well, for which we have only a lower bound to the optimum solution. We want to note that computing Euclidean shortest paths amidst obstacles in R3 is NP-hard [Canny and Reif, 1987]. So approaches like ours which provide (provably good) approximations might be the preferred solution strategy compared to exact solutions in dimensions higher than 2.