CGAL 5.4  Triangulated Surface Mesh Shortest Paths

This package provides an algorithm to compute geodesic shortest paths on a triangulated surface mesh.
The motion planning of a robot across the surface of a 3dimensional terrain is a typical application of the shortest path computation. Using a 2dimensional approximation would fail to capture anything interesting about the terrain we are trying to cross, and would give a poor solution. The problem is often called the Discrete Geodesic Problem. Although the more general version of this problem, shortest paths in 3D in the presence of obstacles, is NPHard, when the motion is constrained to the 2D surface of an object it can be solved efficiently.
The algorithm implemented in this package builds a data structure to efficiently answer queries of the following form: Given a triangulated surface mesh \(\cal{M}\), a set of source points \(S\) on \(\cal{M}\), and a target point \(t\) also on \(\cal{M}\), find a shortest path \(\lambda\) between \(t\) and any element in \( S \), where \(\lambda\) is constrained to the surface of \(\cal{M}\).
The algorithm used is based on a paper by Xin and Wang [3], a fast and practical algorithm for exact computation of geodesic shortest paths. It is an extension of earlier results by Chen and Han [1] and Mitchell, Mount, and Papadimitriou [2] .
This package is related to the package The Heat Method. Both deal with geodesic distances. The geodesic shortest path package computes the exact shortest path between any two points on the surface. The heat method package computes for every vertex of a mesh an approximate distance to one or several source vertices.
The main class of this package is Surface_mesh_shortest_path
. In the following we describe the typical workflow when using this class
The shortest paths are computed on a triangulated surface mesh, represented by a model of the FaceListGraph
concept. There is no restriction on the genus, connectivity, or convexity of the input surface mesh.
For efficiency reason, index property maps for vertices, halfedges and faces are internally used. For each simplex type the property map must provide an index between 0 and the number of simplices. We recommend to use the class CGAL::Surface_mesh
as model of FaceListGraph
. If you use the class CGAL::Polyhedron_3
, you should use it with the item class CGAL::Polyhedron_items_with_id_3
, for which default property maps are provided. This item class associates to each simplex an index that provides a \(O(1)\) time access to the indices. Note that the initialization of the property maps requires a call to set_halfedgeds_items_id()
.
The access to the embedding of each vertex is done using a point vertex property map associating to each vertex a 3D point. Defaults are provided for CGAL classes.
If the traits class used holds some local state, it must also be passed to the class when constructing it (the default one provided does not).
The set of source points for shortest path queries can be populated one by one or using a range. A source point can be specified using either a vertex of the input surface mesh or a face of the input surface mesh with some barycentric coordinates. Given a point \(p\) that lies inside a triangle face \((A,B,C)\), its barycentric coordinates are a weight triple \((b_0,b_1,b_2)\) such that \(p = b_0\cdot~A + b_1\cdot~B + b_2\cdot~C\), and \(b_0 + b_1 + b_2 = 1\).
For convenience, a function Surface_mesh_shortest_path::locate()
is provided to construct face locations from geometric inputs:
p
living in 3D space, this function computes the point closest to p
on the surface, and returns the face containing this point, as well as its barycentric coordinates;r
living in 3D space, this function computes the intersection of the ray with the surface, and (if an intersection exists) returns the face containing this point, as well as its barycentric coordinates;Usage of this function is illustrated in the example Surface_mesh_shortest_path/shortest_path_with_locate.cpp.
A time consuming operation for shortest path queries consists in building an internal data structure used to make the queries. This data structure is called the sequence tree. It will be built automatically when the first shortest path query is done and will be reused for any subsequent query as long as the set of source points does not change. Each time the set of source points is changed the sequence tree needs to be rebuilt (if already built). Note that it can also be built manually by a call to Surface_mesh_shortest_path::build_sequence_tree()
.
As for specifying the source points, the target point for a shortest path query can be specified using either a vertex of the input surface mesh or a face of the input surface mesh and some barycentric coordinates.
There are three different kinds of query functions that can be called using the class Surface_mesh_shortest_path
. Given a target point, all these functions compute the shortest path between that target point and the set of source points:
Surface_mesh_shortest_path::shortest_distance_to_source_points()
provides the closest source point to the target point together with the length of the shortest path.Surface_mesh_shortest_path::shortest_path_points_to_source_points()
provides all the intersection points of the shortest path with the edges and vertices of the input surface mesh (including the source and the target point). This function is useful for visualization purposes.Surface_mesh_shortest_path::shortest_path_sequence_to_source_points
gives access to the complete sequence of simplices crossed by the shortest path using a visitor object model of the concept SurfaceMeshShortestPathVisitor
.Some convenience functions are provided to compute:
CGAL::AABB_tree
.In short, we recommend to use a CGAL kernel with exact predicates such as CGAL::Exact_predicates_inexact_constructions_kernel
.
If you need the constructions to be exact (for the shortest path point computation for example), you should use a kernel with exact constructions. Although the algorithm uses square root operations, it will also work on geometry kernels which do not support them by first converting the kernel's number type to double
, using the std::sqrt
, and converting it back. Note that it would be preferable to use a kernel with directly supports square roots to get the most precision of the shortest path computations.
Using a kernel such as CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt
with this package will indeed provide the exact shortest paths, but it will be extremely slow. Indeed, in order to compute the distance along the surface, it is necessary to unfold sequences of faces, edgetoedge, out into a common plane. The functor SurfaceMeshShortestPathTraits::Construct_triangle_3_to_triangle_2_projection
provides an initial layout of the first face in a sequence, by rotating a given face into the xy
plane. SurfaceMeshShortestPathTraits::Construct_triangle_3_along_segment_2_flattening
unfolds a triangle into the plane, using a specified segment as a base. Since this results in a chain of constructed triangles in the plane, the exact representation types used with this kernel (either CORE::Expr
or leda_real
) will process extremely slow, even on very simple inputs. This is because the exact representations will effectively add an \(O(n)\) factor to every computation.
The following example shows how to get the shortest path to every vertex from an arbitrary source point on a surface. The shortest path class needs to have an index associated to each vertex, halfedge and face, which is naturally given for the class Surface_mesh
.
File Surface_mesh_shortest_path/shortest_paths.cpp
The following example shows how to get the shortest path to every vertex from an arbitrary source point on the surface. Note that this example uses the Polyhedron_items_with_id_3
item class. The shortest path class needs to have an index associated to each vertex, halfedge and face. Using this item class provide an efficient direct access to the required indices.
File Surface_mesh_shortest_path/shortest_paths_with_id.cpp
Although it is better to have an index built into each simplex, you can also use a surface mesh without internal indices by using external indices. The following example shows how to proceed in this case.
File Surface_mesh_shortest_path/shortest_paths_no_id.cpp
This example shows how to compute the sequence tree from multiple source points, using an iterator range of Surface_mesh_shortest_path::Face_location
objects generated at random.
File Surface_mesh_shortest_path/shortest_paths_multiple_sources.cpp
This example shows how to implement a model of the SurfaceMeshShortestPathVisitor
concept to get detailed information about the sequence of simplicies crossed by a shortest path.
File Surface_mesh_shortest_path/shortest_path_sequence.cpp
These benchmarks were run using randomly generated source and destination points over multiple trials. The measurements were executed using CGAL 4.5, under Cygwin 1.7.32, using the Gnu C++ compiler version 4.8.3 with options O3 DNDEBUG
. The system used was a 64bit Intel Core i3 2.20GHz processor with 6GB of RAM
Model  Number of Vertices  Average Construction Time (s)  Average Queries Per Second  Peak Memory Usage (MB) 

ellipsoid.off  162  0.00258805  1.21972e+06  0.39548 
anchor.off  519  0.0580262  230461  3.88799 
rotor.off  600  0.0386633  326175  3.10571 
spool.off  649  0.0418305  299766  3.75773 
handle.off  1165  0.0976167  227343  7.66706 
couplingdown.off  1841  0.138467  246833  10.1731 
bones.off  2154  0.0101125  1.31834e+06  0.865896 
mushroom.off  2337  0.206034  202582  22.5804 
elephant.off  2775  0.136177  313785  14.0987 
cow.off  2904  0.259104  206515  17.4796 
knot1.off  3200  0.279455  207084  25.314 
retinal.off  3643  0.255788  247617  29.8031 
femur.off  3897  0.25332  264825  21.4806 
knot2.off  5760  0.295655  309593  22.5549 
bull.off  6200  0.513506  209994  34.983 
fandisk.off  6475  0.609507  198768  71.3617 
lionhead.off  8356  1.23863  145810  86.6908 
turbine.off  9210  2.23755  93079.5  172.072 
man.off  17495  1.59015  187519  148.358 
Model  Number of Vertices  Average Construction Time (s)  Average Queries Per Second  Peak Memory Usage (MB) 

ellipsoid.off  162  0.00321017  911025  0.245674 
anchor.off  519  0.03601  353062  3.19274 
rotor.off  600  0.015864  805416  1.97554 
spool.off  649  0.0165743  802701  2.09675 
handle.off  1165  0.0294564  646057  4.62122 
couplingdown.off  1841  0.126045  272465  7.80517 
bones.off  2154  0.055434  536646  4.0203 
mushroom.off  2337  0.139285  290425  11.462 
elephant.off  2775  0.167269  285076  11.2743 
cow.off  2904  0.15432  328549  13.0676 
knot1.off  3200  0.114051  454640  16.1735 
retinal.off  3643  0.233208  287869  18.6274 
femur.off  3897  0.128097  457112  16.8295 
knot2.off  5760  0.413548  260195  33.484 
bull.off  6200  0.371713  297560  30.522 
fandisk.off  6475  0.545929  223865  39.5607 
lionhead.off  8356  0.70097  229449  59.6597 
turbine.off  9210  1.35703  157301  90.7139 
man.off  17495  1.75936  185194  122.541 
The following figures track the construction time, query time, and peak memory usage for the various test models as the number of source points increases. Notice that none of the values increases significantly as the number of source points increases. In fact, in most cases, the running time and memory go down. This is because a larger number of source points tends to result in a more flat sequence tree, which translates to reduced runtime and memory costs.
A geodesic curve is a locally shortest path on the surface of some manifold, that is, it cannot be made shorter by some local perturbations. On a surface mesh, this translates to a curve where, when the faces crossed by the curve are unfolded into the plane, the curve forms a straight line. Another way of describing it is that there is exactly \(\pi\) surface angle to both sides at every point along the curve (except possibly at the curve's endpoints).
A geodesic curve between two points is not necessarily a shortest path, but all shortest paths on surface meshes are formed by sequences of one or more geodesic paths, whose junction points are either vertices on the boundary of the mesh, or saddle vertices. We call such a curve on the surface of the mesh a potential shortest path between its two endpoints.
A visibility window (or visibility cone) is a pair of geodesic curves which share a common source point and enclose a locally flat region of the surface mesh. Locally flat means that between every pair of points inside the window, there is exactly one geodesic path between them which also stays inside the bounds of the window. Thus, operations, such as distance calculations, can be done with normal 2D Euclidean operations while inside the window. When a visibility window encounters a vertex (a nonflat part of the surface), a branch occurs, forming a subwindow to either side.
A saddle vertex on a surface mesh is a vertex \(v\) where the sum of surface angles of all faces incident at \(v\) is greater than \(2 \pi\), or, in simpler terms, one cannot flatten all the faces incident to \(v\) into the plane without overlap. Identifying and dealing with saddle vertices are important in shortest path algorithms because they form blind spots which cannot be reached by a single geodesic curve.
In order to deal with this, we must create a new set of child visibility windows which branch out around the saddle vertex. The paths through these child windows would first arrive at the saddle vertex, and then follow a new visibility window (forming a kind of polyline on the surface). Note that similar behavior is required when we reach a boundary vertex of a nonclosed surface mesh.
In order to compute shortest paths, we build a sequence tree (or cone tree) from each source point. The sequence tree describes the combinatoric structure of all potential shortest paths which originate from a single source point, by organizing them into a hierarchy of visibility windows.
Whenever a vertex of the surface mesh is encountered, a branch occurs in the sequence tree. If the vertex is a nonsaddle vertex, then only two children are created, one for each edge incident to that vertex on the current face. If the vertex is a saddle vertex, in addition to the two children mentioned above, a special type of node, called a pseudosource, is created which branches out from that vertex.
Once a sequence tree is built, the potential shortest paths from the source to every point inside a given visibility window can be computed. The sequence of faces along each branch of the tree are laid out edge to edge, into a common plane, such that the geodesic distance from any point on the surface to its nearest source point can be obtained using a single 2D Euclidean distance computation. Note that if the window belongs to a pseudosource, the distance is measured from the target to the pseudosource, and then the distance from the pseudosource back to its parent is measured, and so on back to the original source.
The size of the sequence tree from any source point is theoretically infinite, however we only ever care about trees which are of depth at most N, where N is the number of faces in the surface mesh (since no shortest path can cross the same face twice). Even then, the size of this truncated sequence tree is potentially exponential in the size of the surface mesh, thus a simple breadthfirst search is not feasible. Rather, we apply techniques to eliminate entire branches which are provably unable to contain shortest paths from the source point(s). The techniques used are given in greater detail in a paper by Xin and Wang [3], which itself expands an earlier work by Chen and Han [1] and Mitchell, Mount, and Papadimitriou [2] .
Handling multiple source points is simply a matter of constructing multiple sequence trees concurrently, using a method similar to the multisource Dijsktra's algorithm.
Continuous Dijkstra is simply the application of the graphsearch algorithm to a nondiscrete setting. As we build the search tree, newly created nodes are tagged with a distance metric, and inserted into a priority queue, such that the shortest distance nodes are always first.
This observation by Chen and Han states that out of all the branches that occur at any given vertex of the surface mesh, only a limited number have more than one child which can define shortest paths. This is accomplished by maintaining, for each vertex, all nodes of the sequence tree which can contain that vertex inside their visibility window.
This method alone can decrease the running time for construction of the sequence tree construction to polynomial time.
An additional distance filter proposed by Xin and Wang helps prune the search tree even further by comparing the current node's distance to the closest distance so far of the three vertices on the current face. Details on this method can be found in their paper [3].
In order to locate the shortest path from a target point to a source point, we must select the correct visibility window. A simple method is to keep track for each face \(f\) of all windows which cross \(f\). In practice, at most a constant number of windows will cross any given face, so for simplicity this is the method we employ. An alternative is to construct a Voronoilike structure on each face, where each cell represents a visibility window. We did not attempt this method, however it would seem likely that it would be of no computational benefit.
In this section we give a brief outline of the pseudocode for this algorithm. More details can be found in [1] and [3].
  Global Values  G : FaceGraph(V,E,F)  V  the set of vertices  E  the set of edges  F  the set of planar faces Q : PriorityQueue  A priority queue ordered using the metric given by Xin and Wang   Types  type VisbilityWindow: f : a face of F, the current face of this window, we say this window 'crosses' face f s : a point on the surface of F, the source point of this window d : the 'base distance' to s, only nonzero if s is a pseudosource l : the leftside bounding ray of this window, with its origin at s r : the rightside bounding ray of this window, with its origin at s p : its parent VisibilityWindow   Methods  method XinWangDistanceFilter: Input: w : a VisibilityWindow Output: filter : true if w passes the distance filtering metric given by Xin and Wang, false otherwise method PropagateWindow: Input: w : a visibility window e : an edge on face w.f Output: w' : A new visibility window on the face opposite w.f across edge e Begin: Let f' be the face on the opposite side of e as w.f Lay out face f' along e, such that it shares a common plane with w.f Create a new VisibilityWindow w', with  w its parent  the same source point and base distance as w  its boundary rays clipped to the subsegment of e covered by w return w' method CreateFaceWindow: Input: f : a face of F v : a vertex of f w : a VisibilityWindow which intersects f and contains v Output: w' : a new VisibilityWindow, with  w its parent  its source point s = v  its two bounding rays along the edges incident to v  face f as its crossed face  its base distance being the distance of window w to v method CreatePseudoSource: Input: w : the parent window v : a saddle vertex of V Begin: For each face f incident to v: w' = CreateFaceWindow(f, v, w) Q.insert(w') method TreeDepth: Input: w : a VisibilityWindow in some sequence tree T Output: The depth of node w in its current sequence tree (this would typically be cached in w itself) method ShortestPathTree: Input: s[1..n] : a set of source points on the surface of G. For simplicity (and without loss of generality), we will assume they are all vertices of G. Output: T[1..n] : a set of sequence trees for the source points Declare: O : a map of (f,v) => VisibilityWindow, which gives the 'vertex occupier' for (f,v), that is the window which crosses face f and whose source is nearest to vertex v S : a map of v => VisibilityWindow, which gives the window whose source is nearest to v. Note that this is a strict subset of O Begin: for i in 1..n: Let r[i] be the root of T[i], with distance 0 to s[i] CreatePseudoSource( r[i], s[i] ) While Q is not empty: w = Q.take() if XinWangDistanceFilter(w) and TreeDepth(w) <= F: if w contains a vertex v of w.f: if w is closer to v than O[w.f,v]: O[w.f,v] = w if v is a boundary vertex or a saddle vertex, and w is closer to v than S[v]: S[v] = w CreatePseudoSource(w, v, w.dist(v)) let {e_0, e_1} be the edges of f incident to v for i in [0,1]: w' = PropagateWindow(w, e_i) Q.insert(w') else: let e_n be the edge 'closer' to window w w' = PropagateWindow(w, e_n) Q.insert(w') else: let e_o be the only edge crossed by window w w' = PropagateWindow(w, e_o) Q.insert(w') return T[1..n]
To perform shortest path distance queries to each vertex, we can simply use the results stored in S
after the completion of ShortestPathTree, as it contains a map from each vertex to the VisibilityWindow
which has the shortest path to that vertex. Performing shortest path computations to any arbitrary face location is slightly more complex. As eluded to above, after completion of the algorithm, we traverse each of the sequence trees, and for each face f
we store all the VisibilityWindows which cross f
in a lookup structure. Then, to find the shortest path to a point on the surface of face f
, we look up that prestored set of VisibilityWindows associated with it, and among those windows we select the one which contains the query point and has the shortest path back to the origin. Though it may seem slow since it involves a linear search but it is efficient in practice since the number of faces crossing any single face is typically limited (this is due to the additional filtering method given by Xin and Wang [3].
The actual surface paths can be reconstructed by backtracking from the VisibilityWindow, through its parents in the tree up to the root, and keeping track of each face that was crossed.
This package is the result of the work of Stephen Kiazyk during the 2014 season of the Google Summer of Code. He has been mentored by Sébastien Loriot and Éric Colin de Verdière who also contributed to the documentation and the API definition.