Fast exact and approximate geodesics on meshes

The computation of geodesic paths and distances on triangle meshes is a common operation in many computer graphics applications. We present several practical algorithms for computing such geodesics from a source point to one or all other points efficiently. First, we describe an implementation of the exact "single source, all destination" algorithm presented by Mitchell, Mount, and Papadimitriou (MMP). We show that the algorithm runs much faster in practice than suggested by worst case analysis. Next, we extend the algorithm with a merging operation to obtain computationally efficient and accurate approximations with bounded error. Finally, to compute the shortest path between two given points, we use a lower-bound property of our approximate geodesic algorithm to efficiently prune the frontier of the MMP algorithm. thereby obtaining an exact solution even more quickly.


Introduction
In this paper we present practical methods for computing both exact and approximate shortest (i.e.geodesic) paths on a triangle mesh.These geodesic paths typically cut across faces in the mesh and are therefore not found by the traditional graph-based Dijkstra algorithm for shortest paths.The computation of geodesic paths is a common operation in many computer graphics applications.For example, parameterizing a mesh often involves cutting the mesh into one or more charts (e.g.[Krishnamurthy and Levoy 1996;Sander et al. 2003]), and the result generally has less distortion and better packing efficiency if the cuts are geodesic.Geodesic paths are used in segmenting a mesh into subparts, as done in [Katz and Tal 2003;Funkhouser et al. 2004].Mesh editing systems such as [Kobbelt et al. 1998] also use geodesics to delineate the extents of editing operations.Simulating fire on a mesh [Lee et al. 2001] also benefits from geodesics.In addition, geodesic paths establish a surface distance metric, which is an essential building block for many other techniques.For example, radial-basis interpolation over a mesh requires calculation of geodesic distances, and is used in numerous applications such as skinning [Sloan et al. 2001], mesh watermarking [Praun et al. 1999], and the definition of surface vector fields [Praun et al. 2000].Shape classification algorithms such as [Hilaga et al. 2001] use Morse analysis of a geodesic distance field.Parameterization metrics based on isomaps [Zigelman et al. 2002;Zhou et al. 2004;Peyré and Cohen 2005] are also driven by geodesic distances.In this paper we explore the problem of producing both exact and approximate solutions for geodesic paths (and hence distances) on triangle meshes (Figure 1).We present three contributions: Exact algorithm We first present an efficient implementation of the exact geodesic algorithm by Mitchell, Mount, and Papadimitriou (MMP) [1987].Using a simple parameterization of the dis- tance function over the edges, the implementation is actually practical even though, to our knowledge, it has never been done previously.We demonstrate that the algorithm's worst case running time of O(n 2 log n) is pessimistic, and that in practice, the algorithm runs in sub-quadratic time.For instance, we can compute the exact geodesic distance from a source point to all vertices of a 400K-triangle mesh in about one minute.Approximation algorithm We extend the algorithm with a merging operation to obtain computationally efficient and accurate approximations with bounded error.In practice, the algorithm runs in O(n log n) time even for small error thresholds.Exact geodesic path between two points We show how to efficiently obtain the exact solution to the "single source, single destination" problem, by using a lower-bound property of our approximation algorithm to prune the frontier of the MMP algorithm.In practice, we compute the shortest path between two points on a 1M-triangle mesh in just a few seconds.

Related work
The MMP algorithm [Mitchell et al. 1987] provides an exact solution for the "single source, all destination" shortest path problem on a triangle mesh.Their algorithm partitions each mesh edge into a set of intervals (windows) over which the exact distance computation can be performed atomically.These windows are propagated in a "continuous Dijkstra"-like manner.They prove a worst case running time of O(n 2 log n).Unfortunately, as far as we know the MMP algorithm has not been implemented previously and thus has not made its way into practice.An exact geodesic algorithm with worst case time complexity of O(n 2 ) was described by Chen and Han [1996] and partially implemented by Kaneva and O'Rourke [2000].We show that our MMP implementation runs many times faster than that implementation.Kapoor [1999] describes an algorithm for the "single source, single destination" geodesic path between two given mesh vertices, in O(n log 2 n) time.This is a complicated method which calls as subroutines many other complicated computational geometry algorithms; it is unclear if this algorithm will ever be realized.Approximate geodesics with guaranteed error bounds can be obtained by adding extra edges into the mesh and running Dijkstra on the one-skeleton of this augmented mesh [Lanthier et al. 1997].
Many extra edges are required to obtain accurate geodesics.Algorithms such as [Kanai and Suzuki 2001;Martinez et al. 2004] rely on iterative optimization, and the results therefore depend significantly on the initial approximate path.If this approximation is poor, the method might converge to an (incorrect) locally shortest path, or might require a large number of iterations.Mitchell [2000] presents a broad survey of approximate algorithms for graph and geodesic search.Kimmel and Sethian [1998] employ a variant of the fastmarching method to compute approximate geodesics on meshes in O(n log n) time.Several methods [Novotni and Klein 2002;Kirsanov 2004;Reimers 2004] explore an improved update rule for geodesic computation.Many of these methods require special processing of triangles with obtuse angles.We show that our MMPbased approximation algorithm yields more accurate solutions than the fast-marching method when applied to meshes.Moreover in Appendix B, we also demonstrate benefits of our algorithm when applied to meshes that approximate smooth manifolds.Polthier and Schmies [1998] explore a different definition of geodesic path on meshes using a notion of "straightest" instead of "shortest".Because these straightest geodesics are not always defined between pairs of points on a mesh, this notion may be inappropriate for many applications.Pham-Trong et al. [2001] explore geodesic paths over smooth parametric surfaces.

Exact algorithm
A window in the wall admits light into the room and its borders define the illuminated regions on the other walls of the room.

Common knowledge
Given a piecewise planar surface S defined by a triangle mesh, and a source vertex vs ∈ S, the MMP algorithm computes an explicit representation of the geodesic distance function D : S → R. For any point p ∈ S, this function D(p) returns the length of the geodesic path from p back to the source vs. Once a complete representation for D has been computed, one can quickly apply a "backtracing" algorithm to compute the shortest path from any query point to the source.This distance function can also be used to calculate isolines of constant distance (Figure 1).Shortest paths can be visualized as rays emanating from the source vertex vs in all tangent directions.These shortest paths are governed by the following three properties.Interior to a triangle, a shortest path must be a straight line.When crossing over an edge, a shortest path must correspond to a straight line if the two adjacent faces are unfolded into a common plane.As proven in [Mitchell et al. 1987], the only vertices that a shortest path can pass through (besides the source and destination) are (1) boundary vertices, (2) saddle vertices which are vertices with total angle greater than 2π (i.e. also called hyperbolic vertices), and (3) parabolic vertices whose total angle equals 2π.Parabolic vertices typically do not require special treatment since their neighborhoods unfold isometrically.Algorithm overview The basic idea of the MMP algorithm is to track together groups of shortest paths that can be parameterized atomically.This is achieved by partitioning each mesh edge into a set of intervals that we call windows.We will show that all the shortest paths within a window can be encoded locally using a 6tuple (b 0, b1, d0, d1, σ, τ ).The windows are then propagated across faces of the mesh in Dijkstra-like sweep.

Distance field along a window
Consider a specific shortest path from the source vertex vs to some point p on an edge e, and let us assume that this path does not pass through any saddle or boundary vertices.In this case, when all the faces intersecting the path are unfolded in a single common plane, the path forms a straight line.Consider the set of neighboring points on e whose shortest paths back to the source pass through the same sequence of faces.These paths are also straight lines in the same unfolding.In particular, the paths form a pencil of lines emanating from the unfolded source vertex.We represent this group of shortest paths atomically over a window w of the edge e (see Figure 2(a)).
The distance field D(p) over the window w is represented compactly as follows.We first store the endpoints of the window using two scalar valued parametric coordinates b0, b1 ∈ [0, e ] measuring distance along the edge.Next, we encode the position of the source vertex (relative to the window in the planar unfolding) using its distances d0, d1 to the window endpoints.Finally, we record a binary direction τ specifying the side of the edge on which the source lies.From this tuple (b0, b1, d0, d1, τ ), it is easy to position the source in the planar unfolding adjacent to the edge, by intersecting two circles as shown in Figure 2 σ Suppose now that the shortest path passes through one or more saddle/boundary vertices on its way to the source, and let s be the nearest such vertex to w.Again, consider the set of neighboring points on e whose shortest paths go through the same strip of faces back to s.In the unfolding of the strip between e and s, these shortest paths will form a pencil of lines emanating from s as seen in Figure 3.Because this set of shortest paths share the same path from s back to the source vertex vs, the distance field is now characterized by (1) the position of this pseudosource vertex s relative to the edge, and (2) the length σ = D(s) of the path from s back to the source vs.We refer to σ as the pseudosource distance.To summarize, the distance field D(p) over the window is expressed as a tuple (b0, b1, d0, d1, σ, τ ), where b0, b1 define the endpoints of w, d0, d1 are the corresponding distances to the pseudosource, σ is the geodesic distance from s to the source vs, and τ encodes the direction of s from the directed edge e.

Window propagation
Given a window w on an edge e1, we propagate its distance field across an adjacent face f to define new potential windows on the two opposing edges e2, e3.Essentially, we compute how the pencil of straight lines would extend across one more unfolded face in the strip.However, the edges e2 and e3 may already contain previously propagated windows, so we must "intersect" these previous windows with the new potential windows, to capture their combined minimum distance field.In other words, we only keep the part of a newly computed potential window that has a smaller geodesic distance than that already associated with points on e2 or e3.
Let w be the new potential window on one of the two opposing edges (Figure 4(a,b)).To define the distance field over w , we extend the rays from the pseudosource s through the endpoints of w and intersect them with the new edge, to obtain the new interval We then compute the new distances d 0 , d 1 from these new endpoints to the pseudosource s.The pseudosource distance σ = σ is unchanged, and the direction τ is assigned to point into face f .We have a special case when w is adjacent to a saddle/boundary vertex v, since shortest paths may pass through v, i.e. v may act as a new pseudosource.Suppose v lies at the left endpoint q0 = p0 of w.In this case we add extra windows on the edges − − → p1p2 and − − → p2p0.These additional windows cover the parts of the edges that lie to the left of the ray (s, p0) and are not already "illuminated" by s through w (Figure 4(c)).These new additional windows will have a pseudosource at v with σ = D(v).The case with v at q1 = p1 is treated in a symmetrical manner.

Intersection of overlapping windows
Suppose that the newly created window w0 and at least one existing window w1 on edge e have a nonempty intersection region δ = w0∩w1.We must decide which of the windows defines the minimal distance function for each point in δ, and update the windows along e appropriately.In particular, if one of the windows defines a larger distance everywhere over δ, then we simply cut δ away from its interval.The more interesting case is when w0 is minimal on part of δ, while w1 is minimal on the remaining part of δ.To correctly partition δ, we must then find the point p ∈ δ where the distance functions defined by w0 and w1 are equal, i.e. s0 − p + σ0 = s1 − p + σ1.If we define our planar coordinate system to align e with the x axis as shown in Figure 5(a), this can be expressed as: This equation can be reduced to a quadratic with a single solution in the required range, i.e. px ∈ δ such that Finally, we adjust the boundary b1 of the left window and b0 of the right window to the location of px as shown in Figure 5(b).

Continuous Dijkstra
The MMP algorithm propagates distance information out from the source in a Dijkstra-like fashion.When windows are created, they are placed in a priority queue sorted by distance back to the source.When a window is popped off the queue, it is propagated outward across a face.We next present this algorithm in more detail.
The priority queue is initialized with a window for each edge adjacent to the source vs.The distance fields for these initial edges are trivial, e.g.vertices adjacent to vs are assigned distances D equal to the corresponding edge lengths.The basic step is to select (and remove) a window from the queue and propagate it as in Section 3.2.Note that the propagation step can add, modify, or remove existing windows, and the queue is updated accordingly.We repeat the process until the queue is empty.Interestingly, the algorithm will generate the correct solution regardless of the order in which windows are removed from the queue.However, selecting windows in arbitrary order leads to an extremely slow result.The optimal approach is to propagate the windows as a wavefront, by ordering them in the queue according to minimal distance from the source vertex, i.e. minp∈w D(p) for window w.In practice, other reasonable criteria are acceptable.For example, the minimal distance at the window endpoints is faster to compute and the overall algorithm performance is nearly the same.Because of the limited precision of the computations, small gaps or overlaps can be generated between windows on an edge.We check and fix such problems by either extending or reducing the window extents.Even for our largest models (700K faces) we have not encountered numerical difficulties.Note that the total path lengths are simple sums of distances across finitely many triangle faces (O( √ n) in practice), so numerical errors do not accumulate exponentially as they might with ODE integration algorithms.

Geodesic path construction: Backtracing
Once all edges are covered by windows representing geodesic distance, it is easy to trace a shortest path from any surface point p back to the source.First, in the general case that p lies in a face interior, we consider all windows on the three edges bounding the face, and minimize p − p + D(p ) over all points p within these windows.Having jumped to a first window, we can then iteratively hop to previous windows all the way back to the source.That is, given location p on a current window, we find the adjacent face f according to the direction τ .We reconstruct the location of the pseudosource s in the plane of the face using the window parameters, and intersect the line − → p s with the two opposite edges of face f .This intersection point gives us a new point on a new window, from which we repeat the process.When reaching a pseudosource s itself (i.e. a saddle or boundary vertex), we iteratively move to windows adjacent to s, circumnavigating either clockwise or counter-clockwise, until finding a window with a new pseudosource.

Performance analysis
Let n be the number of the mesh edges.When the windows are propagated as a wavefront, it is shown in [Mitchell et al. 1987] that each edge may have O(n) windows and therefore the total number of windows can be O(n 2 ).This results in a worst case complexity of O(n 2 ) space and O(n 2 log n) time.The log n factor is due to the need for a priority queue and for the binary searching required to find overlapping windows during window propagation.However, for more typical meshes, we observe that edges have an average of O( √ n) windows.Intuitively, given a mesh with uniformly distributed vertices, the number of edges can be thought of as being proportional to area, while the number of edges crossed by a shortest path is proportional to diameter.This intuitive reasoning gives us the expected O( √ n) window-per-edge complexity.To confirm this intuition, we constructed a series of subdivision meshes for several simple surfaces using the Loop subdivision scheme [1987], and ran our implementation.Let Wi and Wi+1 be the total number of resulting windows on two subsequent levels  of mesh subdivision.We estimate the exponent p in the O(n p ) window complexity by evaluating log 4 (Wi+1/Wi), since the number of mesh edges increases by 4 at each level.From Figure 6 we see that the window complexity grows as approximately O(n 1.5 ).Surprisingly, the window complexity is even less when the mesh surface has a rough texture, as shown in the experiments of Figure 7. Intuitively, the bumpy features in the surface cause adjacent windows to overlap and thereby annihilate each other.

Approximation algorithm
In this section we introduce a method to compute an approximation D to the geodesic distance function D, which requires less time and space.The method works just like the exact algorithm, except for one key difference -before propagating a window, we try to merge it with an adjacent window on the same edge, as illustrated in Figure 8(a).The merging replaces two windows w0, w1 in the queue with a new window w that covers w0 ∪ w1.This involves choosing a new pseudosource s that effectively replaces the previous approximate distance function with a new one.There are constraints as described in the next section, and merging is only performed if these constraints can be satisfied.The windows w0 and w1 are then deleted from the queue, and w is inserted into the queue.
In addition, we prove in Appendix A that our approximation D is a lower bound, namely D(p) ≤ D(p), ∀p ∈ S.This property is employed later in Section 5.

Constraints on window merging
Some conditions must be satisfied before we merge two windows.Directionality: The two windows w0, w1 must have direction values τ in agreement.
Visibility: We define the visibility region of a window to be the area spanned by the rays exiting the window from the direction of the pseudosource, as illustrated in Figure 8(a).To guarantee that D is defined without gaps, the visibility region of the new window w must at least cover the union of those of the merged windows.Continuity: To maintain distance field continuity along the edges and at the vertices of the mesh, we must preserve distances at the endpoints of the merged window w .Accuracy: If window merging is performed indiscriminately, it usually results in slightly more than one window per edge on average, and is thus fast and consumes little memory.However, in practice it is desirable to bound the error max p=(x,0)∈w where D is the original and D the merged approximate distance function.The difference ∆D reaches its maximum either at the endpoints of wi, i = 0, 1, or where ∂ ∂x ∆D = 0, which is equivalent to the following quadratic equation for each of the two wi: To measure the error D − D between our approximation and the exact distance function, we accumulate error by storing a scalar value ξ on each window, assigning ξ = max(ξ0, ξ1) + ∆D(p).This is rather conservative but still works well in practice.Then, one possible test is to disallow the merge unless ξ ≤ ε abs where ε abs is an absolute error threshold.Instead, we prefer to bound the relative error (the error measured as a fraction of geodesic distance) at the same point p, by testing if ξ /D(p) ≤ ε rel .However, if we only bound the global (accumulated) error, many merges will occur near the source vertex vs (until the error threshold is reached), leaving little opportunity for later merges farther from the source, and thus possibly resulting in an excessive number of windows.Our solution is to additionally bound the local error x y 0 q 0 q 1 s s0 s1 x y ∆D(p) (actually, the relative local error ∆D(p)/D(p)) by a fraction of the global error tolerance.We set this fraction to be 10% in all our experiments, although ideally this value should depend on the size of the mesh.Thanks to this heuristic, we are able to satisfy the global error bound using significantly fewer windows overall.

Monotonicity:
We must be careful that the distance function over the edge is not smaller than the distance function over the parent windows from which it was propagated, where correspondence is defined by the propagation method.However, it is difficult to analytically check this property in practice.When we do not check this, loops in window propagation are possible.Because we maintain consistency of the source direction τ , our algorithm cannot produce any "bouncing" infinite loops, propagating back and forth between two adjacent edges.It is conceivable, though, that more exotic infinite loops could occur.Of course, we could explicitly maintain a directed graph representing the propagation evolution, explicitly check for loops in this graph, and disallow any attempted propagation steps that would complete a loop.In practice, we have never encountered such loops in any of our experiments, and so we do not explicitly perform these checks.

Finding the merged window pseudosource
We attempt to find a pseudosource position for the merged window that satisfies all the constraints of the previous section.Denote the two adjacent windows wi, with pseudosources si, for i = 0, 1.Using the same local coordinate frame used previously, the merged window w will have two endpoints q i ≡ (b i , 0).Our goal is to find a new pseudosource s = (s x , s y ) and pseudosource distance σ for the merged window that satisfy the following constraints.
To maintain continuity, we require that the geodesic distances at its endpoints Di ≡ D(q i ) = si − q i + σi are preserved by the merge.This can be expressed as It follows that This constrains s to lie on a conic curve γ.
To maintain directionality, we must impose the inequality s y ≥ 0.
To satisfy the visibility constraint we require our solution to lie inside the sector between the two lines Li, i = 0, 1 that pass through the q i to si (Figure 8 To constrain σ and the d i to be non-negative we add the inequalities σ ≥ 0, q i − s = Di − σ ≥ 0. It follows from (4.1) that these inequalities are equivalent to: If all the above constraints are not simultaneously satisfiable we disallow the merge.Otherwise we pick the s with minimal σ value.This must occur when one of our inequality constraints is "tight".

Backtracing
The geodesic path for our approximation algorithm is traced similarly to the algorithm in Section 3.5 but with one essential difference.When window w is the result of merging two original windows, its pseudosource position is different from the pseudosource positions of those original windows.If we were to trace back in the direction of the merged pseudosource, the resulting path would be different from that represented in the forward propagation, and its overall length might exceed the computed error bound.
Our approach is to obtain the original pre-merge pseudosource by maintaining a list of references to the windows that were successively merged into w (together with their endpoints).The average length of these lists is only about 2 in all our experiments.Another benefit of using the correct pseudosource is that we can trivially guarantee that the source will be reached without any loops.

Exact geodesic between two vertices
Our goal is to find the geodesic path between a source vertex vs and a target vertex vt on the mesh.Note that it is inefficient to run our exact algorithm on the entire mesh (or even until reaching vt).
In this section we present an algorithm that performs a sequence of pruned searches, exploiting progressively tighter lower and upper bounds on geodesic distance, so that the final, exact algorithm need only be run on a "thin" region surrounding the solution.
Our approach can be seen as a "continuous A* search", in that it adapts the traditional edge-based A* algorithm [Pohl 1971] to meshes.A similar pruning approach is also explored in [Floater et al. 2002]  This inequality is the core of the following algorithm: Step 1: Using Dijkstra search on edges only, compute an upperbound distance Ust(Dijkstra) by searching from vs until vt is reached.This step is made almost twice as fast using a bidirectional search, which runs two simultaneous Dijkstra searches from vs and vt until they both retire a common vertex.
Step 2: Start our approximation search (Section 4) from vt until vs is reached, which computes a lower-bound distance function Lt(•).
During the search, we use the inequality (5.1) to prune the search by only propagating windows that have at least one point p satisfying Lt(p) + p, vs ≤ Ust(Dijkstra), where •, • measures Euclidean distance (not on the mesh) and is therefore a trivial lower bound on Ds(p).
Step 3: Using the windows provided by the previous step, apply backtracing (Section 4.3) to form a path from vs back to vt.
The length of this path defines a tighter upper-bound distance Ust(backtrace).
Step 4: Start our exact search (Section 3) from vs until vt is reached, which computes exact distance D(vs, •).During the search, we again use the inequality (5.1) to propagate only windows that have at least one point p satisfying Ds(p) + Lt(p) ≤ Ust(backtrace).
Step 5: The geodesic distance between vs and vt has now been computed as Dst = Ds(vt).To obtain the geodesic path vt back to vs, apply backtracing (Section 3.5) using the windows provided by the previous step.
As future work, it would be useful to obtain even tighter A* bounds by precomputing distances to a set of landmark points on the mesh, as explored in [Goldberg and Harrelson 2005] for graph search.

Experimental results
We tested the algorithms on a Pentium M 1.6GHz PC with 1GB RAM.As shown in Table 1, our exact algorithm is useful even for large models.For instance, the exact geodesic distance from a source point to all vertices of the 400K-triangle David model is computed in 75 seconds.By comparison, the implementation of Chen and Han [1996] by Kaneva and O'Rourke [2000] runs successfully only on the 30K-triangle Buddha S model from our dataset, with a computation time of over 28 hours.In practice the main bottleneck of our exact algorithm is the memory space required to store all the windows.For 1GB memory, we are able to process a mesh of up to 700K faces.This space complexity constraint provides strong motivation for our approximate algorithm.With a 0.1% relative error bound, our approximate geodesic algorithm runs significantly faster and uses much less memory than the exact algorithm.Table 1 reports  Absolute errors are reported as percentages of the object diameter.
We compare the accuracy of our approximate algorithm with the fast marching (FM) algorithm of [Kimmel and Sethian 1998].Specifically, we used the FM implementation of Peyré and Cohen [2003;2005], and verified that two other recent FM implementations [Reimers 2004;Sifri et al. 2003] produced identical distance results and had similar speed.Table 1 shows that our approximation algorithm has similar running times to the FM algorithm, but more importantly it has significantly better accuracy.The error distribution graph of Figure 11 shows that our algorithm also has much better accuracy than the improved update rule of [Kirsanov 2004].

Path results
Figure 10 shows our point-to-point exact shortest path algorithm of Section 5 applied to a 1M-triangle Buddha XL model.It takes about 4 seconds to compute the path crossing half the model.Shortest paths between relatively closer vertices can be computed at interactive rates.For small models, paths between two arbitrary vertices can be computed in a matter of milliseconds.Our mesh geodesic algorithm also provides a practical solution to the problem of path planning in the presence of obstacles [Hershberger and Suri 1999].
In this setting, the mesh is a triangulation of a planar region, with holes in the mesh encoding the geometry of the obstacles.A simple example is shown in the inset figure.

Conclusions and future work
We have presented both exact and approximate algorithms for computing geodesic paths and distances on triangle meshes, and showed that these algorithms are fast enough to be practical on complicated meshes.For point-to-point geodesic computation, we developed a technique that uses a sequence of pruned searches with narrowing distance bounds to efficiently home in on the exact shortest path.The bottleneck of the exact geodesic algorithm is the memory requirement due to the large number of windows.One idea to remedy this problem is to instead store a "pseudosource direction vector" at each mesh vertex.This direction vector is sufficient to repeatedly trace the geodesic path back (through successive face unfoldings) to preceding pseudosources.With this new data structure, windows can be deleted once they are retired from the queue.Other future directions are to consider weighted or anisotropic distance metrics, and to generalize to higher-dimensional or smooth manifolds.

Acknowledgments
We thank Ron Kimmel, Martin Reimers, Michael Floater and Leonard McMillan for fruitful discussions.We also thank Gabriel Peyré, Martin Reimers and Oren Sifri for providing us with their implementations of the fast marching method.We are very grateful to the anonymous reviewers for their detailed constructive comments.This research was partially supported by the Research Council of Norway (BeMatA project).Vertices are sorted on the horizontal axis by their exact geodesic distance from the source.Flat-exact is the method of [Kirsanov 2004].

Figure 12:
The first row shows three of the approximating meshes.The second row shows errors for the meshes on seven subdivision levels.The table shows the errors and convergence rates.scheme, with vertices projected onto the sphere; the coarsest mesh has 18 vertices and the finest one has 64K vertices.
To analyze the algorithms we measure errors using the maximal absolute difference between G(v) and the distances D(v) provided by the algorithms, as well as the maximal and average relative differences |D(v) − G(v)|/G(v) between the distances.
In addition, we numerically compute the convergence rate R = log 2 (Ei−1/Ei) where Ei are errors at successive subdivision levels.Figure 12 compares our exact and approximate algorithms with the FM method.The results indicate that the exact algorithm has second-order convergence, which is in agreement with the results of [Floater 2005] for approximation of curves.
For the results in Figure 12, we maintained an average of 1.5 windows per edge (WPE) by manually adjusting the relative error bound ε rel across subdivision levels.Note that our approximation algorithm is adaptive with respect to the number of windows per edge, since it allocates more windows to "difficult" mesh regions.Table 2 demonstrates that when WPE increases, the convergence rate grows and approaches that of the exact algorithm.These results also indicate that setting WPE to a value as small as 1.5 already increases convergence rate substantially.

Figure 1 :
Figure 1: Geodesic paths from a source vertex, and isolines of the geodesic distance function.

Figure 2 :
Figure 2: (a) Pencil of rays from the source vertex s passing through window w over the strip of unfolded triangles.(b) The position of the source vs (or more generally pseudosource s) is parameterized relative to the window w that lies on edge − − → p 0 p 1 .
(b), and thereby recover the distance field within this window.

Figure 3 :
Figure 3: (a) Orthogonal projection of a 3D mesh near a red saddle vertex.Part of the edge in the upper triangle is not visible by rays from the source vertex vs. (b)Unfolding the triangles into the plane of the upper triangle reveals that the total angle is greater than 2π, resulting in two different "images" of vs in the unfolding.All shortest paths from vs to the red window w pass through the red pseudosource vertex s.

Figure 4 :
Figure 4: (a) Window propagation resulting in one window.(b) Window propagation resulting in two windows.(c) Special case of window propagation (saddle vertex at p 0 ); two additional (red) windows are added to the left of the ray (s, p 0 ).

Figure 5 :
Figure 5: (a) Two overlapping windows w 0 and w 1 (with unfolded pseudosources s 0 and s 1 ) and their intersection δ = [b 0 , b 1 ].(b) The formation of two disjoint windows; in the special case that σ 0 = σ 1 , then s 0 − p = s 1 − p = d as illustrated.

Figure 6 :Figure 7 :
Figure 6: These various subdivision meshes show how the resulting number of windows increases as the mesh is subdivided.WPE is the average number of windows per edge; 'Exp' is the exponent p when estimating O(n p ). Knot Time WPE Smooth 5.16 20.6 +10% noise 3.67 15.1 +20% noise 2.52 10.7 +40% noise 1.48 6.4 +80% noise 0.81 3.5 Figure 7: Performance of the exact algorithm when adding random geometrical noise to the smooth knot surface.The noise moves each vertex in a random direction by a uniform random distance from 0 to the indicated percentage of average edge length.The model and zoom-in shown have "20%" noise.

Figure 8 :
Figure 8: (a) Merging two windows into one such that the visibility region is not reduced.(b) Selecting the pseudosource for the new window.The visibility constraint V corresponds to the yellow region while the constraints of (4.3) define the pink region.
(b)).(If the intersection point of Li has positive y-coordinate, then the allowed region V is a triangle.Otherwise V is open.) although their scheme lacks true distance bounds.Denote by Pst the geodesic path between vs and vt.Let Ds(p) and Dt(p) be the geodesic distances from point p to vs and vt respectively, and let Dst = Ds(t) = Dt(s) be the length of Pst.Then, it is obvious that any point p on Pst satisfies Ds(p) + Dt(p) = Dst.If Ls(p) and Lt(p) are lower-bound functions of Ds(p) and Dt(p) respectively, and Ust is an upper-bound value for Dst, then any point p on Pst also satisfies Ls(p) + Lt(p) ≤ Ust.(5.1) both (1) the maximal absolute difference |D(v) − D(v)| between the approximate and exact distances, and (2) the average relative difference |D(v)−D(v)|/D(v).

Figure 10 :
Figure 10: Illustration of our exact algorithm for finding the shortest path between pairs of vertices on the Buddha.The colored regions correspond to the pruned searches of successive search steps, and the final exact paths are shown in red.Timings in seconds are listed for both the 1M-face Buddha XL and 200K-face Buddha M meshes.

Figure 11 :
Figure 11: Error distribution on the rockerarm model (40K vertices).Vertices are sorted on the horizontal axis by their exact geodesic distance from the source.Flat-exact is the method of[Kirsanov 2004].

Table 1 :
Comparison of our exact and approximation algorithms with fast-marching for the models of Figure 9. Times are in seconds; WPE indicates average number of windows per edge.Figure 9: Models used for the tests inTable 1.

Table 2 :
Convergence rates of our approximation and exact algorithms for the sphere model.See also Figure12.When the average number of windows per edge (WPE) increases, the convergence rate approaches the exact algorithm rate.(Missing cells denote values identical to the exact algorithm.)