A Survey of Algorithms for Geodesic Paths and Distances
AA Survey of Algorithms for Geodesic Paths and Distances
KEENAN CRANE,
Carnegie Mellon University, USA
MARCO LIVESU,
CNR IMATI, Italy
ENRICO PUPPO,
University of Genoa, Italy
YIPENG QIN,
Cardiff University, UKNumerical computation of shortest paths or geodesics on curved domains, as well as the associated geodesicdistance , arises in a broad range of applications across digital geometry processing, scientific computing, com-puter graphics, and computer vision. Relative to Euclidean distance computation, these tasks are complicatedby the influence of curvature on the behavior of shortest paths, as well as the fact that the representation ofthe domain may itself be approximate. In spite of the difficulty of this problem, recent literature has developeda wide variety of sophisticated methods that enable rapid queries of geodesic information, even on relativelylarge models. This survey reviews the major categories of approaches to the computation of geodesic pathsand distances, highlighting common themes and opportunities for future improvement.CCS Concepts: •
Mathematics of computing → Graphs and surfaces ; •
Theory of computation → Computational geometry ; •
Computing methodologies → Shape modeling ; Scientific visualization .Additional Key Words and Phrases: datasets, neural networks, gaze detection, text tagging
ACM Reference Format:
Keenan Crane, Marco Livesu, Enrico Puppo, and Yipeng Qin. 2020. A Survey of Algorithms for Geodesic Pathsand Distances. 1, 1 (July 2020), 56 pages. https://doi.org/10.1145/nnnnnnn.nnnnnnn
The ability to rapidly compute shortest paths and/or distances between pairs or sets of points is afundamental operation in computational science. This survey considers algorithms for computing geodesic paths and distances, i.e. , paths through curved domains that arise in a broad range of tasksacross digital geometry processing, scientific computing, computer graphics, and computer vision.Relative to computing distance on Euclidean domains, this problem is complicated by the influenceof curvature, as well as the fact that the domain itself may not be known exactly.Several problems arise in this context, which have different applications, and are best tackledwith different techniques: from tracing a geodesic line given a starting point and a direction; tocomputing the shortest path between a pair of points; to computing the distance function overa whole surface, with respect to a source point or set; to provide a framework to rapidly extractthe shortest path between any pair of points. Broadly speaking, there are two major classes ofmethods: those rooted in computational geometry, which view a polyhedral surface as an exactdescription of the geometry, and those rooted in scientific computing, which view a polyhedron asan approximation of a smooth surface. Neither class of approaches is universally “better”: each maybe best-suited to a particular task ( e.g. , finding geodesics on a CAD model, versus approximatinggeodesic distance on a scanned surface), and in a particular setting, according to a wide variety oftrade offs in terms of accuracy, storage cost, run time performance, and scalability.In general one might wish to compute geodesics and geodesic distances on many differentdata structures (point clouds, voxelizations, etc. ), though in this survey we focus primarily on
Authors’ addresses: Keenan Crane, Carnegie Mellon University, Pittsburgh, USA, [email protected]; Marco Livesu,CNR IMATI, Genoa, Italy, [email protected]; Enrico Puppo, University of Genoa, Genoa, Italy; Yipeng Qin, CardiffUniversity, Cardiff, UK.2020. XXXX-XXXX/2020/7-ART $15.00https://doi.org/10.1145/nnnnnnn.nnnnnnn , Vol. 1, No. 1, Article . Publication date: July 2020. a r X i v : . [ c s . G R ] J u l Crane, et al. polyhedral surfaces , represented as triangle meshes. The survey is organized according to thedifferent geodesic distance problems and the attendant classes of approaches to their solution(see Table ?? for a summary). After introducing the basic problems in Sec. 1.1 and providing thenecessary background in Sec. 2, we review major classes of algorithms in Secs. 3, 4 and 5. In Sec. 6,we examine the relationship between mesh quality and geodesic computation. In Sec. 7, we providea partial evaluation of the reviewed methods, on the basis of available implementations. Finally, inSec. 8, we draw some conclusions and we highlight common themes and opportunities for futureimprovement in geodesic algorithms. Previous surveys.
Our focus in this survey is on practical algorithms and their behavior onreal-world datasets. A relatively recent survey by Bose et al. [15] deals primarily with theoreticalaspects of one particular problem (the polyhedral shortest path problem), with a focus on asymptotictime complexity and approximation bounds, as well as special cases such as convex polyhedra.Less attention is devoted to broader problems (such as geodesic distance transforms) or issues suchas real-world performance, mesh quality, etc. ; in this respect, such a survey is complementary toours. Patané [76] provides a detailed survey on the specific topic of Laplacian spectral kernelsand distances, which is likewise complementary to our account of PDE-based methods. Peyré andcolleagues [79, 80] provide a nice overview of several aspects of PDE-based methods, though thelast few years have seen major advancements in PDE-based methods, which we discuss in Sec. 3.
Tasks in geometry processing may require a variety of different queries about geodesics andgeodesic distances; though seemingly similar, these queries must be answered via very differentalgorithms. We review several important cases here.
Initial Value Problems.
Perhaps the simplest type of problem is the initial value problem , whichtraces out a geodesic starting at a given point in a given direction. In the Euclidean case, thesolution is simply a straight ray through space, though in general one must of course account forthe influence of curvature; in the discrete setting one also encounters some difficulty in definingwhich ray is most “straight” (see Section 2.2). Computationally this problem, which we will refer toas
Geodesic Tracing (GT) , tends to be among the easiest of geodesic queries; we examine it in detailin Sec. 5.
Boundary Value Problems. In boundary value problems , a set of points on the domain is given asinput, and one seeks (for instance) geodesic paths between these points, or geodesic distances to allthese points. This type of problem is generally not as easy as an initial value problem since, forinstance, one cannot simply shoot a geodesic ray from one point and hope to hit a particular targetpoint. We consider in particular the following problems: • Point-to-Point Geodesic Path (PPGP): given two distinct points p , q , find any shortest geodesic γ between them. (Note that in general the solution may not be unique—see Sec. 2.) • Single Source Geodesic Distance / Shortest Path (SSGD/SSSP): given a source point p , computea scalar function φ ( q ) that provides the length of the shortest path from p to each point q . InSSGD, only the distance function is provided; in SSSP, the shortest paths are also explicitlycomputed. • Multiple Source Geodesic Distance / Shortest Path (MSGD/MSSP): in this case the source isno longer just a single point p , but rather a collection S of points, curves, etc. . The result is afunction ϕ ( p ) which for each point p gives the distance to the closest point in the set; thisquery is also sometimes called a distance transform . In MSSP, paths are also provided. • All-Pairs Geodesic Distances / Shortest Paths (APGD/APSP): given a source set S , computethe shortest distance/paths between all pairs of points in S . , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 3
In Sections 3 and 4, we review the methods for solving these problems, categorized according totheir basic algorithmic approach they take. Most methods that solve the SSGD/SSSP problems alsoprovide a solution to PPGP as a byproduct; several easily generalize to MSGD/MSSP; only fewmethods address the APGD/APSP problems. In Section 4.2, we review some algorithms that arespecifically developed to solve PPGP.Note that most methods in the literature assume that initial points or boundary conditions arespecified at vertices of an input triangulation, though some methods allow direct evaluation ofgeodesic queries at arbitrary points of the domain. For algorithms that only support queries atvertices, a pragmatic approach is to locally refine the mesh by splitting faces or edges at the querypoint.
When thinking about algorithms for computing geodesics, it is important to consider what ourdomain represents: does it exactly represent the geometry of interest, or is it merely an approxima-tion of the true domain? The answer to this question depends on context. For instance, in problemslike animation or CAD/CAM, where surfaces are designed by artists or engineers, we may havean exact boundary representation and want to compute exact paths along this surface. On theother hand, when a mesh is obtained via, say, 3D scanning or numerical simulation, it can onlyprovide an approximation of the true domain, i.e. , the real physical surface of interest. In this case,the accuracy of geodesic computation is fundamentally limited by the accuracy of the domainrepresentation itself. The choice of application therefore influences how we talk about error, andalso which algorithms are best suited to a given task. abc d
In the context of polyhedral surfaces, a common misconception is that the“correct” answer is obtained by considering piecewise linear paths through thedomain – though of course these paths may only be approximations of true(smooth) geodesic curves. A simple 1D example is illustrated in the inset. Whatis the “correct” distance between the point a and the point c ? If we view thepolygon as an exact representation of the geometry ( i.e. , if we wish to computedistance along the square itself), then the geodesic distance is obtained bysumming the two edge lengths | ab | + | bc | = √ ≈ .
82. On the other hand,if we view the polygon as an approximation of the smooth circle, the distance should just be thearc length between a and c , i.e. , π ≈ .
14. Not surprisingly, the polygonal geodesic distance is anunderestimate of the smooth geodesic distance, since each straight segment takes a “short cut” frompoint to point. This same phenomenon occurs on triangulated surfaces: the polyhedral geodesicdistance generally underestimates the smooth geodesic distance (for instance, the distance alonga cube provides a poor approximation of the distance along the circumscribed sphere). To makethe distinction clear in this survey, we therefore distinguish between two versions of the geodesicproblem: • Polyhedral Geodesic Problem — view a discrete surface as an exact description of thegeometry, and aim to compute exact geodesics or exact distances along this polyhedralsurface. • Smooth Geodesic Problem — view a discrete surface as an approximation of the truegeometry, and aim to compute the most accurate possible approximation of geodesics orgeodesic distances.The polyhedral problem captures the standard viewpoint of computational geometry; the smoothproblem encapsulates the standard viewpoint of scientific computing. These two problems interactof course: for instance, one can use the exact polyhedral distance to approximate the true smooth , Vol. 1, No. 1, Article . Publication date: July 2020.
Crane, et al.
Fig. 1. When a polyhedral surface is an approximation of a smooth surface, the “exact” polyhedral distancestill does not recover the true distance of the smooth surface. For instance, even for very nice triangulationsof the smooth sphere, algorithms from computational geometry improve accuracy by only one order relativeto PDE-based methods (from linear to quadratic). distance. Surprisingly enough, however, the best strategy for approximating the smooth geodesicdistance may not be to simply compute the exact polyhedral distance. For instance, even on a nicelytriangulated sphere, the polyhedral distance gives only an O ( h ) approximation of the smoothdistance (Fig. 1); greater accuracy can be achieved by considering the distance along a spline orsubdivision surface that approximates the same sampling of the domain [30, 72]. The fact thatone can do better than the exact polyhedral distance should not come as a surprise: higher-ordergeometric elements ( e.g. , spline or subdivision patches) use a much larger stencil of information thanindividual triangles, which correspond to linear elements (consider, for instance, approximating thecircle by four Bézier curves rather than four straight segments). On the other hand, higher-orderinterpolation provides no benefit for the polyhedral geodesic problem, where one is interested inthe exact distance along a particular piecewise linear surface. In this section we provide some basic concepts and definitions that will facilitate our discussion ofalgorithms—for a more thorough introduction, see [36] in the smooth setting (especially Chapter4, Sections 4–4, 4–6, and 4–7), and [26] for the discrete setting. At a high level, geodesics canbe characterized as curves that are in some sense straightest and (locally) shortest , though onemust be careful about the relationship between these two characterizations (Sec. 4.1), especially onpolyhedra where they may lead to different discrete definitions [81]. We first give some standardbackground in the smooth setting (Sec. 2.1), followed by a discussion of geodesics on polyhedralsurfaces (Sec. 2.2).
Intrinsic vs. Extrinsic Geometry.
When studying geodesics, how should we describe the shapeof the domain? An important feature of geodesics is that they depend only on distances along thesurface, and not at all how the surface sits in space. As a concrete example, consider a piece ofpaper rolled up into a tube (Fig. 2): a straight line ℓ drawn between two nearby points p and q onthe flat piece of paper is still a shortest curve on the tube (see Fig. 2). One can find an even shorterpath by drawing a straight line through space, but this line is no longer a path along the surface.Likewise, there are many other ways we could bend the tube without changing the shape or lengthof shortest paths—for this reason, when studying geodesics we really want to ignore the extrinsic , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 5
Fig. 2. A geodesic is any “straight” curve between two points p and q , and are a feature of the intrinsicgeometry of the surface: they do not change if we bend the surface, so long as there is no stretching, shearing,or ripping. Though the shortest path will always be a geodesic (here, ℓ ), a geodesic is not necessarily theshortest path (consider ℓ ). geometry of the surface (how it sits in space), and focus purely on its intrinsic geometry (onlythose things that can be measured by walking along the surface). The desire to compute shortestintrinsic (rather than extrinsic) paths is of course quite natural: for instance, the “shortest” routebetween two cities is the one that best avoids mountains and valleys—not the one that tunnelsstraight through the Earth! Shortest vs. Straightest Curves.
The tube example also helps to illustrate another feature ofgeodesics: shortest implies straight, but straight does not necessarily imply shortest. Consider, forinstance, a long straight line ℓ from p that “goes off the edge of the map” before returning to q .This curve is still a geodesic, but it is not a shortest geodesic —in particular, it is not as short as ℓ .Finally, the word “shortest” does not imply that there is a unique geodesic of minimal length. Forinstance, there may be two equally short ways to go around a hill; in fact, there may be infinitelymany shortest paths between two given points, such as the north and south pole of a sphere.Unfortunately, finding geodesics is not as simple as just “unrolling” a smooth surface into theplane and finding straight lines (as we did with the tube), since most surfaces cannot be flattenedwithout distorting distances. Consider for instance the many different map projections used bycartographers (Mercator, Robinson, etc. ), none of which preserve geodesic distances. However, wecan nonetheless reason about shortness and straightness of curves: for instance, a curve on anembedded surface is “straight” if there is no acceleration as we travel along the curve, except in thenormal direction—in other words, if we only turn by the bare minimum amount needed to remainon the surface. Alternatively, we can adopt the perspective of Riemannian geometry , which allowsone to reason about intrinsic geometry without thinking about how it is embedded in space.
In the intrinsic picture, our main object of study is a smooth surface M with Riemannian metric д . The fact that M is smooth implies that at each point p ∈ M we havea tangent space T p M , where each tangent vector X ∈ T p M specifies a vector pointing “along” thesurface, i.e. , all possible directions we can travel away from p . Since we have no information abouthow the surface sits in space, the one and only way to measure lengths and angles of tangentvectors is via the Riemannian metric, which for tangent space T p M provides a positive-definiteinner product д p : T p M × T p M → R ≥ . For instance, the length of any tangent vector X ∈ T p M is given by | X | : = (cid:112) д p ( X , X ) ; the angle between any two tangent vectors X , Y ∈ T p M is givenby arccos ( д p ( X , Y )/| X || Y |) , just as in ordinary Euclidean R n . In fact, whenever the surface M is sitting in space, the Riemannian metric and Euclidean inner product coincide. When the meaningis understood from context, or when working with a vector field X ( i.e. , a choice of vector in eachtangent space), we can drop the subscript p . , Vol. 1, No. 1, Article . Publication date: July 2020. Crane, et al.
The Riemannian metric enables us to easily define geodesic curves in termsof conditions on length. In particular, consider any differentiable map γ from an interval [ a , b ] ofthe real line into the surface M , and let (cid:219) γ ( t ) : = ddt γ ( t ) ; we say that γ is arc-length parameterized if | (cid:219) γ ( t )| = t ∈ [ a , b ] . We can express the length of any curve γ as L ( γ ) : = ∫ ba | (cid:219) γ ( t )| dt = ∫ ba (cid:113) д γ ( t ) ( (cid:219) γ ( t ) , (cid:219) γ ( t )) dt . The geodesic distance ϕ ( p , q ) between any two points p , q ∈ M is the infimum of length over allcurves γ such that γ ( a ) = p and γ ( b ) = q . An arc-length parameterized curve γ is a shortest geodesic if it achieves this infimal length. An arc-length parameterized curve γ is a geodesic if it is “locallyshortest,” i.e. , if there is some δ > γ is a shortest geodesic when restricted to any interval [ t , t ] ⊂ [ a , b ] such that t − t ≤ δ . Intuitively, a shortest geodesic is a curve that cannot be madeany shorter (without moving its endpoints); a geodesic is a curve that cannot be made shorter byadjusting any small piece of it. We can also characterize geodesics in terms of straightness rather than shortness . This perspective is best understood from a dynamic point of view: imagine we are drivingalong a road at constant speed, and never turn our wheel left or right. Although we may encounterhills and valleys along the way, our trajectory is in some sense as straight as it possibly could be. Inthe extrinsic setting the only acceleration we experience is in the direction normal to the surface, i.e. , the minimal possible acceleration needed to keep our vehicle on the ground. In the intrinsicsetting, we do not need to worry about this normal acceleration since our surface is not sitting inspace—instead, we just say that a geodesic is a curve of zero acceleration.
How can we measure the acceleration of a curve γ ( t ) in a surface M ? In the Euclidean plane, wecan just take the second derivative of γ with respect to t . In the Riemannian setting, however, thevelocity vectors γ ′ ( t ) and γ ′ ( t ) at two different times t , t will in general belong to two differenttangent spaces T γ ( t ) M , T γ ( t ) M , and hence cannot be compared directly. Instead, we can use the Levi-Civita connection or covariant derivative ∇ X Y , which (intuitively) uses the metric д to measurehow much Y is turning as we move along the direction X . In particular, a curve γ is a geodesic if itstangent (cid:219) γ exhibits zero turning as we move along the tangent direction, i.e. , if ∇ (cid:219) γ ( t ) (cid:219) γ ( t ) = t ∈ [ a , b ] . (Since the covariant derivative is local, it does not help us define shortest geodesics.) More generally, for any arc-length parameterized curve γ ( s ) we have ∇ (cid:219) γ ( s ) (cid:219) γ ( s ) = κ д ( s ) n ( s ) , where for each time t , the vector n ( t ) ∈ T γ ( t ) M is the unit normal to the curve ( i.e. , a unit vectororthogonal to the tangent), and the scalar κ д ( t ) ∈ R is the geodesic curvature , which measures therate at which the curve is bending. A geodesic is then a curve with zero geodesic curvature .We can also understand geodesic curvature from an extrinsic point of view. Suppose we have amap f : M → R assigning each point of the surface M a location in space, and let N : M → S ⊂ R be the corresponding Gauss map giving the unit normal at each point. A curve γ : [ a , b ] → M on thesurface can now be realized as a curve ¯ γ : = f ◦ γ in Euclidean space. Assuming ¯ γ is parameterizedby its arc-length s , its unit tangent is given by ¯ T ( s ) : = dds ¯ γ ( s ) ; we will use n : = T × N to denotethe normal to the curve (as opposed to the normal N of the surface). The derivative of T in turndescribes the curvature of the curve, which can be decomposed into two scalar components: κ д : = ⟨ n , dds T ⟩ , κ N : = ⟨ N , dds T ⟩ . , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 7
Fig. 3. Extrinsically, the curvature of a curve ¯ γ can be decomposed into the normal and geodesic curvatures κ N , κ д , which measure the change in the tangent in the direction of the surface normal N and the curvenormal n , respectively (left). A “wiggly” but planar curve will have zero normal curvature and nonzero geodesiccurvature (center), whereas a great arc on a sphere has zero geodesic curvature but nonzero normal curvature(right). Note that there is no derivative in the T direction, since ¯ γ is arc-length parameterized and hence itstangent doesn’t change length. In other words, the geodesic curvature κ д describes how much thecurve γ is bending “in plane”; the normal curvature κ N describes how much γ is bending “out ofplane.” Geodesics are again curves of zero geodesic curvature, i.e. , those that go straight along thesurface (and bend purely to remain on the surface). An illustration is provided in Fig. 3. At any point p ∈ M and any unit vector X ∈ T p M , there is a unique geodesictraveling away from p in the direction X , i.e. , such that (cid:219) γ = X (this perspective leads to the tracingalgorithms discussed in Sec. 5). More generally, the exponential map exp p : T p M → M takes anytangent vector X at p to the point q of M reached by walking in the direction X /| X | a distance | X | .In general, however, this map is not invertible: there can be distinct points q (cid:44) q that are reachedby starting at p and walking along geodesics in different directions (or for different distances). Anyneighborhood U of p over which exp p is invertible is called a normal neighborhood of p ; on any suchneighborhood, the inverse of the exponential map defines the logarithmic map log p : U → T p M .The logarithmic map provides local coordinates around p called normal coordinates , by effectively“flattening out” a small patch of the surface onto its tangent space. Polar coordinates on the tangentspace are sometimes called geodesic polar coordinates ; here, geodesics can be characterized (locally)as lines of constant angle. Circles in this tangent plane are mapped by the exponential map togeodesic circles, which orthogonally intersect geodesic rays from p (Gauss’ lemma).Away from a normal neighborhood, geodesics may start behaving in less intuitive and undesirableways. We review some global aspects of geodesics in Sec. 4.1. Many of the algorithms discussed in this survey aim not to computeindividual geodesic curves, but rather the geodesic distance function ϕ : M × M → R over the entiresurface. In particular, for any two points x , y ∈ M , ϕ ( x , y ) is the smallest length of any geodesicbetween x and y . This function satisfies all the properties of a distance metric, i.e. , nonnegativity( ϕ ( x , y ) ≥ ϕ ( x , y ) = ⇐⇒ x = y ), symmetry ( ϕ ( x , y ) = ϕ ( y , x ) ) and triangleinequality ( ϕ ( x , y ) + ϕ ( y , z ) ≥ ϕ ( x , z ) ); these properties will be important to keep in mind whenstudying numerical approximations of geodesic distance. A geodesic ball B p ( r ) is the set of all points q ∈ M such that ϕ ( p , q ) < r . Cut locus.
For a given point p , the injectivity radius is the radius r > B p ( r ) such that there is a unique shortest path from any point q ∈ B p ( r ) back to p . Outside thisradius, there are points q with two or more shortest paths to p . The collection of all such points is , Vol. 1, No. 1, Article . Publication date: July 2020. Crane, et al. called the cut locus . More generally, for any source set Ω ( e.g. , a collection of isolated points, or anetwork of curves, surfaces, etc. ) the cut locus is the set of points q for which there is not a uniqueshortest path to some point in Ω . The geodesic distance function ϕ ( x , y ) is smooth away from thecut locus. In most situations, it is not particularly important that geodesic distance algorithmsaccurately reconstruct the cut locus: only that the distance values or paths computed at each pointare accurate. However, some applications ( e.g. , computing the medial axis of a shape) may requirean accurate cut locus. Geodesics on polyhedral surfaces behave somewhat differently from geodesics on smooth surfaces—especially in the vicinity of vertices, where the surface fails to be smooth. Methods coming fromcomputational geometry (Sec. 4) must therefore carefully consider the specific behavior of polyhe-dral geodesic paths; methods based on the finite element viewport (Sec. 3) largely side-step thesequestions by approaching geodesic computation from the perspective of function approximation rather than explicit path tracing.A polyhedral surface is, roughly speaking, a collection of Euclidean (planar)polygons glued together to form a surface. Since any planar polygon canbe triangulated (and the choice of triangulation has no effect on geodesics),we will assume that the combinatorics of any polyhedron are encoded by asimplicial complex M with vertices V , edges E , and faces F . For simplicity wealso assume that this complex is manifold , i.e. , the link Lk ( i ) of every vertex i ∈ V interior vertex is a single closed loop; the link of every boundary vertexis a single path (see inset figure).The geometry of a polyhedral surface can be specified in one of two ways: either extrinsically ,using vertex coordinates f : V → R n which are linearly interpolated over each triangle, or intrinsically , by specifying positive edge lengths ℓ : E → R > that satisfy the triangle inequalityin each face: ℓ ij + ℓ jk ≥ ℓ ki for all ijk ∈ F . Of course, one can always obtain edge lengths from Fig. 4. Intrinsically, a vertex i of a polyhedral surface looks like either a circular cone, a flat disk, or a saddle,depending on the sign of the angle defect Ω i . , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 9 vertex coordinates ( ℓ ij : = | f j − f i | ), though the vast majority of geodesic algorithms can in principlebe implemented without reference to vertex coordinates—reflecting the fact that geodesics are apurely intrinsic object. Instead, quantities like triangle areas A ijk and interior angles θ jki at trianglecorners can be computed directly from edge lengths (using expressions like Heron’s formula and thelaw of cosines). This fact is especially useful given that in some geometric problems ( e.g. , manifoldlearning) one may have distances between points, but not an explicit embedding in R n .For each vertex i ∈ V , the angle defect is the quantity Ω i : = π − (cid:213) ijk ∈ F θ jki , where θ jki denotes the interior angle at vertex i of triangle ijk , and the sum is taken over all triangles ijk containing vertex i . Intuitively, this quantity captures the “flatness” of the vertex, and is oftenviewed as a discrete analogue of the Gaussian curvature. An important mental image is provided inFig. 4: imagine that the triangles around a vertex i ∈ V are flat pieces of paper glued together alongedges. Depending on the value of Ω i , these triangles can be bent into a smooth circular cone, a flatcircular disk, or a saddle-like figure, all without changing geodesic distance. We will therefore referto vertices with positive, zero, and negative angle defect as spherical or cone-like , Euclidean , and hyperbolic or saddle-like , respectively. This local picture nicely encapsulates the intrinsic geometryof any polyhedron: it is smooth and intrinsically flat away from the vertices—even the location of edges is irrelevant when it comes to thinking about geodesics and geodesic distance, since the edgesare effectively invisible to an intrinsic observer. The sign of Ω plays an especially important rolein the context of polyhedral geodesics, since shortest geodesics will often pass through verticeswhere Ω i <
0, but can never pass through vertices where Ω i > i.e. , by finding an arrangement of vertices in R thatagrees with the intrinsic edge lengths—several examples are shown in Fig. 4. This picture makesit clear that, locally , geodesics on polyhedral surfaces can be constructed by simply drawing linesegments in the Euclidean plane. The main computational challenge, therefore, is answering more global questions: for example, which sequence of triangles must we unfold to find the shortest suchline? An analogous perspective is not generally not available for smooth surfaces, since any localflattening will invariably distort lengths ( i.e. , geodesics are rarely straight lines in local coordinates).On the other hand, the fact that our surface is no longer smooth makes the definition of polyhedralgeodesics somewhat more nuanced—especially in the vicinity of vertices. In the smooth setting we had two equivalent characterizations ofgeodesic curves: they are both straightest (Sec. 2.1.3) and locally shortest (Sec. 2.1.2). As studied byPolthier & Schmies [81], these two characterizations are no longer equivalent in the polyhedralsetting. This situation reflects a common “no free lunch” scenario in the discretization of objectsfrom differential geometry [27], namely that one typically cannot find a single definition (in thiscase, for discrete geodesics) that exactly captures all the key properties of the original smoothobject (in this case, both straightest and locally shortest).Locally, polyhedral geodesics essentially behave the same as in the smooth setting. Consider forinstance a pair of points p , q contained in the same triangle—here geodesics are just ordinary linesegments, which are both shortest and straightest. Likewise, for two points p , q close to a commonedge (and away from any vertex) we can simply unfold the two incident triangles into the planeand connect them by the unique shortest, straightest segment (see Fig. 5, left). Globally, however,the situation is more complicated due to non-smooth points at vertices. , Vol. 1, No. 1, Article . Publication date: July 2020. Straightest.
To find the straightest path leaving a point p ∈ M in a direction u ∈ T p M , we cansimply apply the local observations made above: inside a given triangle the shortest path is foundby extending a straight ray along u ; to continue this path into the next triangle we can imagineunfolding neighboring triangles into the plane and extending this ray into the next triangle. Theresulting path corresponds to a straight line along a strip of planar triangles, as depicted in Fig. 5,right. (Note that for very long paths we may encounter the same triangle more than once, inwhich case we would have multiple copies of this triangle in the unfolding.) This tracing operationeffectively defines a discrete version of the exponential map discussed in Sec. 2.1.4.What should we do if our path enters a vertex i ∈ V ? In particular, whichoutgoing direction describes the “straightest” curve? Unless the angle defect Ω i isequal to zero, we cannot simply unfold triangles into the plane. An idea consideredby Polthier & Schmies [81] is to instead pick the outgoing direction such that thereis “equal angle” on either side. More precisely, let u be the incoming direction; wecan define the outgoing direction v as the one such that the total angle between u and v is exactly half the sum Θ i : = (cid:205) ijk θ jki of interior angles θ jki at vertex i (seeinset figure). Equivalently, one can work in a local polar coordinate system whereangles θ jki are normalized to sum to 2 π ; this viewpoint has been carefully studiedby Troyanov [95], and was later adopted in geometry processing for problemsinvolving polyhedral geodesics [81] and tangent vector fields at vertices [111];Sharp et al. [91, Section 5.2] provides a concise description. As in the smooth case,this definition of straightness yields a unique solution to the initial value problem, even for pathsthrough vertices with nonzero angle defect. Shortest.
In contrast, the behavior of shortest curves on a polyhedral surface depends critically onthe sign of the angle defect Ω i . Consider for instance two points a , b directly opposite a cone-likevertex ( Ω i > γ from a to i , then from i to b .However, one can find an even shorter path by walking “around” the cone—just as one might find ashorter path by walking around a hill, rather than walking over it. In particular, if we cut the conefrom the boundary through the point a and up to the vertex i , then the straight line segment from b to either copy of a gives us a shortest path. Hence, straightest curves on a polyhedral surface arenot always shortest. In this case, γ is not even locally shortest, since even in a tiny region around i would can make it slightly shorter by going around vertex i , instead of through it. In general, it is never be advantageous (in terms of path length) to pass through a spherical vertex. Therefore, theonly shortest geodesics passing through a positive vertex i are those terminating at i .On the other hand, if i is a saddle vertex ( Ω i < ℓ from a point p ∈ M to a saddle vertex i ,there will be infinitely many outgoing directions v that yield a shortest path;these directions form a wedge-like region of angle | Ω i | around the incomingdirection u (see inset figure). The union of the incoming path ℓ with the wedgeis sometimes called a funnel , and is the starting point for a large family ofalgorithms reviewed in Sec. 4.1. Again the “shortest” and “straightest” picturesdisagree: a straightest geodesic passing through a saddle vertex must bisect thefunnel, whereas a shortest geodesic can include any from i contained in thefunnel. , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 11
Fig. 5. Any pair of adjacent triangles can be unfolded into the plane without distorting distance (left). Ageodesic on a polyhedral surface is therefore equivalent to a straight line along some planar triangle strip—solong as it does not pass through any vertices.Fig. 6. The triangles around any vertex i with positive angle defect Ω i > (left) are intrinsically no differentfrom a smooth cone (center). Therefore, the shortest path between two points a , b on a polyhedral surface willnever pass through a positive vertex: we can always cut open the cone in a way that lets us draw a straightline from a to b without passing through i . When i sits directly between a and b , this path will not even beunique (right). The analysis above has important implications for the exponentialmap on a polyhedral surface. Consider for instance tracing straightest geodesics in every possibledirection from a vertex i ∈ V . When paths hit a spherical vertex j , they split into two groups whichmeet discontinuously on the opposite side of j . When paths hit a saddlevertex k , they again split into two groups, which do not completely coverthe funnel opposite k . In either case, the exponential map fails to be injectiveas soon as we hit any non-flat vertex; in other words, the injectivity radiusis simply the distance to the closest vertex. Moreover, every spherical vertexof a polyhedral surface is contained in the cut locus. As explored by Itoh &Sinclair [47], this means that algorithms which exactly compute distance ona polyhedral surface will (surprisingly enough) do a very poor job of approx-imating the cut locus of the smooth surface it approximates. For instance, thecut locus on a polyhedral sphere will contain every single vertex , whereas the cut locus of a smoothsphere consists of just a single point. A large number of methods for computing geodesic distance are based on formulating the problemin terms of partial differential equations (PDEs) on a smooth manifold, then discretizing and solvingthese PDEs via, e.g. , finite element methods (FEM) or other numerical techniques. These methodsare generally suitable for computing the single or multiple source geodesic distance (SSGD/MSGD);explicit geodesic can subsequently be obtained by, e.g. , tracing curves through the gradient ofthe distance function (though such tracing requires careful consideration in order to achieveaccurate results; see Sec. 5). Some of these methods are also quite attractive for solving all-pairsgeodesic distance problems (APGD), since their solutions are well-approximated by precomputing a , Vol. 1, No. 1, Article . Publication date: July 2020. relatively small collection of eigenfunctions (Sec. 3.3.1). PDE-based methods are attractive becausethey are built on top of well-established techniques from scientific computing (such as FEM), asopposed to computational geometry methods (Sec. 4) which must be built “from the ground up.”As a result, such techniques often come with a clearer picture of, e.g. , convergence rates underrefinement, making them particularly well-suited for the smooth geodesic problem (Section Sec. 1.2).Moreover, they easily generalize to problems involving multiple or curved sources, and can oftenbe implemented on data structures beyond triangulations (as shown in Fig. 7). On the other hand,since these methods compute only the geodesic distance ϕ , additional work must be done if onewishes to extract geodesic paths , as discussed in Sec. 5.PDE-based methods can be categorized according to the type of equation used to characterizegeodesic distance. Different starting points will lead to different numerical treatments, whichsubsequently have different computational trade offs ( e.g. , different mesh or solver requirements).At a high level there are two basic classes of methods: • Wavefront-based.
The basic principle behind wavefront-based methods resembles (in spirit)classical methods like Dijkstra’s algorithm, or the window-based methods discussed in Section4: information about distance propagates outward from a given source point. In the continuoussetting, this perspective corresponds to hyperbolic PDEs , i.e. , those that describe the evolutionof a wavefront emanating from the source. • Diffusion-based.
Diffusion-based methods more closely resemble problems arising in, say,spectral graph theory: information about distance is obtained by way of functions associatedwith a discrete
Laplace operator , computed via a process that looks more like repeated localaveraging rather than wavefront propagation. In the continuous setting, this perspectivecorresponds to elliptic and parabolic PDEs such as Poisson equations and heat diffusion.
Trade offs.
Historically, wavefront-based methods were developed prior to most diffusion-basedalgorithms; as a result, a wide variety of higher-order accurate strategies have been developed forregular grids on Euclidean domains ( e.g. , [1]), motivated in large part by accurate numerical solversfor level set equations [74]. However, most of these techniques do not immediately generalize tothe setting of curved surfaces, where one typically does not have a regular, uniform tessellationof the domain needed to support (for instance) larger finite difference stencils. Diffusion-basedalgorithms tend to have some nice performance advantages, since they are largely based on solvingeasy linear elliptic equations rather than more difficult nonlinear hyperbolic ones, though solutionscan sometimes be over-regularized ( e.g. , near the cut locus of the geodesic distance function). Aparticularly nice feature of linear methods is that one can often pre-factor the associated matrices,which in practice yields about an order of magnitude improvement in amortized performance usingmodern direct solvers [21, 89].Meshing requirements for elliptic methods also tend to be less stringent than for hyperbolicones; see Sec. 6 for further discussion. As discussed in Sec. 1.2, the accuracy of all PDE-basedmethods (as well as all methods from computational geometry) is fundamentally limited by theaccuracy of the domain representation itself ( e.g. , the use of linear elements to approximate curvedsurfaces). The only real resolution to this issue is to compute geodesic distance on higher-ordersurface representations such as splines and subdivision surfaces; several authors have recentlyconsidered this approach [30, 72].Overall, there is a rather nice viewpoint that unifies all current PDE-based methods: they are alleffectively trying to minimize the residual of the eikonal equation ∫ M (|∇| − ) dA , albeit by verydifferent means (see especially the discussion at the end of Sec. 3.3.3. Since this energy is nonlinearand nonconvex, it is not surprising to find that many different computational strategies have beenproposed: hyperbolic methods like fast marching deal with nonlinearity by taking advantage of , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 13
Fig. 7. Unlike methods from computational geometry, which are specially tailored to triangulations, PDE-based methods can easily be implemented on a variety of data structures by applying standard discretedifferential operators. Here, we show a variety of implementations of the heat method using different discreteLaplace operators. special update orderings, whereas elliptic methods like the heat method decompose the problem intotwo linear pieces connected by a nonlinear change of variables. Very little work has been done onsynthesizing these perspectives ( e.g. , combining fast linear approximations with carefully-orderedupdates), though is likely to be fruitful, especially given the diverse range of application problemsand machine architectures on which geodesic distance queries are needed.
PDE-based methods for geodesic distance computation all have a close relationship to the
Laplace-Beltrami operator ∆ , which is discretized by a (weighted) graph Laplacian . In particular, for a graph G = ( V , E ) with edge weights w : E → R , the graph Laplacian is encoded by a real symmetricmatrix L ∈ R | V |×| V | with off-diagonal entries L ij = (cid:40) − w ij , ij ∈ E , , otherwise , and diagonal entries L ii = (cid:213) ij ∈ E w ij . A simple choice of edge weights is just w ij =
1, though of course these weights do not capturemuch geometric information. When the vertices i ∈ V have associated vertex positions f i ∈ R , theedge weights can be determined by some function of the edge lengths; one common choice is touse Gaussian edge weights w ij = e −| f j − f i | / t (for some small parameter t > heat kernel (seeSec. 3.3). Finally when the edges of the graph come from the edges of a triangulated surface, one , Vol. 1, No. 1, Article . Publication date: July 2020. typically adopts the cotangent weights w ij = ( cot α ij + cot β ij ) , where α ij , β ij are the interior angles opposite edge ij (or zero, for edges on the domain boundary).More broadly, these discrete Laplacians are only the tip of the iceberg: discrete Laplace operatorshave been developed for a large variety of geometric data structures, providing a wide variety ofoptions for implementing PDE-based geodesic schemes (see especially the discussion in Sec. 3.3.3).Some methods, such as those based on diffusion (Sec. 3.3) easily generalize to these settings byjust “plugging in” a different Laplace matrix; other methods, such as fast marching (Sec. 3.2.1)do not immediately generalize since they may depend on particular features of a discretization( e.g. , the fact that values at the vertices of a simplex uniquely determine an interpolating affinefunction). Note that since the Laplace-Beltrami operator depends only on the intrinsic geometryof the domain ( i.e. , its metric) and not the way it sits in space, any distance algorithm defined interms of the Laplacian (including all the methods we consider here) will automatically be isometryinvariant. The basic prototype for wavefront-based methods is the linear hyperbolic wave equation d dt u t = ∆ u t , (1)where for each time t , u t is a real-valued function on a manifold M , and ∆ is the Laplace-Beltrami operator on M (for M = R n , ∆ is just the standard Laplacian). The intuitive connection to geodesicdistance is that a perturbation at some initial point ( e.g. , a small stone dropped in a pond) will sendout a wavefront that maintains a constant distance from the source point. Hence, the arrival timeof the wavefront is correlated with the geodesic distance. Though a few methods extract distanceinformation from the classical wave equation [92] or the quantum mechanical (Schrödinger) waveequation [43], by far the most common approach is to consider the nonlinear hyperbolic eikonalequation |∇ ϕ | = , on Mϕ = ∂ M (2)which directly characterizes the distance function ϕ : M → R in terms of the norm of its gradient |∇ ϕ | . Intuitively, this equation says something very straightforward: on the boundary of the domain ∂ M ( i.e. , at any point in the source set), the distance ϕ should be zero. At every other point of thedomain, the change in distance per unit distance along the direction of greatest increase shouldsimply be 1. However, actually solving this equation is not as straightforward since it is nonlinearand hence cannot be approximated using standard linear finite element methods. Instead, thegeneral strategy is to iteratively update the solution using techniques like (nonlinear) Gauss-Seidel;for carefully-crafted update orders (and suitable conditions on the mesh geometry) such strategiescan actually converge in a single iteration, corresponding to well-established algorithms such as fast marching (this perspective is nicely explained by Bornemann and Rasch [14]). The fast marching method was originally developed for distance transformson Euclidean domains discretized as regular grids; the basic principles of this method and its variants( e.g. , fast sweeping [65]) are shared by a broad class of schemes used in level set methods [75] andto solve Hamilton-Jacobi equations (of which the eikonal equation is one example). Kimmel andSethian developed a level set method for curved domains represented as triangulated surfaces [51].The basic strategy of this method is very similar to Dijkstra’s shortest path algorithm: set thedistance at the source point (or set) to zero, and set the tentative distance to infinity; then use a , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 15 n –1 n
42 3
Fig. 8. Methods based on solving hyperbolic equations may fail to respect causality in the presence of obtuseangles. For instance, in the left figure a source at 1 will get propagated to 2 and 3 before 4, even though 4is closer. In the worst case this phenomenon can be highly nonlocal, as seen on the right where the node n closest to the source node 1 is the last one to be updated. region-growing strategy to update the remaining distances in an “upwind” order, i.e. , consider thenode with smallest distance first. Like Dijkstra, the running time is therefore O (| V | log | V |) , sincefor a triangulation, | E | ∈ O (| V |) .Relative to Dijkstra, the key modification is that distances are not up-dated according to paths along edges; instead, one updates the distancevalue at a vertex by solving for the linear function that satisfies the eikonalequation [51, Section 4.1]. In particular, if the values ϕ i , ϕ j at two cornersare known, one simply needs to pick a third value ϕ k so that the slope ofa triangle |∇ ϕ | passing through all three values equals + |∇ ϕ | is the integrand of the Dirichlet energy , which in turn can be expressed in terms of the cotan Laplace operator . In particular, consider a triangle ijk where the distance values at vertices i and j are known, and the distance at vertex k remains to be determined. Suppose we encode thesedistance values as a column vector ϕ = [ ϕ i ϕ j ϕ k ] T , and let L ∈ R × be the local cotan matrix givenby 2 L ii = cot α j + cot α k , L ij = − cot α k (see for instance [26, Section 6.2]). A linear function ϕ interpolating the distance values at the threevertices will then satisfy the eikonal equation if ϕ T Lϕ =
1, which we can view as an ordinaryquadratic equation for the unknown distance ϕ k . Solving this equation yields (in general) twosolutions corresponding to positive and negative slope; the update will always use the larger value,since distance increases monotonically as we move away from the source.A basic difficulty with wavefront propagation methods is that the order in which vertices areupdated may violate causality , i.e. , even if a vertex i is closer to the source than a vertex j , thedistance at j may be finalized prior to the distance at i ; hence, the solution ϕ will fail to be monotonic(two examples of how this can happen are shown in Figure 8). Sethian and Kimmel suggest toresolve this issue by an “unfolding” procedure, but this procedure is nonlocal and may not terminatebefore reaching the boundary. In general, one therefore has to either apply an iterative strategy(as discussed by Bornemann and Rasch [14]), or perform remeshing, as discussed in Sec. 6. From acomputational point of view, the basic challenge with any method based on wavefront propagation(including Dijkstra’s algorithm) is that it is difficult to parallelize: the upwind strategy effectivelyinduces a serial ordering on computation. Since the order of computation is data-dependent ( i.e. , itdepends on earlier distance values), it can also lead to significant dynamic branching and incoherentmemory access compared to linear methods (Sec. 3.3) which more typically depend on a fixedtraversal order ( e.g. , using a fixed matrix factorization). , Vol. 1, No. 1, Article . Publication date: July 2020. The basic prototype for diffusion-based methods is the linear parabolic heat equation ddt u t = ∆ u t , (3)where for each time t the function u t is a real-valued function on the domain. Despite the verysimilar form of Eqn. 3 and Eqn. 1, they have very different behavior and hence lead to very differentcomputational methods. Whereas small perturbations to u yield high-frequency oscillations insolutions to the wave equation, such perturbations are diffused or smoothed out by the heatequation, hinting at greater stability ( e.g. , small computational errors are not propagated forwardin time). In particular, if δ x is a Dirac delta centered at a point x , then the heat kernel at x is thesolution to the heat equation ddt k t = ∆ k t , k = δ x . (4)More generally we will use k t ( x , y ) to denote the heat kernel at time t , centered at the point x , andevaluated at the point y . In the Euclidean case the heat kernel is just a Gaussian of increasing widthand decaying magnitude, though in general it has no closed form (see for instance [91, Section 3.2]).Any other solution to the heat equation can be expressed as a convolution with the heat kernel, i.e. , a gradual “blurring out” of the initial data. The heat equation also has an important statisticalinterpretation: if one views the initial function u as describing the spatial distribution of a largenumber of discrete particles, then the solution u t describes the distribution after these particleshave taken “random walks” ( i.e. , Brownian motion) for time t [57].Diffusion-based methods were initially used to compute smooth distance-like functions (Sec. 3.3.2)or functions that satisfy the axioms of a distance metrics, but do not actually provide the geometricdistance between points (Sec. 3.3.1); more recently, a diffusion-based strategy was introduced thatprovides the true geodesic distance (Sec. 3.3.3). Computationally, the appeal of all such methodsis they amount to sparse linear systems that can be efficiently solved using standard techniquesfrom numerical linear algebra. As a result, one does not have to develop solvers specially tailoredto the task of computing geodesic distance, as with the hyperbolic methods discussed in Sec. 3.2 orthe polyhedral strategies in Sec. 4. Instead, any development in numerical linear algebra (such asmore accurate linear solvers, or fast parallel implementations) can immediately be used to improvecomputation of geodesic distances. Diffusion-based methods have also become popular due to theirease of implementation, which generally involves only three basic ingredients: • Building a weighted graph Laplacian L associated with the mesh; • Solving one or more sparse linear algebra problems involving L ; • Evaluating simple per-triangle or per-vertex expressions.Per-element operations might be performed either before or after solving linear problems (or both),but generally amount to evaluations of, e.g. , simple closed-form sums (this pattern is in fact sharedby all methods we will discuss in this section). The bulk of the complexity is therefore taken care ofby the linear solver, which can be treated as a “black box”: one does not have to implement complextopological data structures, or even maintain a priority queue. On the flip side, this formulationassumes that a suitable mesh is already provided as input; see Sec. 6 for further discussion.
The original motivation for diffusion-based distances comes from a desireto embed an abstract surface or manifold into a Euclidean space – such an embedding allows oneto approximate distances on the original manifold by measuring the ordinary Euclidean distancebetween points in R n , though this approximation may only roughly resemble the true geodesicdistance. This point of view was originally studied by Bérard et al. for a different purpose: to study , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 17 the precompactness of smooth manifolds M with bounded diameter and curvature [9, Chapter VIE.53],[10]. Later, this same point of view became a common theme in the context of data analysis,machine learning, and dimensionality reduction, where the manifold M and Laplace-Beltramioperator ∆ are replaced with a discrete graph G = ( V , E ) and graph Laplacian L ∈ R | V |×| V | . Inparticular, let ψ i ∈ R | V | and λ i ∈ R be the eigenvectors and eigenvalues of L , so that Lψ i = λ i ψ i , i ∈ V . (5)Belkin & Niyogi propose an embedding by Laplacian eigenmaps that maps each vertex v ∈ V tothe coordinates f i : = ( ψ ( v ) , . . . , ψ k ( v )) (for some choice of n ≤ | V | ), arguing that this embeddingminimizes the ℓ distortion of edge lengths [7]. A notion of distance is then given by ϕ E ( i , j ) : = | f j − f i | , where | · | is just the Euclidean norm. Gobel & Jagers study a notion of distance that isclosely connected to the random walk interpretation of the heat equation – in particular, theyinterpret a graph as a Markov chain with probabilities determined by edge weights, and show thatthe expected time to walk from a vertex i to a vertex j and then back to i determines a metric on thisgraph [44]. This commute time distance ϕ C is equivalent to the effective resistance between i and j inan electrical network [52], which is also connected to the Laplacian L [41, 87]. Lipman et al. define aclosely related biharmonic distance ϕ B based on the bi-Laplacian ∆ rather than the Laplacian ∆ [59].Coifman & Lafon instead define a diffusion distance ϕ D which captures the amount of informationthat diffuses from a vertex i to a vertex j after time t [25]. In fact, all of these distances can be linkedback to a discrete diffusion process – consider in particular the (discrete) harmonic Green’s function G x : V → R > , which is the solution to the equation LG x = δ x (where δ x is the Kronecker deltacentered at x ∈ V ). This function can also be viewed as the stationary solution to the discrete heatequation with fixed (Dirichlet) boundary conditions at vertex x . Letting H x be the correspondingfunction for the bi-Laplacian, we can write the three distances mentioned above as ϕ C ( x , y ) = G x , x − G x , y + G y , y ϕ B ( x , y ) = H x , x − H x , y + H y , y ϕ D ( x , y ) = || k t , x − k t , y || , where k t , x : V → R > denotes the solution to the discrete heat equation at time t and with initialconditions u = δ x . Evaluating each of these distances for a given pair of vertices x , y ∈ V thereforeamounts to solving a small number of linear equations. If, however, one wishes to compute thedistance between a large number of vertex pairs ( e.g. , to solve the APGD problem), it is often moreefficient to use a spectral expansion of these distance functions, i.e. , to write them in terms of theeigenvalues and eigenvectors of the discrete Laplace operator (Equation 5), leading to a unifiedview of all four distances defined above: (eigenmap) ϕ E ( x , y ) = (cid:205) ni = ( ψ i ( x ) − ψ i ( y )) , (commute time) ϕ C ( x , y ) = (cid:205) ni = λ i ( ψ i ( x ) − ψ i ( y )) , (biharmonic) ϕ B ( x , y ) = (cid:205) ni = λ i ( ψ i ( x ) − ψ i ( y )) , (diffusion) ϕ D ( x , y ) = (cid:205) ni = e − tλ k ( ψ i ( x ) − ψ i ( y )) . Note that the global point signature of Rustamov [86] yields an identical distance to the commutetime distance ϕ C . At this point it becomes apparent that all of these distances are variations on acommon theme: embed the surface in R n according to the eigenfunctions; compute a weightedEuclidean distance in the embedding space. The different choice of weights will impact the regularityand other features of the resulting distance function—for instance, the diffusion distance exhibitsa degree of anisotropy around the source that depends on the diffusion time t , as shown in [59,Figure 2]. , Vol. 1, No. 1, Article . Publication date: July 2020. Notably, none of these distances give a good approximation of geodesic distance: on the wholethey look more like highly regularized versions of the true geodesic distance, and cannot resolvethe cut locus. However, they do provide an extremely efficient way to obtain a rough proxy fordistance: since all of these functions are quite smooth, one can obtain a good approximation byusing a truncated series where n is much smaller than the number of vertices | V | (in practice,around 150–200). In this case one can precompute the eigenfunctions, at which point it becomesextremely efficient to compute the point-to-point geodesic distance (by just evaluating a smallsum), and by extension, solving the APGD problem for a small subset of vertices ( e.g. , a collectionof landmark points) becomes relatively inexpensive. However, computing the distance to a large setof vertices (such as the boundary of the domain) may be expensive, since one needs to evaluate alarge sum for each point; here, other algorithms discussed in this section may be more appropriate.Finally, since these distances are ultimately derived from distances in Euclidean R n , they exactlysatisfy key properties of a metric: nonnegativity, symmetry, and triangle inequality. Note that any embedding into R n will immediately satisfy these properties; likewise, as with any algorithm builton top of the Laplace-Beltrami operator, these distances are also isometry invariant.Further connections between distance and spectral approximation are detailed in the recentsurvey by Patané [76]. There is a large class of methods for generating smoothed distance-like functions to the boundary ∂ M of a domain M by solving linear elliptic equations; these methodsare largely motivated by problems from image processing and computer vision, such as extractingmedial axes or other “skeletons” from image regions. For instance, Tari et al. [94] solve the screenedPoisson equation ( ∆ − ρ id ) u = M , u = ∂ M , (6)and Gorelick et al. [42] solve a Poisson equation ∆ u = − M , u = ∂ M ; (7)derivatives of such functions can then be used to extract an approximate cut locus/medial axis toserve as a “skelton” of the region M . Many variants of this strategy have been considered that yieldmore distance-like functions, such as applying various pointwise transformations or “normaliza-tions” to the solution u of PDEs like Equation 6 or 7; such as the Spalding-Tucker transformation u (cid:55)→ (cid:112) |∇ u | + u − |∇ u | [8, Section 6]. Much like the functions studied in Sec. 3.3.1, however, noneof these smoothed functions correspond to the actual geodesic distance ϕ . Moreover, unlike thefunctions from Sec. 3.3.1, they do not provide a well-defined distance metric d : M × M → R , butrather just compute the distance transform of the domain boundary ∂ M . The diffusion-based perspective can also be used to compute the true geodesicdistance, rather than a distance-like function, via the heat method [29]. This method is ultimatelyconnected to both the eikonal equation discussed in Sec. 3.2 as well as the Poisson-based approachesfrom Sec. 3.3.2, though the starting point is an important observation about the close relationshipbetween the geodesic distance function ϕ (Sec. 2.1.5) and the short-time heat kernel k t (Equation 4)known as Varadhan’s formula [96]: ϕ ( x , y ) = lim t → (cid:112) − t log k t ( x , y ) . (8)This formula effectively says that if a “pin prick” of heat centered at the point x diffuses for a very short time t , then the resulting heat distribution looks nearly identical to the geodesic distance , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 19 function, up to a simple transformation of the value at each point. However, it is quite challengingto compute a numerical approximation of the heat kernel that decays at exactly the right rate;Crane et al. sidestep this issue by recalling the eikonal equation (Equation 2), which says thatthe gradient must have unit magnitude ( |∇ ϕ | = u that is a monotonicfunction of distance can be used to obtain the distance itself, by normalizing the gradient andrecovering the corresponding scalar potential. In particular, to compute the distance to a point x one can apply the following procedure:(1) Solve the heat equation ddt u t = ∆ u t , u = δ x ;(2) Evaluate the vector field X = −∇ u /|∇ u | ;(3) Solve the Poisson equation ∆ ϕ = ∇ · X .Step 1 is approximated by a single step of backward Euler, i.e. , by solving the linear elliptic equation ( id − t ∆ ) u t = δ x (9)for a small time step t >
0. The Poisson equation in the third step corresponds to minimization ofthe energy ∫ M |∇ ϕ − X | dA , i.e. , it finds the function ϕ whose gradient is as close as possible to X (in the L sense). Note that all of these steps are described as operations in the smooth setting; aswith many PDE-based methods, a final algorithm is obtained by picking a particular discretizationof the domain and corresponding discrete differential operators. For instance, the heat methodhas been implemented on triangle meshes, polygonal meshes, point clouds [28], images [107],subdivision surfaces [31], tetrahedral meshes [8], spline surfaces [72], and voxelizations [18]. Avariant of the heat method that replaces the Laplace operator ∆ with the connection Laplacian ∆ ∇ also computes parallel transport of vectors along shortest geodesics [91], as well as the logarithmicmap discussed in Sec. 2.1. An interesting consequence of using higher-order elements of splinesand subdivision surfaces is that one can in principal obtain estimates of geodesic distance thatare even more accurate than the exact polyhedral distance—see [72] and in particular [30], whoexplore higher-order schemes.One can use the heat method to connect ideas from the various PDE-based methods consideredso far. From the perspective of Sec. 3.3.2, Equation 9 is a screened Poisson equation (where theboundary is just a single point ∂ M = { x } ), and Varadhan’s formula can be viewed as one possible“normalization.” However, as shown in Crane et al. [28, Figure 6], this simple transformation doesnot produce an accurate distance approximation, motivating the need for a more sophisticatednormalization strategy. The heat method can also be given a variational interpretation whichconnects it to fast marching methods [8]. In particular, consider the energy E ( u ) : = ∫ M (|∇ ϕ | − ) dA , i.e. , the L residual of the eikonal equation; the corresponding Euler-Lagrange equations are givenby the nonlinear equation ∆ ϕ = ∇ · (∇ ϕ /|∇ ϕ |) . The heat method can therefore be viewed as the first iteration of a fixed point scheme, where thesolution to the heat equation provides an initial guess ϕ for ϕ . Belyaev and Fayolle show thatimproved accuracy can be achieved by applying successive iterations ∆ ϕ k + ← ∇ · (∇ ϕ i /|∇ ϕ i |) ,and more broadly, by applying other optimization strategies to minimize E ( u ) , albeit at highercomputational cost. This point of view provides a clear connection between the heat method andthe fast marching method, which can likewise be viewed as a numerical method for minimizing , Vol. 1, No. 1, Article . Publication date: July 2020. E ( u ) (as discussed in Sec. 3.2.1). Finally, Litman & Bronstein note that the heat kernel k t is well-approximated via a spectral expansion akin to those seen in Sec. 3.3.1, enabling them to dramaticallyaccelerate the heat method via precomputation of Laplacian eigenvectors [60]. The methods reviewed in this section are aimed at resolving the boundary problems in the polyhedralsetting. We review both exact and approximated methods, by classifying them according to theapproach taken to the solution.
Global methods compute globally shortest paths. We distinguish between: methods that work onthe whole polyhedral domain and provide an exact solution (Section 4.1.1); approximations of suchmethods that trade-off accuracy for speed (Section 4.1.2); and methods that work on a discretizationof the domain on graphs, thus providing necessarily approximated solutions (Section 4.1.3).
Methods for solving the polyhedral geodesic distance problemare built on the piecewise flatness of polyhedral surfaces. This property enables the planar unfoldingof triangle strips, which simplifies the computation from 3D polyhedral surfaces to 2D unfoldedplanes (see Section 2.2). Within the unfolded triangle strips, locally shortest paths are computedas straight line segments. By enumerating all possible locally shortest paths between two points,globally shortest paths can be obtained by finding the ones with minimum length. However, withoutan efficient pruning strategy, the number of such locally shortest paths grows exponentially withthe size of Σ and quickly becomes infeasible to compute [6]. Therefore, the key challenge is toremove the computational redundancy as much as possible.Addressing this challenge, O’Rourke et al.[73] first proposed an algorithm to solve the PPGPproblem in O ( n ) time, where n is the number of vertices on Σ . Their method can be viewed asa theoretical milestone from that it shows the PPGP problem can be solved in polynomial time.However, their method is still too time-consuming for practical applications. Thus, our discussionstarts from [71], which laid the foundation of practical polyhedral geodesic algorithms. 𝑠 𝜎 𝑝 𝑑 𝑑 𝑎 𝑎 𝑤 𝜏 Fig. 9. Illustration of window structure. The shortest paths within an unfolded triangle strip can be encodedas a window w = ( a , a , d , d , σ , τ ) , where a , a are the endpoints of w ; d , d are the correspondingdistances from a , a to the pseudosource p (saddle vertex); σ denotes the geodesic distance from p to source s ; and τ denotes the propagation direction. Adapted from [83, 93]. Mitchell et al. [71] proposed the first practical algorithm for geodesic computation on polyhedralsurfaces, which is commonly referred to as the MMP (Mitchell-Mount-Papadimitriou) algorithm.Their main idea can be summarized as the continuous Dijkstra technique, which extends the famousDijkstra’s algorithm [35] from graphs to polyhedral surfaces. As an analogy, they propose to viewedges of polyhedral surfaces as nodes of a graph. However, an edge contains infinite points and its , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 21 distance cannot be represented by a scalar value. To handle this problem, a dedicated data structurenamed window is introduced, which encodes all locally shortest paths in an unfolded trianglestrip with a tuple (Figure 9). Then, the polyhedral geodesic distances can be computed by findingthe optimal windows on edges of Σ . The optimization of windows on an edge is accomplishedby trimming overlapping windows into disjoint ones according to the smaller distance in theoverlapping part (Figure 10). More specifically, windows are trimmed in a pairwise manner afterbeing propagated to edges and stored in an ordered list according to their positions. 𝑥 𝑠 𝑠 𝑏 𝑏
0𝑦 𝑥𝑠 𝑠 𝑝𝑑 𝑑(𝑎) (𝑏) Fig. 10. Illustration of window trimming. (a) Suppose there are two overlapping windows whose unfoldedpseudosources are s and s and their intersection δ = [ b , b ] . (b) By solving a quadratic equation, the twowindows are trimmed into disjoint ones. An illustration of special case σ = σ and then ∥ s − p ∥ = ∥ s − p ∥ = d is shown. From [93]. It can be observed that MMP’s time cost is positively correlated to the number of windowsarriving each edge. To minimize it, the MMP algorithm borrows the wavefront propagation paradigmappearing in Dijkstra’s algorithm and fast marching (Sec. 3.2.1) and propagates windows across Σ from near to far by maintaining a priority queue. Such paradigm ensures the redundant windowswill be trimmed at the earliest possible stage. When propagation is finished, each edge of Σ willbe subdivided into a list of end-to-end linked windows containing the geodesic distance field toany point on Σ . Mitchell et al. proved that the algorithm creates at most O ( n ) windows. Thus, itcan be easily derived that MMP can solve the SSGD problem in O ( n log n ) time and O ( n ) space,where n is the number of vertices of Σ . The key component of their proof is the analysis of themaximal number of windows arrived at each edge, which is O ( n ) . Nevertheless, such analysis istoo pessimistic and is inconsistent with MMP’s practical performance.The MMP algorithm [71] is commonly viewed as a landmark in the research of polyhedralgeodesic algorithms. Its distinct contribution is the window propagation framework , which containsthree major components: window propagation, window pruning (e.g. trimming) and windowmanagement (e.g. priority queue). This framework is employed by most of the following works[20, 84, 93, 104, 106] and they differ from their unique techniques used in the three components.Although it is straightforward for “continuous Dijkstra” technique to use a priority queue tomanage windows, an extra O ( log n ) time cost is introduced, which slows down the computation.Addressing this issue, Chen and Han [20] proposed the CH (Chen-Han) algorithm. Their algorithmmanages windows with a First-In-First-Out (FIFO) queue, whose overhead is O ( ) . For windowpruning, they proposed a very simple rule named “one angle, one split” to remove redundant In [71], Mitchell et al. used the term interval . While in the following works [84, 93, 104, 106], an alternative term window gains popularity from its intuitiveness. Thus, the term window is employed in this survey., Vol. 1, No. 1, Article . Publication date: July 2020. windows around vertices of Σ (Figure 11). This rule only cares about the redundancy around Improving Chen and Han’s Algorithm on the Discrete Geodesic Problem • s t d (a) A pseudo-source window s R e a,b [ ] (b) An interval window Fig. 2. A pseudo-source window and an interval window. w ws vv v I I Fig. 3. The principle of one angle one split (see Lemma 2.2) implies thatonly three of the four possible edge sequences in this example can pro-vide shortest paths. The three green windows belong to edge sequences thatwill continue, and the red one corresponds to an edge sequence that muststop because the shortest paths that it contains cannot continue to pass thistriangle. has at most n levels, where n is the total number of the faces. The CHalgorithm [Chen and Han 1990] avoids exponential explosion basedon the following Lemma 2.2, namely the property of “one angle onesplit” [Chen and Han 1990]. In Figure 3, two windows w , w , areentering (cid:3) v v v from edge v v and both cover the vertex v . “Oneangle, one split” implies that at most one of w , w , can have twointerval-window children during window propagation. We say thatthe two-child window occupies the angle (cid:4) v v v over the one-childone. Of course, w i ( i = ,
2) can also have a pseudo-source-windowchild at vertex v only if it can provide the shortest distance to v .L EMMA
Let w and w be two windows on the same directededge v v of (cid:3) v v v , shown in Figure . If both of the windows coverthe vertex v , then at most one of them can have two children whichcould be used to define a shortest sequence.
3. THE IMPROVED CH ALGORITHM
The CH algorithm [Chen and Han 1990] is of an O ( n ) timecomplexity, but in practice it runs very inefficiently, comparedwith the O ( n log n )-time MMP algorithm [Mitchell et al. 1987].This is demonstrated by previous experimental results [Kaneva andO’Rourke 2000; Surazhsky et al. 2005]. The inefficiency of the CHalgorithm, we think, is because it generates considerably many use-less windows that don’t contribute to any shortest path. That is tosay, the key idea in Chen and Han [1990], “one angle one split,”is still too loose to effectively check the validity of a new win-dow. We will describe the CH algorithm in Section 3.1. Section3.2 presents a filtering theorem for checking a new window morestrictly and Section 3.3 suggests that we should also maintain a pri-ority queue throughout the algorithm. Finally, in Section 3.4 we propose a method for backtracing a shortest path from the source toany vertex, or even any surface point. We first need to conduct initialization: we set the distances at ver-tices, except the source, to be +∞ ; then, we associate each anglewith a null interval window and associate each vertex with a nullpseudo-source window; finally, we introduce a first in first out queue Q to store pending events.A LGORITHM
The CH algorithm
Assign the source s with distance 0, create a pseudo-source win-dow w for s , and put w into the queue Q ; While Q is not empty and the level size doesn’t exceed the facenumber n Take out the head window w from Q ; If w is a pseudo-source window, say, w = ( d , v ) If d is less than the current distance estimate at vertex v Update the distance at v ; If v is a saddle vertexDelete the old pseudo-source window at v and its sub-trees;For each edge opposite to v , add a child window( d , v , e , [0 , Q ;Update the distance of each vertex v (cid:6) incident to v with w and add a pseudo-source window ( d + (cid:7) v v (cid:6) (cid:7) , v (cid:6) ) to Q if d + (cid:7) v v (cid:6) (cid:7) is less than the current distance at v (cid:6) ; Else /* w is an interval window, say, w : = ( d , I , e , [ a , b ]).*/ If w has only one child on the left (right) edge, or w failsto occupy the opposite angle over the existing window w (cid:6) according to Lemma 2.2, thenCompute the only child and push it into Q ; Else /* w occupies the opposite angle over w (cid:6) */Delete the abolished subtree of w (cid:6) ;Compute the two children of w and push them into Q ;Check if w can provide a shorter distance to the vertex v opposite to edge e ; if true, update the distance estimate at v ; and if v is a saddle vertex or a boundary vertex, we needalso generate a pseudo-source window at v and insert itinto the priority queue Q .Chen and Han [1990] proved that the total number of nodesgenerated—the abolished nodes included—is O ( n ), rather than ofan exponential size. But in order to achieve an O ( n ) space complex-ity, they suggested only storing leaf nodes and branch nodes, whilethrowing away the one-child interval windows. So the sequence treeof the CH algorithm is a conceptual object rather than a real one.When a pseudo-source window w v at vertex v is occupied over by ACM Transactions on Graphics, Vol. 28, No. 4, Article 104, Publication date: August 2009.
Fig. 11. Vertex v splits any window spanning it into two halves. The “one angle, one split” rule shows thatfor any two windows passing through v , only three of the four split windows can provide shortest paths andcontinue propagation. As illustrated, the three green windows are valid and continue propagation. The redone contains no shortest paths beyond triangle ∆ vv v and stops propagation. From [104]. vertices and is tolerant of overlapping windows. Thus, it is much less powerful than the trimmingrule used by the MMP algorithm [71]. However, it can still be proved that the number of windowsgenerated in CH algorithm is at most O ( n ) . This analysis also supports that the worst-case analysisof the MMP algorithm is too pessimistic and that its practical performance should be much better.Since the CH algorithm uses a O ( ) FIFO queue to manage windows and has a O ( n ) windowcomplexity, its overall time complexity is O ( n ) . In addition, the CH algorithm does not storepropagated windows on edges to perform window pruning. Thus, its space cost is O ( n ) .Theoretically, the CH algorithm achieves the best asymptotic complexities so far. However, itspractical performance is poor. Surazhsky et al. [93] reported that their implementation of theMMP algorithm runs in sub-quadratic time and is many times faster than Kaneva and O’Rourke’simplementation of the CH algorithm [50]. Since then, several methods have been proposed toimprove the practical performance of the MMP algorithm. Liu et al. [64] observed that floatingpoint error may cause the degeneration of window propagations to frequently occur when applyingthe MMP algorithm to real world models. To make the MMP algorithm more robust, they con-ducted a systematic analysis on all the degenerated cases and proposed techniques to handle themaccordingly. Observing that the half-edge data structure used in Surazhsky et al.’s implementation[93] may generate redundant windows, Liu [62] proposed to implement the MMP algorithm withan edge-based data structure. Experimental results show that on average the edge-based versionruns 44% faster and uses 29% less memory.In summary, the MMP algorithm’s practical success comes from its effective wavefront propaga-tion paradigm which enables the removal of redundant windows at the earliest stage. However,the removal of redundant windows is a complex and expensive process which involves inserting anewly propagated window into a ordered list and trimming overlapping parts by solving quadraticequations.Addressing this issue, Xin and Wang [104] proposed the ICH (Improved Chen-Han) algorithm,which combines the advantages of both the MMP algorithm and the CH algorithm. From the MMPalgorithm, they borrowed the wavefront propagation paradigm and used a priority queue to manage , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 23 window propagations according to their distances. To avoid the MMP’s costly window trimmingoperations, they borrowed the “one angle, one split” rule from the CH algorithm and extended itinto three novel window filtering rules (Figure 12). Experimental results show that these new rules
Fig. 12. Illustration of ICH’s window filtering rules. The key idea is to filter out redundant windows with theminimum-so-far distance d , d , d of vertices. Left: d + ∥ IB ∥ > d + ∥ v B ∥ and its symmetric counterpart d + ∥ IA ∥ > d + ∥ v A ∥ ; Right: d + ∥ IA ∥ > d + ∥ v A ∥ . From [104]. help to remove more than 99% redundant windows during propagation. Theoretically, introducingthe priority queue increases ICH’s time complexity to O ( n log n ) and makes its space complexityno longer O ( n ) . However, its practical performance increased dramatically. Experimental resultsshow that ICH greatly outperforms the original CH algorithm. In addition, although ICH usuallypropagates more windows, it runs comparable to the MMP algorithm [93] while using considerablyless space. This reveals that there exists a trade-off between the effectiveness of window pruningstrategies and their costs.One prominent difference between MMP and ICH’s window pruning rules is that: MMP’strimming rule involves interaction between windows while ICH’s “one angle, one split” and windowfiltering rules do not. More specifically, ICH filters redundant windows by distance comparisonwith the minimum-so-far distances at vertices, which is a relatively independent process and isideal for parallelization. However, there is still one obstacle. To remove redundant windows at theearliest stage, ICH organizes window propagations in a strict order with a priority queue, which issequential. Fortunately, since the correctness of window propagation algorithms is independentof the order of window propagations, this strategy can be loosened for parallelization. Based onthe above observations, Ying et al. [109] proposed the PCH (Parallel Chen-Han) algorithm, whichis a parallel version of the ICH algorithm. In their algorithm, k nearest windows are selected andpropagated in parallel at each iteration, where k is a user-specified parameter (Figure 13). Then, thepropagated windows are filtered in parallel according to the minimum-so-far distances at vertices.To avoid data conflict at vertices, the distance updates are delayed until the propagation of selectedwindows finished. Experiment results show that PCH propagates slightly more windows than theICH algorithm, but runs an order of magnitude faster.The performance of the PCH algorithm shows that slightly relaxing the order of propagation willnot have a big impact on the window pruning effectiveness. Based on this spirit, [106] proposedthe FWP (Fast Wavefront Propagation) technique to accelerate the MMP and ICH algorithm byreplacing the strictly ordered priority queue with a loosely-ordered bucket-based FIFO queue. Morespecifically, windows with similar distances are stored in the same bucket and propagated in aFIFO order. Although straightforward, their method performs well in practice from the observationthat both MMP and ICH spend roughly 70% time on maintaining the priority queue (Figure 14).Experimental results show that their FWP-MMP algorithms runs 3-10 times faster than the MMPalgorithm, and their FWP-CH algorithm runs 2-5 times faster than the ICH algorithm. , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 13. The PCH algorithm selects k nearest windows and assigns them to different GPU threads. Eachthread propagates the selected windows and filters the redundant ones independently. Then, the newlycreated windows are collected and re-organized in the memory pool. From [109]. Both PCH and FWP employ a batch processing strategy to accelerate existing algorithms. How-ever, their selection of window batches depends only on the distances and lacks the considerationof their geometric interrelationship. Observing that the prunings of all propagated windows withina triangle are inter-dependent, Qin et al. [84] proposed the VTP (Vertex-oriented Triangle Propaga-tion) algorithm. Their algorithm organizes window propagations with a triangle-oriented growingscheme. In this scheme, a traversed area is defined as the union of all visited triangles. The boundaryof this traversed area forms the propagation wavefront, which contains all the windows to bepropagated. Then, this traversed area is expanded in a Dijkstra-like style by gradually enclosingadjacent triangles (Figure 15). Their algorithm terminates when the traversed area covers the entire Σ . During expansion, local windows entering the same triangle from the same edge are organizedin a batch and they are propagated simultaneously. In such batches, redundancy checks can beintensively performed between any pair of windows. To remove a maximal number of windows in alow cost way, they proposed three rules which are summarized from an exhaustive list of scenariosfor pairwise window pruning inside a triangle. Note that these pruning scenarios are listed underthe assumption that the window trimming technique [93] is not used due to its relatively high , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 25 speaking, the smaller the variance, the higher quality thewavefront has. Wavefront quality depends on the order ofwindows being processed. As shown in Fig. 4, propagatingwindows in the smallest-label-first order leads to a high-quality wavefront (i.e., smooth and with small variance),whereas using the first-in-first-out order produces a low-quality wavefront (i.e., very rough and with large variance).The MMP and ICH algorithms always take the windowwith the minimal label at each iteration, resulting in a high-quality wavefront. However, the overhead required for eachiteration is expensive. Computational results show that morethan 60 percent of the time is used for maintaining a priorityqueue in the MMP and ICH algorithms. See Fig. 5. On theother hand, the CH algorithm organizes the windows in anFIFO queue, leading to a constant time overhead. However,as its wavefront is of poor quality, the CH algorithm produ-ces many useless windows and converges very slowly.Our idea is to balance the wavefront quality and theoverhead for updating wavefronts. Unlike the MMP/ICHalgorithms that propagate only the window with the small-est label in each iteration, our method propagates at least K smallest-label windows at the same time. We propose abucket data structure to organize windows so that it takesonly O ð Þ time to process each window. In addition, K isadaptive to the wavefront size, and it has a constant upperbound, leading to an O ð n Þ worst-case window complexity.We call our method fast wavefront propagation , which distin-guishes itself from the existing slow-propagating algorithmssuch as MMP and ICH. Let s V be the source vertex and p M a point (notnecessarily a mesh vertex) on M . Denote by d ð p Þ the geo-desic distance between s and p . Obviously, d ð is a con-tinuous function. We partition the polyhedral distanceinto equal-length intervals, ½ ; l Þ , ½ l; l Þ , ½ l; l Þ , etc. Eachinterval is called a bucket , which is used to organize win-dows. Observe that the maximum range of geodesic dis-tances in most real-world models is O ð ffiffiffi n p h Þ , where h is the average edge length. There are O ð n Þ windows andeach bucket contains roughly O ð n Þ windows, so we setthe bucket size l ¼ h= ffiffiffi n p . Algorithm 1.
Fast Wavefront Propagation Algorithm
Input: a mesh M ¼ ð V; E; F Þ and a source point s ; Output: the undiscretized geodesic distance for each vertex;1: for each edge e facing s do
2: generate a window w for e with w .birth ¼ & ;3: insert_window( w );4: end for K ¼ P ¼ c small ¼ iter ¼ , c large ¼ wavefront_size;6: while wavefront_size ¼ do K + ¼ ð c large & c small Þ ;8: K ¼ min ð max ð K; Þ , wavefront_size, K max );9: // find the first non-empty bucket10: while buckets[ P ].is_empty() do P ++;12: end while
13: // find P end so that at least K smallest-label14: windows will be propagated15: P end ¼ P , wc ¼ ;16: while wc < K do wc + ¼ buckets[ P end ].Q.size();18: P end ++;19: end while P ¼ P ;21: while P ’ P end do
22: //propagate the windows born in early iterations23: while buckets[ P ].top().birth < iter do w ¼ extract_window( P );25: propagate w across its adjacent triangle;26: for each child window b w do if b w .dist() ( P end ) l then
28: // b w is in a bucket after P end c large ++;30: else c small ++;32: end if b w .birth ¼ iter ;34: insert_window( b w );35: end for end while P ++;38: end while iter ++;40: end while Fig. 5. Maintaining a priority queue (pq) takes roughly 70% of the runtimeof the MMP algorithm. The ICH algorithm has a similar percentage.Fig. 4. The MMP algorithm maintains a priority queue for the windows,leading to a high-quality wavefront with standard deviation std ¼ : . Replacing the priority queue by an FIFO queue results in a verypoor and slow-moving wavefront with std ¼ : . The model hasbeen scaled to a unit cube.
1. In a few extreme pathological cases, the maximum geodesic rangecan reach O ð nh Þ . However, our bucket strategy still works well in prac-tice with l ¼ h= ffiffiffi n p . XU ET AL.: FAST WAVEFRONT PROPAGATION (FWP) FOR COMPUTING EXACT GEODESIC DISTANCES ON MESHES 825
Fig. 14. The running time breakdown of the MMP algorithm. Note that the maintenance of priority queue isthe most time-consuming part which occupies around 70% of the running time. The ICH algorithm has asimilar percentage. From [106]. S (c) (d) (a) (b) S A B C D EF A B C D E S A B C DE A R R (cid:39) Figure 6:
Triangle-oriented wavefront propagation over a mesh. (a)-(b) Face-sorted wavefront propagation, (c)-(d) Vertex-sorted wavefrontpropagation.
Step 2 . If j == N U LL , finish; otherwise, perform pairwise win-dow pruning between w i and w j using Case 1, Case 2 and Case 3in Section 3. Step 3 . If w j is removed from the list in Step 2, set j = j + 1 and goto Step 2. In the event that w i is removed in Step 2, if i == wl.head , set i = j, j = j + 1 and goto Step 4; otherwise,set i = i and goto Step 2. If neither w i nor w j is removed, set i = j , j = j + 1 and goto Step 4. Step 4 . If j == N U LL , finish; otherwise, goto Step 2.There is a double loop in the above procedure. Index j is associat-ed with the outer loop and index i is associated with the inner loop.This procedure is illustrated in Fig 5(c), where it traverses all win-dows in the outer loop (red arrow) and checks each window againstits preceding windows in the inner loop (black arrow). Its timecomplexity is O ( m ) which will be discussed in the next section. Checking with vertices
During the process of enforcing Rule 2,for each propagated window on AC , we also apply Case 12 in Fig3 by checking the window against the distance to vertices (the sameas the filtering rule in ICH [Xin and Wang 2009]). Suppose we already have a window list wl l = { w l , w l , ..., w lm } on AC , which is propagated from AB . In this subsection, we presentthe following procedure to propagate windows from another list f wl r = { w , w , ..., w n } from BC to AC , and merge the propa-gated windows with wl l . Meanwhile, we perform window pruningusing Cases 1-6 in Fig 3.Procedure P rimeM erge ( wl l , f wl r ) ( Rule 3 ) consists of the fol-lowing steps. First, perform one step of propagation for all win-dows in f wl r . Let w i be the propagated version of w i . Then, foreach window from w to w n , run the following substeps: (i) ap-pend it to wl l ; (ii) set j = wl l .tail and i = j ; (ii) performpairwise checking and window pruning on the updated wl l usingSteps 2-4 in Rule 2 except that in Step 2, instead of consideringCases 1-3 only, we need to check where the two windows are fromand use either Cases 1-3 (if both windows are propagated from thesame edge) or Cases 4-6 (if the two windows are propagated fromtwo different edges).We name this procedure P rimeM erge () because it will be usedfor merging window lists on an edge for the first time. It is comple-mentary to SecondM erge () in Section 5.1. The time complexityof P rimeM erge () is O ( m + n ) , which will be discussed in thenext section. Fig 5(d) shows an illustration of the outer loop (redarrow) and the inner loop (black arrow) of this procedure. Order Preservation.
A window list wl = { w , w , ..., w k } isspatially coherent if w i .a w i +1 .a for all i = 0 , ..., k . Proposition 4.2.
If both wl AB and wl BC are spatially coherent,the window list wl l = wl AB ! AC obtained after applying Rule 1and Rule 2 is also spatially coherent. And the merged list obtainedafter applying Rule 3 is still spatially coherent. The proof of this Proposition is given in Appendix C.
Our geodesic algorithm takes triangles as the primitive for synchro-nizing window and distance propagation. All visited triangles forma single connected region, called the traversed area, over the meshsurface. The boundary of this traversed area is defined as the prop-agation wavefront. Our algorithm expands this traversed area in acontinuous Dijkstra style by gradually enclosing unvisited trianglesabutting the traversed area. During each iteration, our algorithmadds one or more unvisited triangles to the traversed area, and thewavefront is also updated. Let R and R be the existing and ex-panded traversed area respectively. We denote the region outside R but inside R as R , which consists of the newly added trian-gles. In this section, we first present a basic face-sorted propagationalgorithm, and then extend it to vertex-sorted propagation, whichachieves improved performance. As shown in Fig 6(a) and (b), our face-sorted geodesic algorithmexpands the traversed area one triangle face at a time. Its outline isgiven below.
Initialization.
Create a single window for every opposite edge of S in its 1-ring neighborhood (bold blue lines around S in Fig 6(a)),and push all triangles that are outside the 1-ring neighborhood of S and that share at least one opposite edge of S to Q . Set D ( S ) = 0 , D ( P ) = dist ( S, P ) if P is a 1-ring neighbor of S , and D ( V ) = for all other vertices. Wavefront Propagation.
Step 1 . Pop the triangle with the highest priority from Q and add itto R . This single triangle forms R . Step 2 . If this triangle has only one edge on the previous wave-front ( DEF in Fig 6(a)), propagate the window list on DE toboth DF and F E using Rule 1 and Rule 2. Push adjacent trianglessharing either DF or F E with DEF into Q and calculate theirpriority; if any of these adjacent triangles is already in Q , simplyupdate its priority. Step 3 . If the popped triangle has two edges on the previous wave-
ACM Trans. Graph., Vol. 35, No. 4, Article 125, Publication Date: July 2016
Fig. 15. Illustration of triangle-oriented growing scheme. The black and red line segments denote theboundary of the traversed area. The green line segments denote the triangle to be added into the traversedarea. Left: The triangles in vertex A ’s 1-ring neighbourhood is adding into the traversed area. The red linesegments become interior edges. Thus, the windows on them are propagated through the traversed area untilthey reach the updated boundary edges. Right: R is the traversed area before growing. After growing, ∆ R willbe merged into R . From [84]. computational cost of solving quadratic equations. Experimental results show that their algorithmruns 4-15 times faster than MMP and ICH algorithms, 2-4 times faster than FWP-MMP and FWP-CHalgorithms, and also consumes the least memory, which ranks it as the best performing polyhedralgeodesic algorithm to date. Remark:
Most polyhedral geodesic algorithms aim at solving the SSGD problem, while Bala-subramanian et al. [6] proposed an algorithm to solve the APGD problem. The main idea of theiralgorithm is to build a vertex-to-vertex graph by computing the minimal-geodesic distances betweenall pairs of vertices on Σ . The minimal-geodesic distance is the length of the shortest path betweentwo vertices under the condition that the path contains no intermediate vertices. Then, the shortestdistances between any pair of vertex can be computed by searching the vertex-to-vertex graph using , Vol. 1, No. 1, Article . Publication date: July 2020. standard algorithms like the Dijkstra’s algorithm [35]. To compute the minimal-geodesic distance,they employed a triangle chain flattening method and proposed to reduce the redundancy throughvisibility. However, the overall method is essentially the same as the standard window propagationalgorithms without any pruning. Thus, its practical performance is expected to be much worsethan other state-of-the-art geodesic algorithms. Note that in [6], the runtime of Surazhsky et al.’s[93] MMP implementation is estimated but not tested by experiments. Thus, the comparison maybe inaccurate. Shortest Path Construction.
All the reviewed methods employed a straightforward backtracing strategy to construct geodesic paths, thus supporting SSSP, too. In general, there are two relatedapproaches to answer two queries respectively, which can be easily adapted to all polyhedralmethods: • Shortest paths to vertices:
As Figure 16 shows, for a vertex v ∈ Σ whose shortest distanceis obtained from window w , let I be the image of w ’s pseudosource on the plane defined by v and w . Then, the entering point p on edge e can be computed easily by intersecting vI andedge e . The entire path can be constructed by backtracing in a similar way. Note that theimage I for each vertex of Σ can be recorded during window propagation and do not requireany extra data structure, which is space-efficient. • Shortest paths to generic points:
As Figure 17 shows, propagated windows must be storedon edges of Σ as the data structure supporting geodesic path queries to points. For aninterior point p of a triangle, the window containing its shortest path can be computed by w p = arg min w ∈ W ∥ p − p ′ ∥ + D ( p ′ ) , where W is a set containing all windows on the threeedges of the triangle, p ′ is a point in w , D ( p ′ ) is the geodesic distance of p ′ . Then, the shortestpath of p can be computed by backtracing from p to p ′ in w p until reaching the source. Notethat this approach can be space-consuming since it requires storing all propagated windowson edges. Improving Chen and Han’s Algorithm on the Discrete Geodesic Problem • s I e[a,b]p Fig. 5. (cid:7) w (cid:7) = d + min p ∈ [ a , b ] (cid:7) I p (cid:7) . The CH algorithm performs window propagation level by level,rather than from near to far. Since the “by distance” manner is morecompatible with the shortest path problem, maintaining a priorityqueue hopefully further optimizes the CH algorithm. Experimentalresults in Section 4 suggest maintaining a priority queue does furtherimprove the CH algorithm greatly in both time and space. In spiteof this, it must be noted that the theoretical asymptotic runningtime of the improved algorithm becomes O ( n log n ) and the spacecomplexity cannot be proved to be still O ( n ). At the conclusion of the CH algorithm, the information encoded onvertices and angles is enough for backtracing the shortest paths tovertices. Consider a vertex v , whose shortest distance may be givenby either an interval window occupying one of v ’s incident anglesor a pseudo-source window occupying one of v ’s incident vertices,say, v (cid:6) . If the neighbor vertex v (cid:6) provides the shortest distance, thenthe shortest path to v consists of two parts: the shortest path to v (cid:6) and segment v (cid:6) v . Otherwise, it is an interval window w on edge e that provides the shortest distance. According to the source im-age I encoded in window w and the position of vertex v , we caneasily compute the entering point p on edge e and the entering direc-tion, which enables us to continue backtracing the shortest path, asFigure 6 shows. Each shortest path will end with the source vertex s .However, with some slight modification the information producedby our algorithm can be used to find the shortest path to a point inthe interior of a triangle, using an idea similar to that proposedby Surazhsky et al. [2005]. As is shown in Figure 7(a), for point p ∈ f , they considered all windows on the three edges boundingthe face f , and minimize (cid:7) p − p (cid:6) (cid:7) + D ( p (cid:6) ) over all points p (cid:6) withinthese windows. After we choose the best window, the backtracingprocess can continue in the same way. However, Chen and Han[1990] suggested throwing away those one-child interval windowsto ensure the O ( n ) space complexity, making it not informative tocompute the shortest path to a point in a face interior if we onlyimplement the first phase of the CH algorithm. So we suggest thatthe one-child interval windows (useless windows excluded) shouldbe kept on the corresponding edges. Thus for a destination point p ∈ f , we can also backtrace the shortest path in the same way.Although we get different information at the conclusion of theMMP algorithm and the improved CH algorithm, as Figure 7 shows,both of the algorithms can be used to backtrace shortest paths to sur-face points. In fact, for the improved CH algorithm, it is possible thatsome directed edges are only partly covered by interval windows,but the information is still enough to backtrace the shortest path toany surface point. It is at the expense of more space. The correctness I e vp
Fig. 6. Backtracing shortest paths to vertices. is based on two observations: (1) at the conclusion of the original CHalgorithm, each directed edge is covered by a list of interval windowsif we don’t abolish the one-child interval windows; furthermore,these windows encode all shortest edge sequences; and (2) Theo-rem 3.2 does nothing but filter out those totally useless windows.As regards the complexity of backtracing a shortest path, whetherused with the MMP algorithm or the improved CH algorithm,traversing the windows of a triangle might require O ( n ) time tofind the distance from a destination point p to the source s . This isbecause a triangle, in either algorithm, might have O ( n ) windowson its boundary. (Surazhsky et al. [2005] report that most trianglesseem to have O( √ n ) windows on average.) A complete implementa-tion of either the CH algorithm or the MMP algorithm would give adata structure that could report the distance from an arbitrary surfacepoint p to the source s in O (log n ) time. But a complete implemen-tation would probably be impractical and unnecessary. (Here wemust thank an anonymous referee who gave this accurate and acuteobservation.)
4. EXPERIMENTAL RESULTS
In this section, we provide some experimental results. Our experi-ments are made on an HP Compaq dc7800 computer with the fol-lowing configuration:—Intel(R) Core(TM)2 Duo CPU;—E8400 @3.00GHz;—2.99GHz, 2.98GB of RAM;—Microsoft Windows XP Professional SP3.In this research, we compare exact algorithms including(1) the MMP algorithm [Mitchell et al. 1987] implemented byDanil Kirsanov, one of the authors of Surazhsky et al. [2005]:Note that Kirsanov’s code was not used for that paper;(2) the original CH algorithm [Chen and Han 1990];(3) the ICH1 algorithm: with the filtering theorem (Theorem 3.2)to improve the CH algorithm;(4) the ICH2 algorithm: further maintaining a priority queue. ACM Transactions on Graphics, Vol. 28, No. 4, Article 104, Publication date: August 2009.
Fig. 16. Illustration of backtracing from a vertex of Σ . From [104]. Methods for exactly solving polyhedral geodesic prob-lems can be time-consuming for large-scale applications. Thus, some methods approximate thepolyhedral distance, which can be applied in scenarios insensitive to accuracy.To solve the SSGD problem efficiently, Surazhsky et al. [93] proposed an approximated versionof the MMP algorithm to reduce the time and memory costs. This algorithm works just as the MMPalgorithm. The only difference is that: in the approximated version, the algorithm tries to merge a , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 27 • S.-Q. Xin and G.-J. Wang w w w w w w w w w w w w w w w w p I s (a) At the conclusion of the MMP algorithm, eachdirected edge is covered by a list of disjoint intervalwindows. N o w i n do w e n t e r i ng f r o m he r e w w w w w w w w w w w w w p I s (b) At the conclusion of the improved CH algorithm, it ispossible that some edges are only partly covered by intervalwindows. Fig. 7. Backtracing shortest paths to a surface point.Fig. 8. Test models.
ACM Transactions on Graphics, Vol. 28, No. 4, Article 104, Publication date: August 2009.
Fig. 17. Illustration of backtracing from a point of Σ . In this scenario, propagated windows are required to bestored on edges. However, they can be organized in different ways. Left: windows on edges are trimmed intonon-overlapping ones as the MMP algorithm do [93]. Right: redundant windows on edges are filtered as theICH, VTP algorithms do [84, 104]. The remaining windows may overlap with each other. From [104]. xy q (cid:1) q (cid:1) s (cid:1) s s xy q (cid:1) q (cid:1) s (cid:1) s s γγ V ( b ) Figure 8: (a) Merging two windows into one such that the visibilityregion is not reduced. (b) Selecting the pseudosource for the newwindow. The visibility constraint V corresponds to the yellow regionwhile the constraints of (4.3) define the pink region. ∆ D ( p ) (actually, the relative local error ∆ D ( p ) / D ( p ) ) by a frac-tion of the global error tolerance. We set this fraction to be inall our experiments, although ideally this value should depend onthe size of the mesh. Thanks to this heuristic, we are able to satisfythe global error bound using significantly fewer windows overall. Monotonicity:
We must be careful that the distance function overthe edge is not smaller than the distance function over the parentwindows from which it was propagated, where correspondence isdefined by the propagation method. However, it is difficult to ana-lytically check this property in practice. When we do not check this,loops in window propagation are possible. Because we maintainconsistency of the source direction τ , our algorithm cannot produceany “bouncing” infinite loops, propagating back and forth betweentwo adjacent edges. It is conceivable, though, that more exotic in-finite loops could occur. Of course, we could explicitly maintaina directed graph representing the propagation evolution, explicitlycheck for loops in this graph, and disallow any attempted propaga-tion steps that would complete a loop. In practice, we have neverencountered such loops in any of our experiments, and so we do notexplicitly perform these checks. We attempt to find a pseudosource position for the merged windowthat satisfies all the constraints of the previous section.Denote the two adjacent windows w i , with pseudosources s i , for i = 0 , . Using the same local coordinate frame used previously,the merged window w (cid:1) will have two endpoints q (cid:1) i ≡ ( b (cid:1) i , . Ourgoal is to find a new pseudosource s (cid:1) = ( s (cid:1) x , s (cid:1) y ) and pseudosourcedistance σ (cid:1) for the merged window that satisfy the following con-straints.To maintain continuity, we require that the geodesic distances at itsendpoints D i ≡ D ( q (cid:1) i ) = (cid:3) s i − q (cid:1) i (cid:3) + σ i are preserved by themerge. This can be expressed as ( s (cid:1) x − b (cid:1) i ) + s (cid:1) y = ( D i − σ (cid:1) ) . It follows that σ (cid:1) = αs (cid:1) x + β ; (4.1) s (cid:1) y = As (cid:1) x + Bs (cid:1) x + C, (4.2)with A = α − , B = 2 α ( β − D ) + 2 b (cid:1) , C = ( D − β ) − b (cid:1) , α = ( b (cid:1) − b (cid:1) ) / ( D − D ) , β = ( b (cid:1) − b (cid:1) − D + D ) / (2( D − D )) .This constrains s (cid:1) to lie on a conic curve γ .To maintain directionality, we must impose the inequality s (cid:1) y ≥ .To satisfy the visibility constraint we require our solution to lie in-side the sector between the two lines L i , i = 0 , that pass throughthe q (cid:1) i to s i (Figure 8( b )). (If the intersection point of L i has posi-tive y -coordinate, then the allowed region V is a triangle. Otherwise V is open.) To constrain σ (cid:1) and the d (cid:1) i to be non-negative we add the inequali-ties σ (cid:1) ≥ , (cid:3) q (cid:1) i − s (cid:1) (cid:3) = D i − σ (cid:1) ≥ . It follows from (4.1) thatthese inequalities are equivalent to: (cid:2) − βα ≤ s (cid:1) x ≤ D − βα , if α > , D − βα ≤ s (cid:1) x ≤ − βα , if α < . (4.3)If all the above constraints are not simultaneously satisfiable we dis-allow the merge. Otherwise we pick the s (cid:1) with minimal σ (cid:1) value.This must occur when one of our inequality constraints is “tight”. The geodesic path for our approximation algorithm is traced sim-ilarly to the algorithm in Section 3.5 but with one essential differ-ence. When window w is the result of merging two original win-dows, its pseudosource position is different from the pseudosourcepositions of those original windows. If we were to trace back in thedirection of the merged pseudosource, the resulting path would bedifferent from that represented in the forward propagation, and itsoverall length might exceed the computed error bound.Our approach is to obtain the original pre-merge pseudosource bymaintaining a list of references to the windows that were succes-sively merged into w (together with their endpoints). The averagelength of these lists is only about 2 in all our experiments. Anotherbenefit of using the correct pseudosource is that we can triviallyguarantee that the source will be reached without any loops. Our goal is to find the geodesic path between a source vertex v s and a target vertex v t on the mesh. Note that it is inefficient to runour exact algorithm on the entire mesh (or even until reaching v t ).In this section we present an algorithm that performs a sequence ofpruned searches, exploiting progressively tighter lower and upperbounds on geodesic distance, so that the final, exact algorithm needonly be run on a “thin” region surrounding the solution.Our approach can be seen as a “continuous A* search”, in thatit adapts the traditional edge-based A* algorithm [Pohl 1971] tomeshes. A similar pruning approach is also explored in [Floateret al. 2002] although their scheme lacks true distance bounds.Denote by P st the geodesic path between v s and v t . Let D s ( p ) and D t ( p ) be the geodesic distances from point p to v s and v t respec-tively, and let D st = D s ( t ) = D t ( s ) be the length of P st . Then, itis obvious that any point p on P st satisfies D s ( p ) + D t ( p ) = D st .If L s ( p ) and L t ( p ) are lower-bound functions of D s ( p ) and D t ( p ) respectively, and U st is an upper-bound value for D st , then anypoint p on P st also satisfies L s ( p ) + L t ( p ) ≤ U st . (5.1)This inequality is the core of the following algorithm: Step 1:
Using Dijkstra search on edges only, compute an upper-bound distance U st ( Dijkstra ) by searching from v s until v t isreached. This step is made almost twice as fast using a bidirec-tional search, which runs two simultaneous Dijkstra searches from v s and v t until they both retire a common vertex. Step 2:
Start our approximation search (Section 4) from v t until v s is reached, which computes a lower-bound distance function L t ( · ) .During the search, we use the inequality (5.1) to prune the search byonly propagating windows that have at least one point p satisfying L t ( p ) + (cid:3) p, v s (cid:3) ≤ U st ( Dijkstra ) , where (cid:3)· , ·(cid:3) measures Euclidean distance (not on the mesh) and istherefore a trivial lower bound on D s ( p ) . Step 3:
Using the windows provided by the previous step, ap-ply backtracing (Section 4.3) to form a path from v s back to v t . Fig. 18. Merging two windows into a new window. s ′ is the pseudosource of the new window, which coversthe visible regions of both merged windows. From [93] window with adjacent windows on the same edge before every propagation step (Figure 18). Tomake sure that the merged window is valid and with bounded error, the two windows to be mergedmust satisfy five conditions: • Directionality.
The two windows must have the same propagation direction with respectto the edge. • Visibility.
To guarantee that the approximated distance field have no gaps, the new windowmust cover the visible regions of the merged windows. • Continuity.
The distance field along an edge must be continuous. Thus, the distances at theendpoints of the new window must be the same as the corresponding ones of the mergedwindows. • Accuracy.
The algorithm bounds the local error of each merging step by a user-specifiedthreshold. , Vol. 1, No. 1, Article . Publication date: July 2020. • Monotonicity.
To guarantee the algorithm’s termination, the new window must have largerdistances than the merged windows. However, the authors did not encounter any infinite loopsin practice and thus they did not perform the corresponding checks in their implementation.In their experiments, the local error threshold is set to be 10%. The test results show that themerging operation is very effective and reduces the WPE (Window Per Edge) to a low numberslightly larger than 1. Thanks to it, their approximated MMP algorithm runs in O ( n log n ) time andoutperforms the Fast Marching method [51] (see Section 3) in both running time and accuracy.To solve the APGD problem efficiently, Xin et al. [105] proposed an algorithm which uses theICH algorithm [104] as a subroutine. The main idea is to embed pre-computed geodesic trianglesin R so that the geodesic distance between two points can be approximated by the Euclideandistance. Their algorithm contains two steps: • Pre-processing step.
In this step, m points are first sampled on Σ . Then, the ICH algorithm[104] is used to compute the geodesic distance field on Σ with the m samples as sources.Utilizing the computed distance field, Delaunay triangulation is applied on Σ according topolyhedral geodesic metric (Figure 19). Finally, the distances between each pair of the m sample points, together with the distances between each mesh vertex and the vertices of thegeodesic triangle containing it, are saved for later use. • Query step.
For any two points on Σ , the two geodesic triangles containing them are firstfound. Then, the two triangles together with the two points are unfolded onto R whilepreserving the geodesic edge lengths. After the unfolding, the geodesic distance between thetwo points is computed as the Euclidean distance between them.From the time cost of the ICH algorithm, it can be easily concluded that the pre-processingstep consumes O ( mn log n ) time. After it, the distance query between any two points on Σ can beanswered in O ( ) time. Intuitively, there is a trade-off between accuracy and computational cost inthis algorithm. That is, the more sample points, the slower the pre-processing and the larger thespace requirements, but more accurate. Shortest Path Construction.
The approximate MMP algorithm [93] employed a similar backtracingstrategy as the original MMP method with only one key difference. Since the windows in it aremerged, the positions of their pseudosources are inaccurate hence not eligible for backtracing. Thus,the authors proposed to maintain a list of the successively merged windows and use the originalpseudosources for backtracing, which yields a more accurate result. Xin et al. [105] employed aquite different strategy to construct shortest paths: since the data structure of their method containsonly geodesic distances but not windows, they are only equipped with an approximated shortestdistance field. Thus, they proposed to construct shortest paths using a gradient descent strategysimilar to the ones used in [28, 51] (see Section 3).
Graph-based methods rely on the assumption that the shortest ge-odesic distance/path between any pair of points p s and p t can be approximated with a chain ofshortest distances/paths ( p s , v , . . . , v k , p t ) , where v , . . . , v k belong to a finite set V G of points of S such that the shortest distance/path between pairs of points of V G is precomputed and encodedin the edges E G of a graph G = ( V G , E G ) . Methods differ for the choice of points of V G end edges of E G and the strategy to build graph G . Once G is given, the PPGP and SSGD problems are easilyresolved through shortest path queries on G , most frequently with standard Dijkstra search [35].Several methods in this class provide data structures that can also support efficient solutions toSSSP or APGD/APSP.The quest for graph-based methods probably started with the PhD thesis of Lanthier [54]. Thebottom line comes from the observation that the shortest geodesic path between a pair of vertices , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 29 (a) (b) (c) (d) (e)
Figure 2:
Algorithm pipeline. Given the input mesh M , we first sample m points distributed uniformly on M (see (a)); Then taking thesample points as sources, we compute the geodesic distance field on M (see (b)) and construct an approximate Delaunay triangulation (see(c)), which induces a geodesic triangulation by replacing each Delaunay edge with a geodesic (see (d)). Next, the exact geodesic distancesbetween any pair of sample points are saved in a file. Finally, with the pre-computed geodesic distance file, we can approximate the geodesicdistance between any pair of query points in constant time O (1) and find the corresponding geodesic path in linear time O ( k ) , where k isthe number of edges crossed by the path (see (e)). Algorithmic pipeline.
Our algorithm contains two steps, the pre-processing step and the query step, as shown in Figure 2.The pre-processing step first samples the given surface with m points, where m is specified by the user (see (a)) and then computesthe geodesic distance fields with the samples as sources (see (b)),which induces a Delaunay triangulation (see (c)). Next, the Delau-nay triangulation edges are replaced by geodesics (see (d)). Finally,we save two kinds of exact geodesic distances into a file, i.e., thedistances between any pair of sample points, and the distances be-tween each mesh vertex and the three vertices of the geodesic tri-angle containing it. The pre-processing step takes O ( mn log n ) time and outputs a file of size O ( m + n ) .In the query step, given two query points q and q on the surface,we first find the geodesic triangles containing them, say T and T (see (e)). Then, we unfold the geodesic triangles T and T to R such that the Euclidean side lengths coincide with the geodesic dis-tances of the curved geodesic triangle. The query points q and q are also mapped to R using unfolding operations. The 2D imagesare denoted by q ′ and q ′ , respectively. Finally, the geodesic dis-tance between q and q is approximated by the Euclidean distancebetween q ′ and q ′ . Note that this approximation can be done inconstant time O (1) as unfolding is a simple and local operation. As shown in Sec. 5.3, the approximation error is closely relat-ed to the maximal side length of geodesic triangles. Thus, anisotropic sampling, which induces geodesic triangulation of sim-ilar sizes, is desired. In this paper, we adopt a variant of thefarthest point sampling method [Peyr ´e and Cohen 2006], which it-eratively updates the sampling set with the farthest point to theexisting sample points. The difference is that we use the im-proved CH algorithm [Xin and Wang 2009] rather than the fastmarching method [Kimmel and Sethian 1998] for geodesic dis-tance computation. The time complexity for sampling m pointsis ! mi =1 f ( ni ) = O ( n log n ) , where f ( x ) = O ( x log x ) .Taking m sample points as source points, we compute themulti-source geodesic distance field using [Xin and Wang 2009],then construct an approximate Delaunay triangulation using[Xin and Wang 2010]. It is shown that the resultant Delaunay tri-angulation is well defined in the sense that all the Delaunay edges do not intersect [Xin and Wang 2010] under moderate condition.However, the edges are not “straight” in general. Finally, we needto replace each curved Delaunay edge with a geodesic path.Note that it takes O ( mn log n ) time to run the exact geodesicalgorithm [Xin and Wang 2009] with each sample point as sourcepoint. Also, it takes O ( n ) time to replace each curved Delaunayedge with a geodesic. Considering the number of edges in the De-launay triangulation is O ( m ) , the time complexity for constructinggeodesic triangulation is O ( mn log n ) . After constructing the geodesic triangulation, we are ready to gen-erate a file containing the following information, which will be usedour all-pairs geodesic distance query algorithm. • The exact geodesic distances between any pair of samplepoints. • For each vertex v , output the geodesic triangle △ s s s containing v and the exact geodesic distances d ( v, s i ) , i =1 , , .Note that all the required geodesic distances are readily availablewhen we constructed the geodesic triangulation in the previous step.Putting them all together, the entire pre-processing takes O ( mn log n ) time and outputs a file of size O ( m + n ) . We first present two simple unfolding operations that play key rolesin our constant-time query algorithm. Observing that geodesic dis-tances follow triangle inequality, thus, we can flatten the curvedgeodesic triangle to a unique triangle (up to rigid motion) in R .Given two curved geodesic triangles △ psr and △ qrs , we can un-fold them to 2D triangles △ p ′ s ′ r ′ and △ q ′ r ′ s ′ , where the corre-sponding Euclidean edge lengths are equal to those of the geodesicedges. Let u ( p, q | rs ) denote such unfolding operation with respectto edge rs . We call it one-side unfolding if p and q are on the sameside of rs and two-side unfolding, otherwise (see Figure 3(a)). Fig. 19. The Delaunay triangulation of a mesh according to polyhedral geodesic metric. From [105]Fig. 20. A shortest path (black line) and a corresponding distance field (color map) computed on the edgegraph (left); and the exact solution (right). From [19]. of mesh Σ can be approximated by a shortest path on the graph of edges of Σ . This approximation,however, turns out to be poor. In fact, such a path is not allowed to cross the faces of Σ , while it is , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 21. The face graph of a triangle for the fixed scheme with two Steiner points per edge (left). The portionof shortest path ab crossing triangle W f i is approximated with edge cd of E G (right). From [56]. forced to meander about their boundaries, thus resulting in a zigzag walk that most of the timesis far more wiggly than the exact geodesic path (see Figure 20 for an example). Starting at thisobservation, Lanthier et al. [55, 56] present three strategies to extend the edge graph to a graph G that incorporates paths across the faces of Σ . They build the set of nodes V G by distributing Steinerpoints along the edges of Σ , and interconnect them with arcs E G that walk either along edges, oracross faces: • The fixed scheme distributes a fixed number of Steiner points uniformly along each edgeand defines V G as the union of all the vertices of Σ and all the Steiner points. Next, for eachtriangle t , it interconnects all the vertices and Steiner points on the boundary of t to forma so-called face graph , which connects each node (either a vertex or a Steiner point) withall the nodes belonging to the other two edges of t and to its two neighbors along the edgeit belongs to. In practice, the face graph of t is a complete graph, except for omitting arcsbetween collinear nodes that are not immediately adjacent along edges of t (see Figure 21left). Graph G is then obtained by collecting all the face graphs. Now we known that an exactpolyhedral geodesic path is a polyline having its nodes at either vertices or edges of Σ [71].The segment ab of one such path traversing a given triangle t is approximated by an edge cd in the face graph of t , where c and d are the closest nodes in the face graph of t along theedges that contain a and b , respectively (see Figure 21 right). • The interval scheme is similar to the previous one, but Steiner points are distributed atuniform distance along each edge, so that the number of Steiner points per edge dependson edge length. Lanthier at al. show that this scheme can achieve the same accurcay of theprevious one with a smaller number of Steiner points, i.e. , a more compact graph G (seeFigure 22). • The spanner scheme uses the same Steiner points as the interval scheme, but it builds a moresparse graph G . Instead of connecting each node v in the face graph to all other visible nodesin the triangle, a predefined set of cones of given width is considered, which emanate from v ,and v is connected to at most one other node per cone (see Figure 23). In this way, graph G isguaranteed to have a number of arcs that is at most a fixed multiple of | V G | .In all schemes, shortest paths are computed by a Dijkstra search on G : SSGD can be resolvedwith a complete visit of G , while for PPGP search can be pruned as soon as distance at the targetpoint becomes definite. SSSP is supported by recording during search the predecessor of each node , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 31
Fig. 22. Approximated paths computer with the fixed scheme (left) and the interval scheme (right): π ( s , t ) isthe shortest geodesic path, π ′ ( s , t ) is the approximated path. From [56].Fig. 23. Spanner edges from point v j (thick red lines) for cones of width ◦ ; thin red lines bound cones.From [56]. in the path to the source: in this way, shortest paths can be retrieved by back-tracing from anyvertex to the source.Lanthier et al. [56] prove that if the number of Steiner points is large enough, they can approximatethe geodesic shortest path within an additive bound that is a function of the length of the longestedge in Σ . The interval scheme results more accurate, while the spanner scheme results morecompact and faster, at the cost of some loss in accuracy. Theoretical bounds are conservative,though, and reasonable bounds can be guaranteed only with a very large number of Steiner points,which is impractical for the applications. On the other hand, they perform extensive testing of thedifferent schemes on polyhedral terrains and provide empirical evidence that very accurate resultsare achieved with just an average of 5-6 Steiner points per edge. More recent comparative tests givenin [19] suggest that the situation may be less favorable if the distance metric is anisotropic: in thatcase, they report the need of 10-20 Steiner points per edge on average to achieve accurate results.Nevertheless, Lanthier’s techniques remain a reference in the context of graph-based methods,because they join practical effectiveness to ease of implementation.In [70], Mata and Mitchell propose a scheme that is similar in spirit to the spanner scheme ofLanthier. Differently from Lanthier, they do not add Steiner points, but they enrich the edge graph , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 24. A cone determined by two rays (continuous paths) from v is split at u ; the path from v to u (dashed)is added as an edge in graph G . From [70]. of Σ with new arcs connecting vertices of V , in such a way that for each vertex v of Σ a ratheruniform set of directions radially arranged about v is explored. They build k equally spaced conesabout v and they propagate the boundaries of such cones according to an exact polyhedral geodesictracing ( e.g. , as in MMP). Propagation is stopped as soon as the two boundaries of a cone intersectdifferent edges of Σ : this means that there exist a vertex u within the cone, which splits it into two;then they determine the shortest path between v and u and add it as an edge to E G . See Figure 24.Note that in this case an edge of E G is not a simple line segment but rather a polyline. Shortestpaths in G are computed through Dijkstra, as in Lanthier’s method.Following the same line of Lanthier’s approach, Aleksandrov et al. [3, 4] proposed differenttechniques to sample Steiner points on edges in geometric progression ( i.e. , more dense close to thevertices and coarser towards the mid edge). The density – hence the number – of sampled pointsdepends on a user specified parameter ϵ and they guarantee that the length of the approximatedpaths is within a factor ( + ϵ ) of the shortest geodesic path. The same authors, in [5], proposea different strategy that samples Steiner points along the bisectors of triangles, rather than ontriangle edges, with a similar geometric progression. They show that this strategy achieves thesame ( + ϵ ) approximation factor with a better time complexity.The results in [3–5, 70] have mostly a theoretical interest, as the number of Steiner verticesnecessary to warrant a small ϵ is too large for practical purposes. Unfortunately, no experimentalresults were provided to show whether the sampling techniques proposed in such works canprovide better results than the original Lanthier’s techniques in practice.Schmidt et al. [90] build a discrete exponential map (DEM) approximation in order to parametrizenormal patches of surface. In the context of the problems we analyze here, they are interestedin resolving SSGD (not SSSP) within the limited scope of a normal neighborhood of the source.In order to estimate geodesic distances, they rely solely on the edge graph of Σ . Once a shortestpath from source p to a generic vertex v has been found in the edge graph, they consider edgesin the path as displacement vectors and they use parallel transport to bring all such vectors to a , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 33
Fig. 25. The normal coordinates u p , q of the unknown radial geodesic from p to q in (a) can be approximatedusing the known geodesics from p to r and r to q . The vector to u r , q (in normal coordinates at r ) is transferredto the tangent plane at p using a 2D rotation with angle θ p , r , producing the red vector in (b). This vector isan approximation to ( u p , q − u p , r ) and can be added to u p , r (c) to get the approximate result ˆ u p , q . From [90].Fig. 26. The DEM with origin at the center of the checkerboard is computed on different discretizations ofthe sphere (below); sampling rate and mesh quality have small effects on this decal. From [90]. common frame. Parallel transport is implemented with a pair of rotations that align local frames atthe vertices along the path. See Figure 25. Next they estimate the geodesic distance between p and v as the length of the vectorial sum of all such displacements. The length and orientation of theresulting vector in fact provide the polar coordinates of v in the DEM centered at p . In practice,this mechanism attempts to “straighten" the wiggly path from the edge graph. Note that, however,they do not compute the straightened path, but only estimate its length. No theoretical analysis ofaccuracy is provided in [90], but empirical tests show good results and resiliency to the quality ofmeshing on reasonably smooth shapes. See Figure 26. , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 27. Unfolding of edge chains to the plane. Edge lengths and relative 1-ring angles are preserved. Thesum vector E (red) is then subdivided according to the orthoprojection of the individual edges. The resultingportions are measured by the respective norms ? here visualized as tensor ellipses (blue) ? and their lengthssummed to get the weight of E . From [19]. Campen et al. [19] elaborate on the same approach in the presence of an anisotropic metric,aiming at a method to resolve SSGD on the whole surface. They start from the observation that thestraightened distance from a path in the edge graph introduces a shortcut that cannot take intoaccount either the presence of holes (boundaries) in the domain, or geometric variations of thesurface (bumps). Hence, the straightening approach by Schmidt et al. cannot be accurate outsidea normal neighborhood. In order to overcome this limitation, they propose a
Short-Term VectorDijkstra algorithm (STVD) that applies shortcuts locally in the context of a standard Dijkstra searchin the edge graph. During the relax phase of Dijkstra algorithm, distance from source s to vertex w is updated by taking the minimum between the distances of each of the k predecessors of w along the path, summed to the shortcut, computed by edge unfolding, from w to that predecessor.In practice, the vector shortcut of Schmidt et al. is applied just locally, on a sliding window oflength k that moves along the path. Values of k in the range between 5 and 10 are reported in theexperiments, with higher values used for higher anisotropies of the input metric. Shortcuts arecomputed with a technique different from [90]: the edge path is unfolded to the plane by preservingthe length of edges, while mapping the angles between pairs of consecutive edges according tothe exponential map at the common vertex, computed as in Fig. ?? . Note that this procedure isfully intrinsic, as the rescaling of angles depends only on the total angle about each vertex in Σ .The length of a shortcut is measured by the Euclidean distance between the source and the targetin the unfolded path. In the presence of anisotropic distance, that shortcut may be segmented inorder to apply different weights in the segments corresponding to different edges. See Figure 27.Campen et al. [19] compare this approach on a single example with anisotropic metric, with respectto several other methods, among which FMM, HM, OUM, and Lanthier’s, by taking as reference asolution computed with Lanthier’s method with 200 Steiner vertices per edge (which is assumed tobe nearly exact). They report that only Lanthier’s with 10-20 Steiner points per edge beats STDVin terms of accuracy. They also comment that the advantage of STDV is much less evident foran isotropic metric; in particular, they report that for low anisotropies, the FMM applied to theintrinsic Delaunay triangulation of a subdivided version of the input mesh is very competitive; noexperiment is shown for this case.Ying et al. [108] start from a basic fact stated in [71]: each exact (polyhedral) geodesic path is apolyline having its internal joints either on edges, or at saddle vertices of Σ (a saddle vertex is avertex having negative Gaussian curvature K - see Section 2.2). A path is said to be direct if it doesnot contain any intermediate saddle vertex. On the basis of this observation, they build the SaddleVertex Graph (SVG) G having the vertices of Σ as set of nodes V G , and one arc in E G for any direct , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 35 path between a pair of vertices. Note that each arc in E G is not a single line segment but a polygonalpath. Such paths need not be stored as they can be back-traced by unfolding triangles starting atthe destination: only a starting direction per edge needs to be stored. If the SVG contains all directpaths, then the exact shortest geodesic path between any pair of vertices can be found as a shortestpath on the SVG. In this respect, the SVG can be considered a data structure to support the APSPproblem, as well as all the other boundary problems. In principle, the SVG could be a completegraph though (hence intractable in the applications, at least on large meshes). Ying et al. showempirically that real world models contain a large fraction of saddle vertices (between 40% and60%) and most paths ( − . ) are not direct. Nevertheless, a full SVG may still be too large forpractical purposes. Ying et al. overcome this limitation by computing a sparsified SVG as follows:for each vertex v they consider only a geodesic disc containing K other vertices and only directpaths to such vertices contribute to the SVG. In this way, K is the maximum degree of nodes in G .The SVG is built by running an exact method (ICH/MMP) from each vertex, pruning the search assoon as the first K vertices have finalized their distance. Variations of Dijkstra search are discussedin the paper, which resolve different boundary problems after the SVG has been built. They runexperiments with values of K between 8 and 1000 and they report the mean approximation errorof shortest paths from SVG to be about 0 .
1% for K =
50. The construction of the SVG can bequite time expensive for high values of K , but they can achieve reasonable times by using a GPUimplementation that computes independently in parallel paths starting at different sources, bysharing just the mesh in read-only mode among the different threads. Concerning query times,Ying et al. compare SVG search with the heat method (with pre-conditioned matrix): they reportthat the SVG is faster for K < < K <
500 and slower, but more accurate,for larger values of K ; they also observe that the SVG is better scalable to large meshes. They alsocompare SVG to the GTU method [105], reporting that SVG is faster and requires less memory; nocomparison about accuracy is reported. Fig. 28. Window propagation in the exact method (left) and in the construction of DGG (right): window W is small enough to stop propagation and allow the path from s to any vertex t in the cone beyond W to travelthrough vertex u . From [99]. Wang et al. [99] define the
Discrete Geodesic Graph
DGG, which improves over the SVG byallowing approximated paths to go also through non-saddle vertices. They start from the keyobservation that cones resulting from window propagation in the exact MMP/CH method becomeprogressively narrower. If a cone is long enough, than any direct path from the source s to a vertex t in the cone can be approximated by the sum of two shortest paths su and ut , where u is a vertexpreceding t and lying on the boundary of the cone. See figure 28. In particular, they show that ifthe geodesic distance between s and t is larger than a certain threshold, which depends on theparameters defining a window W that bounds the cone and on a tolerance threshold ϵ , then the , Vol. 1, No. 1, Article . Publication date: July 2020. approximation error is within a ( + ϵ ) multiplicative factor. So, while generating direct paths as inthe SVG, they set an early termination of window propagation when the length of the cone reachesthe above threshold. Note that a direct path may be eventually approximated with a chain of pathsthrough non-saddle vertices, therefore the approximation rate of a generic path extracted fromthe DGG will not be within ( + ϵ ) , but it remains O ( ϵ ) in any case. The DGG is built through twofundamental steps, for which detailed pseudo-code is provided in [99]: first the standard MMPwindow propagation, with early termination as described above, is run from all vertices; an arc inthe DGG is generated for each vertex reached by a window; next, the resulting set of candidate arcsis pruned by deleting all arcs that can be approximated with sequences of other arcs in the graphwithin the given tolerance; this latter step is performed by running a standard Dijkstra searchfrom each vertex. For anisotropic meshes, in which the fan of arcs incident at some vertices mighthave gaps, an additional step is run to generate additional arcs (see [99] for details). Shortest pathqueries on the DGG are resolved with a SLF/LLL heuristic [11], combined with a technique thatrestricts the neighbors to visit during the relax phase, based on the fact that the angle betweenconsecutive edges in a shortest path cannot be smaller than π − arcsin √ ϵ (see [99] for details). Thelatter pruning technique is reminiscent of an analogous technique presented in [4], based on Snell’slaw, for the case of weighted distances. Extensive experiments are presented in [99] on meshesup to 3M vertices and with accuracies up to ϵ = . i.e. , both boundary and internal vertices of the patchare connected via shortcuts). The shortcuts in the graph are computed with an exact method (inthe ICH/MMP family). The distance between a pair of vertices is found by a Dijkstra search thatis pruned by visiting the graph in a hierarchical manner: search starts from the patch containingthe source at the lowest level of the hierarchy, until it meet its boundary; then it traverses eitherits sibling patch beyond the boundary, or shortcuts from patches in the upper levels, until a patchcontaining the target is found; the path is concluded by moving down the hierarchy until thepatch at the lowest level, which contains the target, is found. See [2] for details. The authors reportquite slow preprocessing times for building the graph, but fast performances and empirically goodapproximations during queries. Local methods start from an initial path connecting two endpoints and aim to refine it to a geodesicpath. These methods are usually iterative, and work either updating one vertex at a time (Sec. 4.2.2)or the whole path at each iteration (Sec. 4.2.3). Point-to-point problems of this kind (PPGP) couldin principle be solved as a byproduct of global methods that strive for SSGP (Sec. 4.1). Nevertheless,there is a number of algorithms that are specifically tailored for it. One good reason to preferlocal methods to compute PPGP is that they are in general faster, easier to code, and require less , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 37 memory due to the smaller domain they consider. On the negative side, they only converge to alocal minimum, producing a locally optimal geodesic path at best. It follows that their ability to findthe global optimum depends on the initialization, which may come from the Dijkstra’s algorithm,from heuristics, or from user interaction.For all methods, the input is a polyhedral mesh and a piece-wise linear path P = p , p , . . . , p n connecting a source p with a destination p n . The computational domain amounts to all the meshvertices, edges and faces traversed by P . The scalability of local methods does not depend on meshsize, but rather on the number of elements traversed by P , which in turn depends on its length andthe local mesh density. A particular instance of the problem occurs when source and destinationcoincide (i.e., the path is a closed loop). Most of the local methods natively support this special case,but there are also methods which are specifically designed to shrink closed loops on a mesh [102].Shrinking loops is useful in a number of scenarios, for example to refine previously computedhomology bases or the set of handles and loops of a mesh (see [32] and references therein). Key to many methods is the ability to locallyassess whether a path is a geodesic or not. To this end, the most important thing to remember isthat while in the continuous the concepts of straightest and shortest paths coincide, in the discretesetting the two concepts are not always equivalent (Sec. 2.2.1). Polthier and Schmies [81] observedthat a path passing through a vertex p divides its total angle θ into two components ( θ l and θ r ) anddefined geodesic paths in terms of these quantities (Fig. 29). In particular, they state that a path is:- locally straightest, if θ l = θ r - locally shortest, if θ l ≥ π and θ r ≥ π One can easily observe that if the path traverses a flat area or an edge (i.e., if θ is 2 π ) the twodefinitions above are equivalent, hence the parallel with the smooth theory still holds. Consideringnon euclidean vertices, we obtain that: (i) a straightest path passing through a spherical vertexcannot be locally shortest (if θ l + θ r < π , then θ l and θ r cannot be both greater or equal than π );(ii) there are infinite shortest paths passing through a hyperbolic vertex (any solution of θ l + θ r = θ ,with θ l ≥ π and θ r ≥ π defines a shortest path). Motivated by the problems of non existenceand non uniqueness of locally shortest paths, Polthier and Schmies advocate the use of locallystraightest paths as a way to trace geodesics on discrete surfaces, and propose both Euler andRunge-Kutta integration schemes based on them [81, 82]. The beauty of locally straightest geodesicsis that they exist and are uniquely defined at any point of a polyhedral mesh. However, there is aprice to pay for this. In particular, it must be observed that straightest geodesics do not convergeto geodesic paths on smooth surfaces under mesh refinement [53, 58], and that locally straightestdistances do not satisfy the triangular inequality, therefore the straightest geodesic distance is not ametric [61]. Last but not least, we remind the reader that these local geodesic criteria apply only tothe Euclidean embedding, and do not extend to alternative metrics (e.g., weighted or anisotropic),thus limiting the applicability of the algorithms that are based on them. Early methods to transform a general path into a geodesic path areheavily based on the local geodesic criteria introduced in Sec. 4.2.1. Given a path P = p , p , . . . , p n ,these methods work by iteratively flattening the mesh around each point p i , and updating itsposition in order to locally straighten or shorten the path. The whole procedure is repeated untilconvergence, that is, until each point locally satisfies its geodesic criterion. Since at each iterationonly one point is updated, these methods tend to require a high number of iterations to reachconvergence. However, the computational cost of each iteration is usually low, as both the localflattening and the vertex update do not involve time consuming operations. , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 29. Left: the total angle θ at point p is divided into left angle θ l = (cid:205) i = θ i and right angle θ r = (cid:205) i = θ i by the path that traverses it. The path is: (i) locally shortest if θ l ≥ π and θ r ≥ π ; (ii) locally straightest if θ l = θ r . Right: three different paths connecting pairs of vertices in a cuboid. The path passing through thespherical vertex is locally straightest, but cannot be locally shortest as its total angle θ is less than π . Thepath on a facet and the one crossing an edge are both locally shortest and straightest, because θ l = θ r = π .Right image from [81]. CyberTape [98] strives for locally shortest geodesics, flattening the mesh around each point p i and tracing a straight line going from p i − to p i + . If the line crosses new edges or vertices,intersection points are added to the path. Depending on where you are, there are multiple waysto locally shorten a path. Fig. 30 shows all the update operators used by the algorithm to handlepoints at edges, spherical vertices, and hyperbolic vertices. Corner cases are also handled: if thelocal update would move the path outside the area spanned by the flattening, it is automaticallyconstrained to adhere to the border of the local domain. This may happen any time the border ofthe flattening is non convex (see an example at the top right corner of Fig. 30). Vertex operators arechosen by querying a look-up table indexed according to the type (euclidean, spherical, hyperbolic)and the relation between θ l and θ r . The algorithm converges when all the points satisfy the geodesiccriterion, which means that the path is locally shortest everywhere. It is interesting to note that,however, the final path may unnecessarily deviate from the input path P . This is because localoperators move the path away from hyperbolic points even if locally shortest paths may traversethem, thus introducing an unnecessary deviation from P . One year later, Martínez et al. [69]published a similar method, which substitutes the lookup table with a more faithful implementationof the ideas of [81]. Specifically, the path is not locally updated at hyperbolic vertices if both θ l and θ r are greater than π , as the path cannot be locally shortened (in [98] the path was pushed awayfrom the vertex anyways). This allows [69] not only to converge to a path that is locally shortesteverywhere, but also to find the one that deviates from the input path P only if strictly necessary.Finally, Xin et al. [103] strive for locally shortest geodesics, but rather than using the local anglecriteria expressed in [81], they adopt an equivalent concept based on the Fermat principle, whichstates that light always follows the shortest optical path. As in [98] and [69], the visibility-basedmethod is guaranteed to converge to a path that is locally shortest everywhere in the sense of [81]. Remark:
A subtle problem affecting all the methods presented in this section is that they do notguarantee convergence to the closest local minimum. In other words, among all the locally shortestpaths connecting source with destination, the output may not be the path which is closest to theinput. Let us consider a path that traverses a spherical vertex in a way that perfectly halves its total , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 39 p i-1 p i p i+1 p i+1 p i-1 p i p i+1 p i+1 p i-1 p i-1 p i+1 p i >p i-1 p i-1 p i+1 p i < <
Graph-based methods can also be designed to work in a local fashion.Similarly to the methods described in Sec. 4.1.3, local graph-based methods find the shortest pathconnecting two points running Dijkstra’s algorithm on a graph, which is initialized with thevertices and edges of the input mesh, and progressively refined adding Steiner points at mesh edges.Differently from global methods, Kanai and Suzuki [49] refine the graph only around the initialedge path connecting source to destination, thus producing a smaller graph and leading to a fasterand better scalable implementation. Being based on Dijkstra’s algorithm, weighted metrics can beeasily incorporated by associating different weights to the edges in the graph. Lanthier et al. [56]present a local heuristic method to refine a given path to an approximated locally shortest path.Their method is a local variant of the global one presented at the beginning of Sec. 4.1.3, and canbe used with both the fixed and the interval scheme. As for the other local methods, the path is , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 41 restricted to a sub-domain composed of all the mesh faces traversed by the initial path. If the initialpath passes through a mesh vertex, the authors propose a heuristic to decide whether the new pathshould pass it to the left or to the right. As already observed at the end of Sec. 4.2.2, these heuristicchoices impact the output result, possibly resulting in a geodesic line that converged to a bad localminimum. Also exact algorithms can be adapted to the PPGP problem. In the second part of theirarticle, Surazhsky et al. [93] proposed a local variant of their MMP implementation optimized tocompute exact point to point shortest paths connecting a source s and a target t . This variant isbased upon an aggressive pruning strategy, which avoids propagating windows that have at leastone point p not satisfying the inequality L s ( p ) + L t ( p ) ≤ U st , where L s ( p ) , L t ( p ) are the lower bounds on the distances d ( s , p ) and d ( s , t ) measured on the mesh,respectively, while U st is the upper bound on the path length, obtained with Dijkstra. The algorithmworks in two steps: at the first step the pruning strategy is used to obtain a minimum amount onwindows; in the second step the exact algorithm uses the windows to compute an exact shortestpath between s and t (details in Section 4.1.1) . As observed by the authors, alternative pruningschemes could potentially be used, possibly providing more performing implementations of pointto point MMP. In GT, the problem of computing a geodesic curve on a domain is formulated as an initial valueproblem . As discussed in Sec. 2.1, a curve γ on a surface patch S is a geodesic if its geodesic curvature κ д is zero everywhere. Geodesic curvature vanishes when its projection on the binormal vectoris zero, which in turn means that the curvature vector γ ′′ is parallel with the surface normal of S at any point of γ . Given this premise, tracing the geodesic curve that starts from point p ∈ S and proceeds as straight as possible in direction v amounts to solving a second order ordinarydifferential equation, subject to the following boundary conditions [82] γ ( ) = p (10) γ ′ ( ) = v . (11)There is a variety of methods that aim to solve this problem, which mostly differentiate toeach other for the type of domain they admit as input. Early works in the field were designed tocompute geodesic paths on parametric surfaces, such as NURBS and Bézier patches [67, 77]. Thisis not surprising, as this was the dominant representation for curved surfaces in the design andmanufacturing industry. In this survey we do not provide details on these methods, partly becausethe article is focused on polyhedral meshes, but also because the reference literature dates back to’80s and ’90s, and has already been covered in previous surveys and books. We point the reader tothe book of Patrikalakis and Maekawa [78] for further details on the topic. Restricting to polyhedral surfaces, the topic was pionereed by Polthier and Schmies [81, 82], who laidthe theoretical bases for tracing geodesics on polyhedral surfaces (Sec. 4.2.1). In [81] they proposetwo alternative methods to integrate a field on a discrete mesh, one based on Euler integration, andthe other based on the fourth order Runge-Kutta method. It should be noted that tracing a path bynumerical integration on a surface immersed in ambient space R requires an extra effort to keepthe path on the surface. Some sort of projection method must be devised in order to guarantee it.The concept of straightest geodesics implicitly solves this problem, as the path can be integratedrotating around vertices by some angular distance, thus never leaving the polyhedral surface. This , Vol. 1, No. 1, Article . Publication date: July 2020. avoids the tedious and error prone implementation of geometric projection (e.g. shooting rays orcomputing intersections with the surface), resulting in a more robuste software.Kumar and colleagues [53] observed that the angle criterion proposed by Polthier and Schmiesdoes not take into account surface normals, and suffers when there are abrupt orientation changesbetween adjacent facets. They proposed an alternative technique to trace a locally straightestgeodesic on polyhedral surfaces. Given a point p on a surface S , and a tracing direction v , theydefine the locally straightest path as the intersection between S and the plane passing through p and containing both v and the surface normal at p . The procedure is repeated iteratively, untilthe path hits the boundary of the mesh, or a certain length is reached. The article also offers avery informative comparison between this strategy, the method of Polthier and Schmies [81], anda state of the art method for geodesic tracing on NURBS. The outcome of their experiments isthat: (i) discrete methods suffer from the bias introduced by the tessellation and are consistentlyless accurate than geodesic tracing on NURBS; (ii) considering the surface normal orientationincreases the accuracy of discrete methods. In particular, the authors observed that while for simpledevelopable surfaces all three methods perform equally well and converge on the smooth geodesiccurve, for doubly curved non developable surfaces the method proposed in [53] is comparativelycloser to NURBS-based methods, whereas the one of Polthier and Schmies [81] converges to adifferent curve. Fig. 32. Illustration of an iteration of the discrete tracing method described in [22]. Left: starting from point p i and tangent vector T i , the subsequent point p i + is computed on a Bézier patch (yellow) that interpolatesa triangle facet (violet). Point p i + is then projected back on the polyhedral surface (middle-left) and thetangent vector updated (middle-right). Tangent vector T i + is parallel transported from p i to p i + . The hat isused to distinguish between entities living in the polyhedral mesh (no hat) from entities living in the Bézierpatch (angular hat). Image taken from [22]. In order to fill the gap between continuous and discrete approaches, recent literature proposeshybrid solutions, where the discrete geodesic path γ is computed on a piece-wise smooth metamesh, and then projected onto the underlying polyhedral surface. In [22] the authors fit a threesided Bézier patch that interpolates both vertex positions and normals of each facet of a triangularmesh, and use it as a domain to solve the geodesic tracing problem. Patches are glued side by side,guaranteeing G continuity across common edges. This ensures that the path can seamlessly travelacross mesh edges and vertices. Differently from [81], where path integration is intrinsically linkedto the polyhedral surface, the communication between the mesh and the atlas of Bézier patchesis provided by shooting rays along surface normals. While this is in general not very robust, theauthors handle bad cases by shrinking the integration step size any time the new point does notproject onto the same facet as the previous one. The method works by iteratively repeating thefollowing four steps (Fig. 32): • Starting from a point p i and a tangent direction, the next point p i + is computed in the Bézierpatch by moving along the tangent direction with a user prescribed step size; • The new point p i + is projected on the polyhedral surface; , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 43 • A new tangent vector is numerically computed by solving a first order ODE with the Runge-Kutta method; • The new tangent vector is parallel transported from p i to p i + ;The algorithm stops when the path hits the boundary of the domain, or when some prescribed pathlength is reached. Compared to classical tracing methods on parametric surfaces, this method hasthe advantage to solve a first order ODE on the tangent vector, thus requiring only G continuity.Classical methods solve a second order ODE on point positions, and therefore require a G con-tinuous domain [78]. Furthermore, when compared to purely discrete methods such as [81, 82],this hybrid approach avoids the problems about existence and uniqueness of the solution aroundspherical and hyperbolic vertices (Sec. 4.2.1), and about convergence to the smooth geodesic curveunder mesh refinement. Methods to trace streamlines of a general vector field can also be used for GT. Typically, thesemethods start from a field defined at the vertices of a discrete mesh, and extend it inside each facetby linear interpolation. For geodesic tracing, the starting point is a distance field, then transformedin a piece-wise linear vector field by computing its gradient (see [68] for an overview of differenttechniques for gradient field computation). Typical tracing techniques are based on numericalintegration, performed with Euler or the fourth order Runge-Kutta method [81]. These numericalapproaches however are problematic, because the approximation error accumulates along the path,possibly producing drifting streamlines that fail to encode the correct global structure of the field.Recent research has shown thath cumulative error propagation can be avoided by integrating thepath one facet at a time, and imposing proper constraints along edges. Bhatia and colleagues [12]introduced the concept of EdgeMaps, replacing classical integration with linear maps defined atthe edges of each facet. Edges are segmented into portions where the flow can be either inwardor outward. Given an entry point to a facet, the corresponding exit point can be computed byapproximating the flow linearly, which means that the curved paths of the vector field are replacedwith straight lines that robustly indicate where the flow leaves the facet (Fig. 33). In [48] it wasobserved that there exist only 23 possible configurations for a triangle, therefore tracing the flowwithin a facet can be conveniently reduced to a map lookup. EdgeMaps do not completely substitutenumerical integration, which is still used at initialization time to find the split points for edges.Furthermore, edge intervals guarantee no overlap of streamlines within the triangle only if theirmap is monotonic, which is true only up to numerical precision. Ray and Sokolov [85] introduced theconcept of stream meshes, which extends the idea of EdgeMaps. They completely avoid numericalintegration, and rather rely on the field direction only along edges, providing an interval pairingsystem which makes heavy use of arbitrary precision floating point numbers to resolve precisionissues. The result is a more robust tracing method, which is guaranteed to avoid crossing andcollapse between close-by curves. It should be noted that, by avoiding integration within triangles,the lines may deviate from the guiding field, while still preserving the correct global behavior.
Xie et al. [101] offer a different view on geodesic tracing, which is employed to parallel transportdeformations across surfaces, essentially solving a PPGP problem where the starting point p represents the input shape, the final point p the target shape, and the geodesic path γ connectingthem realizes a morphing between the two. The authors propose two alternative methods tocompute γ : the shooting method and the straightening method. Both methods are iterative, andgenerate a sequence of paths γ , . . . , γ n that progressively approximate the exact geodesic path , Vol. 1, No. 1, Article . Publication date: July 2020. Fig. 33. EdgeMaps start from a vector field defined on vertices and linearly interpolated inside each triangle(left), and split triangle edges into slabs where the flow is either incoming or outcoming (middle). Theinformation is eventually transformed into a linear map, which allows to easily pair entry and exit points(right). In [48] it was shown that for triangles there exist only 23 possible flow configurations, therefore allpossible maps can be efficiently encoded into a lookup table. connecting p and p . The shooting method works by shooting a discrete geodesic path γ withfinite number of steps starting from p along a random direction v , and to iteratively update theshooting direction v in order to converge to the target path. Assuming all γ i are parameterizedby arc length, the update amounts to computing the vector w = p − γ i ( ) , parallel transport itfrom γ i ( ) to γ i ( ) , and set v i + = v i + δw . The process is repeated until convergence (i.e. until ∥ w ∥ goes to zero). A visual example of the shooting procedure is given in Fig. 34 (top line). Thestraightening method works by initializing γ as an arbitrary path connecting p and p , anditeratively straightening it using a gradient descent approach, minimizing the energy function E ( γ ) = ∫ (cid:28) ddt γ , ddt γ (cid:29) dt . The inner product inside the integral is defined using the Riemannian metric of the shape space.Given the gradient ∇ E , the path is updated as γ i + = γ i − δ ∇ E , where δ > A key issue that arises in all methods for computing geodesic distance (including those for computingthe exact polyhedral distance) is that the surface must be meshed appropriately in order to beable to carry out the necessary computations. The meshing of the surface will have an impact onthings like numerical accuracy, guarantees about properties of the solution, as well as the timeand memory needed to execute the algorithm. In the context of PDE-based methods (Sec. 3), themesh is explicitly given to the algorithm and should ideally satisfy familiar properties demandedby finite element methods for hyperbolic and elliptic problems. In the context of computationalgeometry methods (Sec. 4), the mesh is implicitly constructed during the course of the algorithm– viewed through this lens, many of the strategies discussed in Sec. 4 can be viewed as differentstrategies for keeping the size and quality of this mesh under control. This viewpoint also helpsto understand the trade offs between different classes of methods: PDE-based methods are oftenmore efficient because they can generate a single, high-quality mesh ahead of time and re-use thismesh for many different distance queries; however, this mesh cannot be adapted to a particular , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 45
Fig. 34. Xie et al. [101] propose two iterative methods that use geodesic tracing to solve a PPGP problem. Topline shows the shooting method. Bottom line shows the straightening method. The user chooses the numberof sample points the path should be composed of ( n = , , ), and the algorithms iteratively update aninitial path (green) to make it converge (blue) to the exact geodesic obtained from analytical expressions (red).All geodesics are traced in the 2D hyperbolic space. Image taken from [101]. distance query a priori . In contrast, methods from computational geometry can provide exactsolutions to the polyhedral geodesic distance problem because the mesh implicitly constructedby the algorithm can be explicitly tailored to a particular query point, i.e. , the surface is meshedexactly along characteristics of the solution; however, this same mesh cannot easily be re-used forsubsequent queries, leading to a large number of caching and approximation schemes (as discussedin Sec. 4). In this section we consider the relationship between meshing and accurate geodesicdistance computation. The two basic classes of methods we have studied (hyperbolic and elliptic) naturally lead to twodifferent basic criteria on meshes needed to generate accurate results: hyperbolic methods generallyperform better on acute triangulations; elliptic methods perform better on
Delaunay triangulations.On the whole, these criteria should be viewed as “rules of thumb” – for instance, there are plentyof non-Delaunay meshes on which elliptic methods still perform well. However, in specific casesusing a mesh with the desired characteristics will provide absolute guarantees about properties ofthe solution, e.g. , absence of spurious local minima. Interestingly, acute meshes are a strict subsetof Delaunay meshes; hence, meshes that work well for hyperbolic methods will also work quitewell for elliptic methods.
As noted in Sec. 3.2.1, methods based on wavefront propagation mustupdate distance values at vertices in an order that respects causality, i.e. , the active node of smallestvalue must always be the one closest to the source. On a triangulated surface, this property canbe violated in the presence of obtuse angles – consider, for instance, executing the fast marchingmethod on the triangulation pictured in Fig. 8, left . If we place a source at vertex 1, then the distance , Vol. 1, No. 1, Article . Publication date: July 2020. at vertices 2 and 3 will be computed (and finalized) prior to vertex 4, even though 4 is closer to 1; asa result, we will get an overestimate of the distance at 4. This phenomenon cannot occur unless themesh has obtuse angles; a solution, therefore, is to mesh the domain such that it has only acute (orat least nonobtuse) angles, i.e. , interior angles no greater than π / geodesic triangulation can be obtained for any polyhedron ( i.e. , a triangulation where edges are requiredonly to be geodesic arcs along the polyhedron), though this triangulation may not include theedges of the polyhedron itself. Saraf [88] later proves the existence of acute triangulations whichcontain the edges by constructing compatible acute triangulations on each face; in either case,no upper bound is given on the number of triangles. More recently, Maehara [66] give an upperbound on the size of the triangulation in terms of the local geometry of the given polyhedron,though no practical algorithm is provided. In higher dimensions ( e.g. , for computing geodesicson triangulated regions of R ) even less is known, and in general acute triangulations may noteven exist; Brandts et al. [16] provide a survey. The upshot of this discussion is that (at the timeof writing of this article) there is no way to generate meshes that are guaranteed to satisfy therequirements of wavefront-based methods on polyhedral models, i.e. , no absolute guarantee thatthe upwind ordering will be monotonic with respect to distance. From here, there are essentiallythree options: (i) apply an unfolding procedure, though as noted in Sec. 3.2.1 this procedure may notalways terminate; (ii) iteratively update the solution via multiple sweeps ( à la Bornemann:2004:FED),though at this point one may simply consider other optimization strategies such as ADMM [8],or (iii) simply ignore the fact that the solution may exhibit spurious local minima. In practice,wavefront-based methods perform quite well even on meshes with a some mildly obtuse triangles,and for many applications ( e.g. , visualization) small violations of monotonicity do not present majorproblems.
For geodesic distance methods based on solving elliptic questions, the keyissue is no longer causality but rather satisfaction of a smooth maximum principle . In particular, let L ∈ R | V |×| V | be a matrix representing any discrete Laplace operator, i.e. , a weighted graph Laplacianwith edge weights w ij ∈ R . We say that a function ϕ ∈ R | V | is discrete harmonic (with respect to L ) if it is in the kernel of L ( i.e. , if Lϕ = ∆ ϕ = w ij associated with interior edges arepositive; in particular, when L is the cotangent Laplacian positivity of edge weights is equivalent tothe intrinsic Delaunay condition α ij + β ij > π , where α ij , β ij are the interior angles opposite an interior edge ij . This condition is necessary (but notsufficient) to ensure a global maximum principle—see Wardetzky et al. [100] for further discussion.Violation of the maximum principle can (as with causality) yield spurious local minima in thedistance function computed by elliptic methods.Fortunately, Delaunay triangulations are far easier to obtain than acute triangulations, and thereis a large body of knowledge about both theoretical guarantees and practical algorithms – see for , Vol. 1, No. 1, Article . Publication date: July 2020. Survey of Algorithms for Geodesic Paths and Distances 47
Fig. 35. A mesh with many non-Delaunay triangles (left) may yield inaccurate solutions that depend on thecotan-Laplace operator (right) , but is easily remedied by using a Delaunay triangulation (right) . Here we showthe heat method as implemented in CGAL [97], which uses an intrinsic Delaunay triangulation to avoidincreasing the mesh size. instance the book by Cheng, Dey, and Shewchuk [24]. Less is known about Delaunay triangulationof polyhedral surfaces, though even here there are a variety of practical algorithms, includingDelaunay refinement [23, 33, 34] and edge splits [37, 63]. Such methods can either exactly triangulatethe input geometry, or reduce the size of the triangulation by allowing small modifications to thegeometry. Alternatively, one can construct an intrinsic Delaunay triangulation [13, 40] where thenew triangulation is represented only by the mesh connectivity and a collection of edge lengths; thisdata is sufficient to construct the Laplace operator L , and allows greater flexibility in the definitionof the triangulation while still preserving the input geometry exactly. Fig. ?? shows one example ofhow Delaunay triangulations can significantly improve the accuracy of elliptic methods – in thiscase, without even increasing the size of the mesh. One appeal of the methods considered in Sec. 4 is that they can work directly on a given inputwithout remeshing, though in reality these methods still effectively re-mesh the surface via theconstruction of windows (or other auxiliary data). In other words, they incrementally modify atopological mesh data structure (often, just a standard half edge mesh [93]) that describes a newtessellation of the original mesh. From this point of view, the basic approach of computationalgeometry methods is not entirely different from PDE-based methods: in either case, one mustconstruct a mesh that supports accurate computation of a solution. However there are some notabledifferences, namely (i) window-based methods remesh the domain during computation of geodesicdistance, rather than doing it ahead of time and (ii) the tessellation implicitly generated by thisprocess depends on a particular choice of query point (or points). In particular, the surface iseffectively “meshed along geodesics,” i.e. , it explicitly encodes geodesics from pseudosources towindow boundaries (see for instance Fig. 9). This strategy presents some clear trade offs. On theone hand, it provides a great deal of accuracy since geodesics are directly encoded in the mesh –from the PDE point of view, this is akin to meshing the domain along characteristics of the eikonalequation. It is also therefore much easier to recover geodesic paths. For this same reason, however,the mesh implicitly generated by a window-based scheme is difficult to re-use for computingdistance to a different query point, which will induce an entirely different collection of windows. , Vol. 1, No. 1, Article . Publication date: July 2020.
Method Implem. URL Type (∗) https://github.com/dgpdec C++ API (∗) https://github.com/mlivesu/cinolib C++ API (∗) (*) supports solving by back-substitution on a pre-factored matrix
Table 1. An overview of available implementations. For each method we report: method’s name, reference tothe original papers, link to the code, and its format.
This situation explains the large number of different variants on the basic MMP and CH algorithmsthat are needed to cache or otherwise accelerate new queries, as discussed in Sec. 4.One topic that has not received much attention is the impact of the input triangulation onthe performance of windowing algorithms. One rather surprising observation is that the MMPalgorithm sometimes runs much faster on noisy inputs – see in particular [93, Figure 7]. Why thisshould be true, and more broadly, what kind of input triangulations lead to the best performance is avery interesting question, possibly leading to a hybrid approach where the mesh is first preprocessed( e.g. , to be Delaunay or acute) prior to initiating the window insertion process.
Implementations of several methods reviewed in the previous sections have been released to thepublic domain, either by their authors, or by others. Table 1 provides a summary of the softwarewe could find. Some such implementations come in the form of API’s, some consist of source codefor stand alone applications, and some others are just executables. Running such software in aconsistent framework that allows us to compare their performances is not always possible. Severalmethods involve a setup phase (preprocessing) which is run just once per dataset, and is separatedfrom the query phase; the costs of such two phases must be evaluated separately and this is notalways possible in the implementations we collected. Moreover, as we already discussed previously,methods working in the smooth and in the discrete setting should not be compared directly, asthey address different problems.In the following two subsections, we provide just some comparisons between methods at thestate-of-the-art, focusing on the SSGD problem and showing separate results for the PDE-basedmethods (smooth setting) and the CG methods (discrete setting).
As we reviewed in Sec. 4, most global methods for resolving the SSGD problem belong to twodistinct classes: the exact methods that propagate windows (Sec. 4.1.1), and the graph-based methods(Sec. 4.1.3). In Table 2 we compare the performances of the author’s implementation of the VTPmethod [84], which represents the state-of-the-art in the first class, and our basic implementationof Lanthier’s method [55, 56] with the interval scheme and SLF/LLL heuristic [11] to computeshortest paths. The state-of-the-art for graph-based methods is probably the DGG method [99]. Unfortunately, we only have an executableavailable for this method, which prevents us from profiling performances to compare with the other methods. The authors, Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 49
VTP Lanthier’s S = Lanthier’s S = Lanthier’s S = Lanthier’s S = Dataset Size Setup SSGD Setup SSGD Error Setup SSGD Error Setup SSGD Error Setup SSGD ErrorHand 8k 0.002 0.027 0.031 0.002 0.013 0.082 0.008 0.005 0.108 0.013 0.003 0.245 0.028 0.001Joker 27k 0.016 0.109 0.130 0.015 0.015 0.263 0.033 0.005 0.418 0.053 0.003 0.877 0.116 0.001Bunny 70k 0.042 0.359 0.250 0.035 0.022 0.626 0.086 0.007 0.857 0.133 0.004 1.984 0.294 0.001David 99k 0.055 0.490 0.394 0.050 0.014 0.839 0.116 0.005 1.312 0.185 0.003 2.842 0.432 0.001Armadillo 346k 0.392 2.942 1.876 0.215 0.015 3.569 0.500 0.006 5.441 0.754 0.003 11.945 1.667 0.001Gargoyle 700k 0.402 5.560 2.849 0.372 0.016 6.399 0.902 0.006 8.096 1.415 0.003 20.956 3.068 0.001Blade 2M 2.098 41.600 9.268 1.298 0.012 18.759 2.992 0.005 26.163 4.538 0.003 (*) (*) (*)Happy 2.6M 2.968 37.449 12.777 2.071 0.015 23.994 4.903 0.006 45.299 8.021 0.004 (*) (*) (*)Neptune 2M 3.600 66.840 17.564 2.660 0.015 35.956 6.649 0.006 (*) (*) (*) (*) (*) (*)Lucy7M 14.5M 18.831 567.655 149.795 14.206 0.014 (*) (*) (*) (*) (*) (*) (*) (*) (*)
Table 2. Comparison between VTP [84] and Lanthier’s [55, 56] on various datasets. From the left: dataset nameand size (number of triangles in thousands and millions); VTP times for setup and SSGD query; Lanthier’stimes for setup and SSGD query, and RMS error wrt ground truth. Experiments with Lanthier’s are repeatedwith 1, 3, 5, and 10 Steiner points per edge (average). Times are in seconds; (*) means that the data structureexceeds physical RAM: the program runs, but with very slow times due to paging.
We have run experiments on meshes of increasing size, from a few thousands to about 14 milliontriangles; all meshes are manifold and watertight. Trials were run on a laptop equipped with 2.9Ghz Intel Core i7 CPU and 16 GB RAM by using a single core. For each dataset and method, wedid the setup once, and then ran queries on 100 random seeds. In Table 2, for each method wereport the time to setup and the average time for a single query. The VTP method, which is exact,provides the ground truth for the discrete problem. We compare the accuracy of the approximatedgraph-based method by measuring the average RMS error, evaluated asError = (cid:118)(cid:117)(cid:116) (cid:213) v ∈ V \{ s } (cid:18) dist VT P ( v ) − dist G ( v ) dist VT P ( v ) (cid:19) where V is the set of vertices, s is the source, and dist VT P ( v ) and dist G ( v ) denote the distancesof a vertex v from s computed with the two methods, respectively. For Lanthier’s method, thereis an obvious trade-off between the time performance and the accuracy, which depends on theaverage number of Steiner points per edge. Thus, to make an exhaustive comparison, we ran theexperiments by varying such number between 1 and 10. Higher values were not considered due totheir inefficiency: although the error becomes tiny, the time performance becomes not competitivewith respect to the exact methods.Overall, given a dataset, we assume that the setup is done once and the SSGD queries areperformed many times, which is the common practice in various applications. Therefore, as longas the setup remains within reasonable time bounds, query times are most relevant.For the VTP method, the query times dominate (Table 2). The setup phase of VTP takes negligibletime since it only involves building the basic data structures to support the window propagation.In spite of the theoretically superlinear complexity, the practical query times of VTP appear toincrease linearly with the size of the dataset. The only exception is the larger time consumed onthe Blade dataset, which contains large flat regions. Such performance degeneration may revealthe inherent challenge of the window-based methods: on a flat region, almost all the propagatedwindows are valid and there is little space for acceleration through redundancy removal. In thiscase, all the window-based methods degenerate and behave similarly.On the contrary, the setup phase dominates in Lanthier’s method, which is devoted to buildingthe graph containing both vertices and Steiner points. Note that with 3 or more Steiner points per said that an improved version of DGG is under minor revision for publication in a journal, and they will release the sourcecode as soon as it is accepted. We plan to include further experiments with their code in the final version., Vol. 1, No. 1, Article . Publication date: July 2020. edge, the cost of this phase is comparable or even more expensive than the whole cost of runningVTP. The times for the setup phase increase linearly with the size of the mesh and quadraticallywith the number S of Steiner points per edge (indeed, the size of the graph is quadratic in S , asshown in Fig. 21). Once the graph has been built, query times for meshes of moderate size arebetween one and almost two orders of magnitude faster than VTP, up to the case of 5 Steiner pointsper edge, while they become comparable or slower than VTP by using 10 Steiner points per edge.Speed is paid at the cost of some error, which is about 1 .
5% with 1 Steiner point and about 0 . S = S = S , is inherent in the method. Algorithms for solving various versions of the geodesic problem are quite mature and ready tobe used in applications. We have reviewed methods at the state of the art, which are dividedinto two broad classes: PDE-based methods that resolve the problem in the smooth setting, andcomputational geometry methods that resolve the problem in the discrete setting. The latter classis subdivided further into methods that propagate windows, all stemming from the seminal MMPmethod, and graph methods that introduce one further level of discretization. We have also reviewedlocal methods that address single point-to-point queries and methods for geodesic tracing, as wellas discussed practical aspects such as the criteria that meshes should satisfy in order to squeeze themaximal power from each of the surveyed methods.Several methods, in both classes, support the separation of the setup phase from the query phase,where the former is performed once after loading the data, and the latter may be executed anarbitrary number of times on the same setup, to amortize the computational cost. Depending onmethods, this approach permits to deplete much of the computational burden during pre-processing,thus achieving high speed-ups at query time. Overall, there is no best method for all purposes, aseach class of methods has its own characteristics.PDE-based methods are based on a global approach to the problem, and are largely aimed atsolving the smooth geodesic problem, i.e. , finding the best approximation of the true geodesicdistance on the sampled surface. These methods involve the resolution of large sparse linear systemsand may benefit from modern numerical solvers, which in turn allows them to easily exploit featureslike parallelism. Matrix pre-factorization can be done during the setup phase, thus reducing queriesto fixed-order back-substitution, which is extremely fast. As with any finite element method, theapproximation quality of PDE-based methods will be influenced by the quality of the input mesh,though as discussed in Sec. 6 quality can often be improved via straightforward meshing strategies(such as use of an intrinsic Delaunay mesh). Use of direct solvers can fail for very large meshes(due to limits on memory); in this case a practical solution is to switch to a large-scale iterativesolver (such as multi-grid or preconditioned conjugate gradient). Looking forward, recent work oncomputing localized solutions to large linear system opens the door to using PDE-based methodsfor local geodesic queries ( e.g. , just at a single point) while still leveraging the benefits of directsolvers [45, 46]. Pushing these methods toward higher accuracy via use of higher-order elements isalso an interesting area for future work [30, 72].Computational geometry-based methods are based on local window propagation, and are mostappropriate for solving the polyhedral geodesic problem. Although the seminal algorithms in this , Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 51 category had poor performance due to a large number of windows, careful management of windowconstruction in recent work has led to algorithms with a small memory footprint, thus supportingefficient geodesic queries on large datasets. On the other hand, such methods cannot easily separatethe setup and query phases, largely due to the phenomenon discussed in Sec. 6: splitting edgesinto windows effectively constructs a different mesh of the surface for each geodesic distancequery. Thus, the amortized performance for repeated queries is not immediately as good as, say,PDE-based methods, leading to a focus in recent work on caching or otherwise pre-computinginformation useful for multiple queries [108]. A key question in future work is how to balanceperformance and accuracy with the typically high degree of implementation complexity associatedwith polyhedral algorithms.Graph-based methods provide an approximated solution to the discrete problem and permitto trade-off accuracy for speed. They excel in numerical stability and they are very efficient inseparating setup from query phase, thus achieving a relevant speed-up with respect to windowpropagation methods and times comparable with the fastest PDE-based methods. On the otherhand, the memory footprint to store graphs may become relevant, and even prohibitive if highaccuracy is required on large meshes, thus making such methods not competitive with the exactones.
REFERENCES [1] Shahnawaz Ahmed, Stanley Bak, Joyce R. McLaughlin, and Daniel Renzi. 2011. A Third Order Accurate FastMarching Method for the Eikonal Equation in Two Dimensions.
SIAM J. Scientific Computing
33, 5 (2011), 2402–2420.https://doi.org/10.1137/10080258X[2] Rosario Aiello, Francesco Banterle, Nico Pietroni, Luigi Malomo, Paolo Cignoni, and Roberto Scopigno. 2015. Com-pression and Querying of Arbitrary Geodesic Distances. In
Image Analysis and Processing — ICIAP 2015 , VittorioMurino and Enrico Puppo (Eds.). Springer International Publishing, 282–293.[3] Lyudmil Aleksandrov, Mark Lanthier, Anil Maheshwari, and Jörg-R Sack. 1998. An ε -Approximation algorithm forweighted shortest paths on polyhedral surfaces. In Scandinavian Workshop on Algorithm Theory . Springer, 11–22.[4] Lyudmil Aleksandrov, Anil Maheshwari, and Jörg-Rüdiger Sack. 2000. Approximation algorithms for geometricshortest path problems. In
Proceedings of the 32nd Annual ACM Symposium on Theory of Computing . ACM, New York,New York, USA, 286–295.[5] L Aleksandrov, A Maheshwari, and J R Sack. 2005. Determining approximate shortest paths on weighted polyhedralsurfaces.
Journal of the ACM (JACM)
52, 1 (Jan. 2005), 25–53.[6] Mukund Balasubramanian, Jonathan R Polimeni, and Eric L Schwartz. 2009. Exact geodesics and shortest paths onpolyhedral surfaces.
IEEE transactions on pattern analysis and machine intelligence
31, 6 (2009), 1006–1016.[7] Mikhail Belkin and Partha Niyogi. 2003. Laplacian Eigenmaps for Dimensionality Reduction and DataRepresentation.
Neural Computation
15, 6 (2003), 1373–1396. https://doi.org/10.1162/089976603321780317arXiv:https://doi.org/10.1162/089976603321780317[8] Alexander G Belyaev and Pierre-Alain Fayolle. 2015. On Variational and PDE-Based Distance Function Approximations.In
Computer Graphics Forum , Vol. 34. Wiley Online Library, 104–118.[9] P.H. Bérard. 1986.
Spectral Geometry: Direct and Inverse Problems . Springer-Verlag. https://books.google.com/books?id=_LoZAQAAIAAJ[10] P. Bérard, G. Besson, and S. Gallot. 1994. Embedding Riemannian manifolds by their heat kernel.
Geometric &Functional Analysis GAFA
4, 4 (01 Jul 1994), 373–398. https://doi.org/10.1007/BF01896401[11] Dimitri P. Bertsekas. 1998.
Network optimization: continuous and discrete models . Athena Scientific.[12] Harsh Bhatia, Shreeraj Jadhav, Peer-Timo Bremer, Guoning Chen, Joshua A Levine, Luis Gustavo Nonato, and ValerioPascucci. 2011. Edge maps: Representing flow with bounded error. In . IEEE, 75–82.[13] A. I. Bobenko and B. A. Springborn. 2005. A discrete Laplace-Beltrami operator for simplicial surfaces.
ArXivMathematics e-prints (March 2005). arXiv:math/0503219[14] F. Bornemann and C. Rasch. 2004. Finite-Element Discretization of Static Hamilton-Jacobi Equations Based on a LocalVariational Principle.
ArXiv Mathematics e-prints (March 2004). arXiv:math/0403517[15] Prosenjit Bose, Anil Maheshwari, Chang Shu, and Stefanie Wuhrer. 2011. A Survey of Geodesic Paths on 3D Surfaces.
Comput. Geom. Theory Appl.
44, 9 (Nov. 2011), 486–498. https://doi.org/10.1016/j.comgeo.2011.05.006, Vol. 1, No. 1, Article . Publication date: July 2020.
Source Name Domain Class Accuracy Metrics Problem [51] FM S D (wavefront) A E SSGD[7] S D (diffusion) A E SSGD[25] Diffusion S D (diffusion) A E SSGD[43] S D (wavefront) A E SSGD[59] Biharmonic S D (diffusion) A E SSGD[65] Fast sweep S D (wavefront) A E SSGD[28 ? ] Heat S D (diffusion) A E SSGD[92] S D (wavefront) A E SSGD[62, 64, 71, 93] MMP P W (priority queue) E E SSGD[20, 50] CH P W (FIFO queue) E E SSGD[93] Approx.MMP P W (priority queue) E/A (local bounds) E SSGD[6] P W (FIFO queue) E E APSP[104] ICH P W (priority queue) E E SSGD[105] GTU P W (priority queue) A (implicit bounds) E APSP[109] PCH P W (parallel) E E SSGD[106] FWP P W (bucket queue) E E SSGD[84] VTP P W (priority queue) E E SSGD[55, 56] Lanthier’s P G (Steiner) A (additive bounds) E, W SSSP[70] P G (vertex-to-vertex arcs) A ( ( + ϵ ) bound) E, W SSSP[3–5] Lanthier’s var. P G (Steiner) A ( ( + ϵ ) bound) E, W SSSP[90] P G (Exp.map, local) A E SSGD[19] STDV P G (edge graph) A A SSGD[108] SVG P G (vertex-to-vertex arcs) E/A E, W APSP[99] DGG P G (vertex-to-vertex arcs) A ( O ( ϵ ) bound) E, W APSP[55, 56] P L (Steiner) A E, W PPGP[49] P L (Steiner) A E,W PPGP[98] CyberTape P L (locally shortest) E (local) E PPGP[69] P L (locally shortest) E (local) E PPGP[103] P L (Fermat’s visibility) E (local) E PPGP[61] P L (numerical opt.) A E, W, A PPGP[81, 82] P T (locally straightest) A E GT[53] P T (locally straightest) A E GT[12] EdgeMaps P T (general fields) A E, W, A GT[85] P T (general fields) A E, W, A GT[101] S T (point to point) A E GT/PPGP[22] S T (Bézier patches) A E GT Table 3. Summary of the methods considered in this survey. For each method we report: the type of domain(S for sampling of a continuous surface or P for polyhedral mesh); its class (D for PDE-based methods, W forglobal methods based on window propagation, G for graph-based methods, L for local methods, T for tracingmethods); its accuracy (E for exact, A for approximated); the distance metric is supports (E for Euclidean,W for weighted, A for anisotropic); The class of problems it addresses (PPGP, SSGD, APGD,GT); and itsasymptotic complexity (if given by the authors). [16] J. Brandts, S. Korotov, M. KÅŹÃŋÅ¿ek, and J. Åăolc. 2009. On Nonobtuse Simplicial Partitions.
SIAM Rev.
51, 2 (2009),317–335. https://doi.org/10.1137/060669073 arXiv:https://doi.org/10.1137/060669073[17] Y. Burago and V. Zallgaller. 1960. Polyhedral Embedding of a Net.
Vestnik Leningrad. Univ.
15 (1960), 66–80.[18] Thomas Caissard, David Coeurjolly, Jacques-Olivier Lachaud, and Tristan Roussillon. 2017. Heat kernel Laplace-Beltrami operator on digital surfaces. In . Walter G. Kropatsch, Ines Janusch and Nicole M. Artner, Springer-Verlag, Vienna,Austria.[19] Marcel Campen, Martin Heistermann, and Leif Kobbelt. 2013. Practical Anisotropic Geodesy. In
Proceedings of theEleventh Eurographics/ACMSIGGRAPH Symposium on Geometry Processing (Genova, Italy) (SGP ’13) . EurographicsAssociation, Aire-la-Ville, Switzerland, Switzerland, 63–71. https://doi.org/10.1111/cgf.12173, Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 53 [20] Jindong Chen and Yijie Han. 1990. Shortest Paths on a Polyhedron. In
Proceedings of the Sixth Annual Symposiumon Computational Geometry (Berkley, California, USA) (SCG ’90) . ACM, New York, NY, USA, 360–369. https://doi.org/10.1145/98524.98601[21] Yanqing Chen, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam. 2008. Algorithm 887:CHOLMOD, Supernodal Sparse Cholesky Factorization and Update/Downdate. (2008).[22] Peng Cheng, Chunyan Miao, Yong-Jin Liu, Changhe Tu, and Ying He. 2016. Solving the Initial Value Problem ofDiscrete Geodesics.
Comput. Aided Des.
70, C (Jan. 2016), 144–152. https://doi.org/10.1016/j.cad.2015.07.012[23] S.-W. Cheng, T. K. Dey, and T. Ray. 2005. Weighted Delaunay Refinement for Polyhedra with Small Angles. In
Proceedings of the 14th International Meshing Roundtable , Byron W. Hanks (Ed.). Springer Berlin Heidelberg, Berlin,Heidelberg, 325–342.[24] Siu-Wing Cheng, Tamal K. Dey, and Jonathan Shewchuk. 2012.
Delaunay Mesh Generation (1st ed.). Chapman &Hall/CRC.[25] Ronald R. Coifman and StÃľphane Lafon. 2006. Diffusion maps.
Applied and Computational Harmonic Analysis
21, 1(2006), 5 – 30. https://doi.org/10.1016/j.acha.2006.04.006 Special Issue: Diffusion Maps and Wavelets.[26] Keenan Crane, Fernando de Goes, Mathieu Desbrun, and Peter SchrÃűder. 2013. Digital Geometry Processing withDiscrete Exterior Calculus. In
ACM SIGGRAPH 2013 courses (Anaheim, California) (SIGGRAPH ’13) . ACM, New York,NY, USA, 126.[27] Keenan Crane and Max Wardetzky. 2017. A Glimpse Into Discrete Differential Geometry.
Notices of the AmericanMathematical Society
64, 10 (November 2017), 1153–1159.[28] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. 2013. Geodesics in heat: A new approach to computingdistance based on heat flow.
ACM Transactions on Graphics (TOG)
32, 5 (2013), 152.[29] Keenan Crane, Clarisse Weischedel, and Max Wardetzky. 2017. The Heat Method for Distance Computation.
Commun.ACM
60, 11 (Oct. 2017), 90–99. https://doi.org/10.1145/3131280[30] Fernando De Goes, Mathieu Desbrun, Mark Meyer, and Tony DeRose. 2016. Subdivision exterior calculus for geometryprocessing.
ACM Transactions on Graphics (TOG)
35, 4 (2016), –11.[31] Fernando de Goes, Beibei Liu, Max Budninskiy, Yiying Tong, and Mathieu Debrun. 2014. Discrete 2-Tensor Fields onTriangulations. 33, 5 (2014).[32] Tamal K Dey, Fengtao Fan, and Yusu Wang. 2013. An efficient computation of handle and tunnel loops via Reebgraphs.
ACM Transactions on Graphics (TOG)
32, 4 (2013), 32.[33] T. K. Dey, J. A. Levine, and A. Slatton. 2010. Localized Delaunay Refinement for Sampling and Mesh-ing.
Computer Graphics Forum
29, 5 (2010), 1723–1732. https://doi.org/10.1111/j.1467-8659.2010.01781.xarXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1467-8659.2010.01781.x[34] Tamal K. Dey and Tathagata Ray. 2010. Polygonal surface remeshing with Delaunay refinement.
Engineering withComputers
26, 3 (01 Jun 2010), 289–301. https://doi.org/10.1007/s00366-009-0162-1[35] E. Dijkstra. 1959. A Note on Two Problems in Connexion with Graphs.
Numer. Math.
Differential Geometry of Curves and Surfaces . Prentice-Hall.[37] Ramsay Dyer, Hao Zhang, and Torsten Möller. 2007. Delaunay Mesh Construction. In
Proceedings of the FifthEurographics Symposium on Geometry Processing (Barcelona, Spain) (SGP ’07) . Eurographics Association, Aire-la-Ville,Switzerland, Switzerland, 273–282. http://dl.acm.org/citation.cfm?id=1281991.1282027[38] H. Erten and A. Üngör. 2009. Computing Triangulations without Small and Large Angles. In . 192–201. https://doi.org/10.1109/ISVD.2009.32[39] H. Erten and A. Üngör. 2009. Quality Triangulations with Locally Optimal Steiner Points.
SIAM Journal on ScientificComputing
31, 3 (2009), 2103–2130. https://doi.org/10.1137/080716748[40] Matthew Fisher, Boris Springborn, Alexander I. Bobenko, and Peter Schroder. 2006. An Algorithm for the Constructionof Intrinsic Delaunay Triangulations with Applications to Digital Geometry Processing. In
ACM SIGGRAPH 2006Courses (Boston, Massachusetts) (SIGGRAPH ’06) . ACM, New York, NY, USA, 69–74. https://doi.org/10.1145/1185657.1185668[41] F. Fouss, A. Pirotte, J. Renders, and M. Saerens. 2007. Random-Walk Computation of Similarities between Nodes of aGraph with Application to Collaborative Recommendation.
IEEE Transactions on Knowledge & Data Engineering
19, 3(March 2007), 355–369. https://doi.org/10.1109/TKDE.2007.46[42] L. Gorelick, E. Sharon, R. Basri, A. Brandt, and M. Galun. 2006. Shape Representation and Classification Using thePoisson Equation.
IEEE Transactions on Pattern Analysis & Machine Intelligence
28 (12 2006), 1991–2005. https://doi.org/10.1109/TPAMI.2006.253[43] Karthik S. Gurumoorthy and Anand Rangarajan. 2009. A Schrödinger Equation for the Fast Computation of Approxi-mate Euclidean Distance Functions. In
Scale Space and Variational Methods in Computer Vision , Xue-Cheng Tai, KnutMørken, Marius Lysaker, and Knut-Andreas Lie (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 100–111., Vol. 1, No. 1, Article . Publication date: July 2020. [44] F. GÃűbel and A. A. Jagers. 1974. Random walks on graphs.
Stochastic Processes and their Applications
2, 4 (1974),311–336. https://EconPapers.repec.org/RePEc:eee:spapps:v:2:y:1974:i:4:p:311-336[45] Philipp Herholz and Marc Alexa. 2018. Factor Once: Reusing Cholesky Factorizations on Sub-Meshes.
ACM Transactionon Graphics (Proc. of Siggraph Asia)
37, 6 (2018), 9. https://doi.org/10.1145/3272127.3275107[46] Philipp Herholz, Timothy A Davis, and Marc Alexa. 2017. Localized solutions of sparse linear systems for geometryprocessing.
ACM Transactions on Graphics
36, 6 (2017).[47] Jin-ichi Itoh and Robert Sinclair. 2004. Thaw: A Tool for Approximating Cut Loci on a Triangulation of a Surface.
Experiment. Math.
13, 3 (2004), 309–325. https://projecteuclid.org:443/euclid.em/1103749839[48] Shreeraj Jadhav, Harsh Bhatia, Peer-Timo Bremer, Joshua A Levine, Luis Gustavo Nonato, and Valerio Pascucci. 2012.Consistent approximation of local flow behavior for 2d vector fields using edge maps. In
Topological Methods in DataAnalysis and Visualization II . Springer, 141–159.[49] Takashi Kanai and Hiromasa Suzuki. 2000. Approximate shortest path on a polyhedral surface based on selective re-finement of the discrete graph and its applications. In
Geometric Modeling and Processing 2000. Theory and Applications.Proceedings . IEEE, 241–250.[50] B KANEVA. 2000. An Implementation of Chen & Han’s Shortest Paths Algorithm. In .[51] R. Kimmel and J. A. Sethian. 1998. Computing Geodesic Paths on Manifolds. In
Proc. Natl. Acad. Sci. USA . 8431–8435.[52] D. J. Klein and M. Randić. 1993. Resistance distance.
Journal of Mathematical Chemistry
12, 1 (01 Dec 1993), 81–95.https://doi.org/10.1007/BF01164627[53] GVV Ravi Kumar, Prabha Srinivasan, V Devaraja Holla, KG Shastry, and BG Prakash. 2003. Geodesic curve computa-tions on surfaces.
Computer Aided Geometric Design
20, 2 (2003), 119–133.[54] Mark Lanthier. 1999.
Shortest path problems on polyhedral surfaces . Ph.D. Dissertation. Carleton University, School ofComputer Science.[55] Mark Lanthier, Anil Maheshwari, and Jörg-Rüdiger Sack. 1997. Approximating weighted shortest paths on polyhedralsurfaces. In
Proceedings of the 13th Annual ACM Symposium on COmputational Geometry . ACM, New York, New York,USA, 274–283.[56] Mark Lanthier, Anil Maheshwari, and Jörg-Rüdiger Sack. 2001. Approximating shortest paths on weighted polyhedralsurfaces.
Algorithmica
30, 4 (2001), 527–562.[57] G.F. Lawler. 2010.
Random Walk and the Heat Equation . American Mathematical Society. https://books.google.com/books?id=ujCIAwAAQBAJ[58] André Lieutier and Boris Thibert. 2009. Convergence of geodesics on triangulations.
Computer Aided GeometricDesign
26, 4 (2009), 412–424.[59] Yaron Lipman, Raif M Rustamov, and Thomas A Funkhouser. 2010. Biharmonic distance.
ACM Transactions onGraphics
29, 3 (June 2010), 1–11.[60] R. Litman and A. Bronstein. 2016. SpectroMeter: Amortized Sublinear Spectral Approximation of Distance on Graphs.
ArXiv e-prints (Sept. 2016). arXiv:1609.05715 [cs.DS][61] Bangquan Liu, Shuangmin Chen, Shi-Qing Xin, Ying He, Zhen Liu, and Jieyu Zhao. 2017. An optimization-drivenapproach for computing geodesic paths on triangle meshes.
Computer-Aided Design
90, Supplement C (2017), 105 –112. https://doi.org/10.1016/j.cad.2017.05.022 SI:SPM2017.[62] Yong-Jin Liu. 2013. Exact geodesic metric in 2-manifold triangle meshes using edge-based data structures.
Computer-Aided Design
45, 3 (2013), 695–704.[63] Yong-Jin Liu, Chun-Xu Xu, Dian Fan, and Ying He. 2015. Efficient Construction and Simplification of DelaunayMeshes.
ACM Trans. Graph.
34, 6, Article 174 (Oct. 2015), 13 pages. https://doi.org/10.1145/2816795.2818076[64] Yong-Jin Liu, Qian-Yi Zhou, and Shi-Min Hu. 2007. Handling degenerate cases in exact geodesic computation ontriangle meshes.
The Visual Computer
23, 9-11 (2007), 661–668.[65] Songting Luo. 2013. A uniformly second order fast sweeping method for eikonal equations.
J. Comput. Phys.
Discrete Mathematics
Journal of mechanicaldesign
SmartTools and Apps for Graphics - Eurographics Italian Chapter Conference . The Eurographics Association.[69] Dimas Martínez, Luiz Velho, and Paulo C Carvalho. 2005. Computing geodesics on triangular meshes.
Computers &Graphics
29, 5 (2005), 667–675.[70] Christian S Mata and Joseph S B Mitchell. 1997. A new algorithm for computing shortest paths in weighted planarsubdivisions (extended abstract). In
Proceedings Symposium on Computational Geometry . ACM, New York, New York,, Vol. 1, No. 1, Article . Publication date: July 2020.
Survey of Algorithms for Geodesic Paths and Distances 55
USA.[71] Joseph S. B. Mitchell, David M. Mount, and Christos H. Papadimitriou. 1987. The Discrete Geodesic Problem.
SIAM J.Comput.
16, 4 (Aug. 1987), 647–668. https://doi.org/10.1137/0216045[72] Thien Nguyen, Kestutis Karciauskas, and Jorg Peters. 2016. C Finite Elements on Non-Tensor-Product 2d and 3dManifolds. (2016).[73] Joseph O’Rourke, Subhash Suri, and Heather Booth. 1984. Shortest paths on polyhedral surfaces. In
STACS 85 ,K. Mehlhorn (Ed.). Springer Berlin Heidelberg, Berlin, Heidelberg, 243–254.[74] S. Osher and R. Fedkiw. 2003.
Level Set Methods and Dynamic Implicit Surfaces . Springer Verlag.[75] S. Osher and R. Fedkiw. 2003.
Level Set Methods and Dynamic Implicit Surfaces . Springer Verlag.[76] Giuseppe Patané. 2016. STAR - Laplacian spectral kernels and distances for geometry processing and shape analysis.In
Computer Graphics Forum . Consiglio Nazionale delle Ricerche, Rome, Italy, Wiley/Blackwell (10.1111), 559–624.[77] NM Patrikalakis and L Bardis. 1989. Offsets of curves on rational B-spline surfaces.
Engineering with Computers
5, 1(1989), 39–46.[78] Nicholas M Patrikalakis and Takashi Maekawa. 2009.
Shape interrogation for computer aided design and manufacturing .Springer Science & Business Media.[79] Gabriel Peyré and Laurent D Cohen. 2009. Geodesic methods for shape and surface processing. In
Advances inComputational Vision and Medical Image Processing . Springer, 29–56.[80] Gabriel Peyré, Mickael Péchaud, Renaud Keriven, Laurent D Cohen, et al. 2010. Geodesic methods in computer visionand graphics.
Foundations and Trends® in Computer Graphics and Vision
5, 3–4 (2010), 197–397.[81] Konrad Polthier and Markus Schmies. 1998. Straightest geodesics on polyhedral surfaces. In
Mathematical Visualization .Springer-Verlag, New York, 135–150.[82] Konrad Polthier and Markus Schmies. 1999. Geodesic flow on polyhedral surfaces. In
Data VisualizationâĂŹ99 .Springer, 179–188.[83] Yipeng Qin. 2017.
Fast and exact geodesic computation using Edge-based Windows Grouping.
Ph.D. Dissertation.Bournemouth University.[84] Yipeng Qin, Xiaoguang Han, Hongchuan Yu, Yizhou Yu, and Jianjun Zhang. 2016. Fast and Exact Discrete GeodesicComputation Based on Triangle-oriented Wavefront Propagation.
ACM Trans. Graph.
35, 4, Article 125 (July 2016),13 pages. https://doi.org/10.1145/2897824.2925930[85] Nicolas Ray and Dmitry Sokolov. 2014. Robust polylines tracing for n-symmetry direction field on triangulatedsurfaces.
ACM Transactions on Graphics (TOG)
33, 3 (2014), 30.[86] Raif M. Rustamov. 2007. Laplace-Beltrami Eigenfunctions for Deformation Invariant Shape Representation. In
Proceedings of the Fifth Eurographics Symposium on Geometry Processing (Barcelona, Spain) (SGP ’07) . EurographicsAssociation, Aire-la-Ville, Switzerland, Switzerland, 225–233. http://dl.acm.org/citation.cfm?id=1281991.1282022[87] Marco Saerens, Francois Fouss, Luh Yen, and Pierre Dupont. 2004. The Principal Components Analysis of a Graph,and Its Relationships to Spectral Clustering. In
Machine Learning: ECML 2004 , Jean-François Boulicaut, FlorianaEsposito, Fosca Giannotti, and Dino Pedreschi (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 371–383.[88] Shubhangi Saraf. 2009. Acute and nonobtuse triangulations of polyhedral surfaces.
European Journal of Combinatorics
30, 4 (2009), 833 – 840. https://doi.org/10.1016/j.ejc.2008.08.004[89] Olaf Schenk, Klaus Gärtner, Wolfgang Fichtner, and Andreas Stricker. 2001. PARDISO: A High-performance Serialand Parallel Sparse Linear Solver in Semiconductor Device Simulation.
Future Gener. Comput. Syst.
18, 1 (Sept. 2001),69–78. https://doi.org/10.1016/S0167-739X(00)00076-5[90] R Schmidt, C Grimm, and B Wyvill. 2006. Interactive decal compositing with discrete exponential maps.
ACMTransactions on Graphics (TOG)
25, 3 (2006), 605–613.[91] N. Sharp, Y. Soliman, and K. Crane. 2018. The Vector Heat Method.
ArXiv e-prints (May 2018). arXiv:1805.09170 [cs.GR][92] A. Sinha and M. Kazhdan. 2016. Geodesics using Waves: Computing Distances using Wave Propagation.
ArXive-prints (Dec. 2016). arXiv:1612.02509 [cs.CG][93] Vitaly Surazhsky, Tatiana Surazhsky, Danil Kirsanov, Steven J. Gortler, and Hugues Hoppe. 2005. Fast Exact andApproximate Geodesics on Meshes.
ACM Trans. Graph.
24, 3 (July 2005), 553–560. https://doi.org/10.1145/1073204.1073228[94] Z. Tari, J. Shah, and H. Pien. 1997. Extraction of Shape Skeletons from Grayscale Images.
Computer Vision and ImageUnderstanding
66, 2 (1997), 133–146.[95] Marc Troyanov. 1986. Les Surfaces Euclidiennes a Singularites Coniques.
Ens. Math
32 (1986), 74–94.[96] S. R. S. Varadhan. 1967. On the behavior of the fundamental solution of the heat equation with variable coefficients.
Communications on Pure and Applied Mathematics
20, 2 (1967), 431–455.[97] Christina Vaz, Keenan Crane, and Andreas Fabri. 2018. Heat Method with Intrinsic Delaunay Triangulation. In
CGALUser and Reference Manual (4.12.1 ed.). CGAL Editorial Board. https://doc.cgal.org/4.12.1/Manual/packages.html [98] Charlie CL Wang. 2004. CyberTape: an interactive measurement tool on polyhedral surface.
Computers & Graphics
28, 5 (2004), 731–745.[99] Xiaoning Wang, Zheng Fang, Jiajun Wu, Shi-Qing Xin, and Ying He. 2017. Discrete Geodesic Graph (DGG) forComputing Geodesic Distances on Polyhedral Surfaces.
Comput. Aided Geom. Des.
52, C (March 2017), 262–284.https://doi.org/10.1016/j.cagd.2017.03.010[100] Max Wardetzky, Saurabh Mathur, Felix Kaelberer, and Eitan Grinspun. 2007. Discrete Laplace operators: No freelunch. In
Geometry Processing , Alexander Belyaev and Michael Garland (Eds.). The Eurographics Association. https://doi.org/10.2312/SGP/SGP07/033-037[101] Qian Xie, Sebastian Kurtek, Huiling Le, and Anuj Srivastava. 2013. Parallel transport of deformations in shape spaceof elastic surfaces. In
Proceedings of the IEEE International Conference on Computer Vision . Florida State University,Tallahassee, United States, 865–872.[102] Shi-Qing Xin, Ying He, and Chi-Wing Fu. 2012. Efficiently computing exact geodesic loops within finite steps.
IEEETransactions on Visualization and Computer Graphics
18, 6 (2012), 879–889.[103] Shi-Qing Xin and Guo-Jin Wang. 2007. Efficiently Determining a Locally Exact Shortest Path on Polyhedral Surfaces.
Comput. Aided Des.
39, 12 (Dec. 2007), 1081–1090. https://doi.org/10.1016/j.cad.2007.08.001[104] Shi-Qing Xin and Guo-Jin Wang. 2009. Improving Chen and Han’s algorithm on the discrete geodesic problem.
ACMTransactions on Graphics (TOG)
28, 4 (2009), 104.[105] Shi-Qing Xin, Xiang Ying, and Ying He. 2012. Constant-time All-pairs Geodesic Distance Query on Triangle Meshes.In
Proceedings of the ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games (Costa Mesa, California) (I3D’12) . ACM, New York, NY, USA, 31–38. https://doi.org/10.1145/2159616.2159622[106] Chunxu Xu, Tuanfeng Y. Wang, Yong-Jin Liu, Ligang Liu, and Ying He. 2015. Fast Wavefront Propagation (FWP) forComputing Exact Geodesic Distances on Meshes.
IEEE Trans. Vis. Comput. Graph.
21, 7 (2015), 822–834.[107] Fang Yang and Laurent Cohen. 2015. Geodesic Distance and Curves Through Isotropic and Anisotropic Heat Equationson Images and Surfaces.
J Math Imaging Vis (2015).[108] Xiang Ying, Xiaoning Wang, and Ying He. 2013. Saddle Vertex Graph (SVG): A Novel Solution to the Discrete GeodesicProblem.
ACM Trans. Graph.
32, 6, Article 170 (Nov. 2013), 12 pages. https://doi.org/10.1145/2508363.2508379[109] Xiang Ying, Shi-Qing Xin, and Ying He. 2014. Parallel Chen-han (PCH) Algorithm for Discrete Geodesics.
ACM Trans.Graph.
33, 1, Article 9 (Feb. 2014), 11 pages. https://doi.org/10.1145/2534161[110] Carol T. Zamfirescu. 2013. Survey of two-dimensional acute triangulations.
Discrete Mathematics