arXiv:1710.01620v1 [cs.DS] 4 Oct 2017

Celestial Walk: A Terminating Oblivious Walk
for Convex Subdivisions
Wouter Kuijper
Nedap N.V.

Victor Ermolaev
Nedap N.V.

Olivier Devillers
Loria, Inria, CNRS, Université de Lorraine, France.
October 5, 2017
Abstract
We present a new oblivious walking strategy for convex subdivisions.
Our walk is faster than the straight walk and more generally applicable
than the visiblity walk. To prove termination of our walk we use a novel
monotonically decreasing distance measure.

1

Introduction

Point location in a convex subdivision is a classical problem of computational
geometry for which several data structures have been designed with good complexities in the worst case [11, 4, 15]. These intricate solutions are often unused
in favor of simpler algorithms based on traversal of the planar subdivisions using
neighborhood relations between faces, also known as walking algorithms [2, 3, 9].
These walking algorithms can also be used as a building block in randomized
data structures for point location [14, 6].

Figure 1: Celestial walk (obtuse angles marked by orange dots).

1

Figure 2: Straight (left) and visibility (right) walks
Amongst convex subdivisions, Delaunay triangulations received a lot of attention because of their practical importance. For Delaunay triangulations,
essentially two walking strategies are used: the straight walk and the visibility
walk [9]. The straight walk visits all faces crossed by a line segment between a
known face and the query point, while the visibility walk goes from a face to
another if the query point is on the side of the new face with respect to the
supporting line of the edge common to the two faces (cf. Figure 2).
The straight walk trivially terminates in the face containing the query point
and generalizes to any planar subdivision but with the inconvenience of not
being oblivious: the starting face of the walk must be remembered during the
whole walk.
The visibility walk is oblivious, but proving its termination requires the use
of particular properties of the Delaunay triangulation, and actually the visibility
walk may loop in other subdivisions [9] (cf. Figure 3).
Regarding
√ performance, both walks may visit all triangles in the worst case
and visit O( n) triangles when the points are evenly distributed [10, 8]. From
a practical point of view, the visibility walk is simpler to implement and a bit
faster in practice because it uses less predicates to walk from a triangle to its

Figure 3: Visibility walk may loop

2

neighbor [7].

Contribution
We propose celestial distance 1 as a new way to measure the proximity between
an edge of the subdivision and a query point. This distance measure allows
to design new walking strategies and prove their termination. By design these
strategies are oblivious: only the current edge of the current face determines
the edges by which the walk progresses to the next face.
Our main contribution is the celestial walk which is a refinement of the
well-known visibility walk that has the additional advantage of terminating not
just on Delaunay triangulations but on arbitrary convex subdivisions. This
is particularly useful for constrained and/or incremental meshing where the
conditions necessary for termination of the visibility walk can be locally and/or
temporarily violated.
Another important feature of the celestial walk is that, like the visibility
walk, it uses only orientation predicates to navigate the mesh. In practice,
checking an orientation predicate reduces to computing the sign of a second degree polynomial, the degree of such polynomial being a relevant measure of the
predicate complexity [13, 1]. As a consequence, it becomes relatively straightforward to implement the celestial walk in an efficient and robust manner.

2

Pre-requisites

Let G = (V, E, F ) be a planar straight line graph (PSLG) consisting of a set of
vertices V , a set of half–edges E, and a set of faces F .
We abstract away the borders of G by assuming that it tiles the entire real
plane. At the same time we rule out dense tessellations by assuming G is locally finite meaning the number of vertices (edges, faces) intersecting a given,
bounded area is always finite.
We assume G is given in half-edge representation. In particular we assume
the following atomic functions [5]:
origin : E → V which maps every half-edge to its start vertex,
target : E → V which maps every half-edge to its end vertex,
edge : F → E which maps every face to some edge on its perimeter,
face : E → F which maps every half-edge to its corresponding (left-hand-side)
face,
next : E → E which maps every half-edge to its successor half-edge in the counterclockwise winding order of the face perimeter,
1 The name celestial distance refers to the practice of celestial navigation where angular
distances between the celestial bodies and the horizon are used for navigation at sea.

3

twin : E → E which maps every half-edge to its twin half-edge running in the
opposite direction, i.e.: twin(twin(e)) = e, origin(e) = target(twin(e)) and
vice versa.
The point location problem in PSLGs can now be formulated as follows: given
some goal location p ∈ R2 and an initial half-edge einit ∈ E, find some goal edge
egoal ∈ E such that p ∈ face(egoal ) using only next(·) and twin(·) to get from one
half-edge to the next, i.e.: there must exist a finite path einit = e0 . . . en = egoal
such that for all 0 < i ≤ n it holds ei = next(ei−1 ) or ei = twin(ei−1 ).

3

Celestial Distance

One problem that we encounter when we try to use Euclidean distance as a
measure of progress for a walking algorithm is the fact that Euclidean distance
is not always strictly decreasing for every step in the walk. The latter means
that it is not possible to prove termination using Euclidean distance alone. For
this reason we define the following augmented distance measure on (half-)edges.
For a given point p ∈ R2 and a half-edge e ∈ E we define the celestial distance
of e to p as a pair cd(e, p) = [d, α] where d is the length of the line segment
from p to the closest point on e and α is the wide angle (α ≥ π2 ) between e and
this line segment or 0 in case d = 0. We now define a lexicographic order on
celestial distances as follows:
[d, α] < [d0 , α0 ] iff d < d0 ∨ (d = d0 ∧ α < α0 )
We illustrate this new distance measure in Figure 4. The figure shows a convex face and the resulting partition of the plane obtained by grouping together
points based on their closest edge. This leads to two superimposed Voronoi
partitions: one based on classical Euclidean distance and one based on celestial
distance. The partition based on celestial distance is a proper refinement of the
partition based on Euclidean distance in the sense that points which were ambiguous under Euclidean distance are now partitionable under celestial distance.
Note that the points in the corner-cones radiating outward from the vertices are
all ambiguous under Euclidean distance (since their closest point on the polygon
is the corner vertex which is shared by two edges) yet under celestial distance
they can be further partitioned. In particular, for the partition based on celestial distance the ambiguous corner-cones are split according to the angular
bisectors.

4

Abstract Walk

With all the pre-requisites and our celestial distance measure in place we are
in a position to present the walking algorithm properly. We first present an
abstract version of the algorithm and give a correctness proof for this abstract
version.

4

closest point inside
an edge of the polygon
angular
bisector

closest point is
a vertex of the polygon

Figure 4: Voronoi diagrams of a polygon. In blue the segment Voronoi diagram
of the open edges and the vertices, in purple the celestial Voronoi diagram of
the edges.
In the next section, we will give an efficient, concrete instantiation of this
abstract version where the computation of the celestial distances will be completely implicit. However, the correctness of the final version will rest on the
correctness proof of the abstract version as given in this section.
Algorithm 1 Abstract Walk
1:
2:
3:
4:
5:
6:
7:
8:
9:

if p strictly right of e then
e ← twin(e)
end if
E 0 ← {e}
while E 0 6= ∅ do
e ← select(E 0 )
E 0 ← {twin(e0 ) | e0 ∈ face(e) ∧ cd(e0 , p) < cd(e, p) ∧ p strictly right of e0 }
end while
return face(e)

The algorithm is rather simple: given a starting edge e and a query point p
to the left of e, we select an edge of the face(e) that improves the distance to p.
The abstract walk is formalized in Algorithm 1.
In lines 1-3 we bootstrap the invariant that the query point is always to
the left of the current edge. In line 4 we bootstrap the current set of successor
candidates. In line 5 we enter the main loop. The loop invariant ensures that
the loop will terminate as soon as there are no more suitable successor edges

5

p
e0

r

e
Figure 5: Illustration for limit case of Theorem 1.
that have lower celestial distance to the query point. In the proof below we
will see how this condition is sufficient to ensure that, at termination, it holds
p ∈ face(e). In line 6 we non-deterministically select one of the candidate
successor edges. In line 7 we compute the next set of candidate successor edges
which are all edges on the current face perimeter that have the query point on
the right and have smaller celestial distance to the query point than the current
edge.
Theorem 1. For any planar subdivision, starting edge, and query point Algorithm 1 terminates with p ∈ face(e).
Proof. For termination note that, due to local finiteness of the mesh, and
monotonicity of our distance measure, an infinitely descending chain of celestial distances is ruled out. It remains to prove that after termination it holds
p ∈ face(e). For this it would suffice to show that for a half-edge e and goal
location p such that p is to the left-of e but not in face(e) there always exists another e0 in the perimeter of face(e) such that p is strictly on the right
of e0 and cd(e0 , p) < cd(e, p). So let ep be the point on e closest to p and let
r be the ray from p to ep . We make a case distinction on r. First assume
that r intersects the interior of face(e). In this case it must hold that r intersects the boundary of face(e) at least one more time in another edge that
is closer to p than the current edge for the Euclidean distance, and thus for
our celestial distance (otherwise p would lie in the current face). Next assume
that r does not intersect the interior of face(e), In this case it must hold that
ep = origin(e) or ep = target(e). Assume, w.l.o.g., that ep = target(e). In
this case we claim that it holds that e0 = next(e) has strictly smaller celestial
distance to p. Assume, w.l.o.g., that e0 provides no improvement over e w.r.t.
Euclidean distance to p. It then immediately follows that the closest point on
e0 to p must be the only point that e and e0 share which, by assumption, is
ep = target(e) = origin(e0 ). Since r does not intersect the boundary of face(e) in
any point other than ep it must hold that e0 is in-between e and r (cf. Figure 5).
Now it holds that: wideangle(e, r) = angle(e, e0 ) + wideangle(e0 , r) which implies
wideangle(e0 , r) < wideangle(e, r). This establishes cd(e0 , p) < cd(e, p) and we
conclude the algorithm would progress to twin(e0 ).

6

5

Celestial Walk

In the previous section, the successor of an edge in the walk can be the twin of
any edge of the current face with a smaller distance to the query point. In this
section, we explain a way to actually select such an edge using only orientation
checks.
The main problems that we are solving in this section are the facts that Algorithm 1 is non-deterministic and relies on the explicit computation of celestial
distances. Both these properties make it less immediately applicable. In this
section we therefore develop a derived algorithm that works for convex PSLGs.
As we shall see, the additional assumption of convexity allows us a significantly
more efficient walk.
So let us first consider the problem of determining a successor edge that
has lower celestial distance than our current edge without having to explicitly
compute these distances.
As an example of a convex face consider once more Figure 4. If we assume
the query point is outside the face this leaves two possibilities. First, the query
point may be located inside one of the orthogonal slabs (indicated in blue in
Figure 4) which means the closest point coincides with the orthogonal projection
of the query point on the edge. Second, the query point may be located in one
of the intermediate corner-cones separating the orthogonal slabs which means
the closest point coincides with the corner vertex.
The latter argument gives us a basic refinement of Algorithm 1 for the convex
case: check if the query point is in one of the orthogonal slabs, if so, pick the
corresponding edge, else the query point must be inside some corner-cone inbetween two orthogonal slabs, in this case we can pick either edge (unless one of
them is our current edge in which case we are forced to pick the other alternative
in order to make progress).
To do even better we can resolve the remaining non-determinism by using
the angle as a tie-breaker to determine which of the two candidate edges is
best (i.e.: which edge is most favorably oriented towards the query point).
Theoretically the best tie breaker is the angular bisector (indicated in purple in
Figure 4). The only problem is that the explicit computation and representation
of the exact angular bisector is at least as hard as the explicit computation and
representation of celestial distances.
Fortunately we can avoid the explicit computation of angular bisectors as
well. In particular we make the following case distinction. If the corner vertex
is not obtuse (more precisely the face has an acute or right internal angle at
that corner vertex) it holds that the extensions of both edges lie inside the
corner-cone (i.e.: the leftmost corner in Figure 4). This means that we may
use either edge itself as a crude approximation to the angular bisector (this
is precisely what the visibility walk does in all cases). On the other hand, if
the corner vertex is obtuse it holds that the normal to the base of the triangle
spanned by the two candidate edges forms a suitable approximation to the
angular bisector (because it lies properly inside the corner-cone). In Figure 6
we illustrate this construction. Note that, in the limit, as the internal angle
7

e0

e00

Figure 6: Construction for approximating angular bisector.

e0

e00

Figure 7: Reducing obtuseness check to orientation check.
approaches 180 degrees, or, alternatively, the ratio between the lengths of the
two edges approaches 1, the base normal converges to the exact angular bisector.
The latter considerations lead us to define approx bisector(e0 , e00 ) for some
pair of consecutive edges e0 , e00 such that e00 = next(e0 ) and angle(e0 , e00 ) is
obtuse. In particular, we let approx bisector(e0 , e00 ) denote the line from target(e0 )
— which is equal by definition to origin(e00 ) — in the direction orthogonal to,
and to the right of, the internal baseline segment that connects origin(e0 ) and
target(e00 ). For an illustration of this construction cf. Figure 6. We also define
obtuse(e0 , e00 ) as a predicate that is true iff target(e00 ) is strictly right of the line
from target(e0 ) — which we assumed is equal to origin(e00 ) — in the direction
orthogonal to, and to the left of, e0 . For an illustration of this construction cf.
Figure 7. With these additional definitions Algorithm 2 formalizes the celestial
walk.
In lines 1-3 we bootstrap the invariant that the query point is always to the
left of the current edge. In line 4 we initialize a perimeter edge variable used
for iterating over the outline of the current face. In line 5 we enter the outer
loop. As can be seen by inspecting the code-paths inside the loop body, the
8

Algorithm 2 Celestial Walk
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

if p strictly right of e then
e ← twin(e)
end if
e0 ← next(e)
while e 6= e0 do
if p strictly right of e0 then
e00 ← next(e0 )
while obtuse(e0 , e00 ) and p left of approx bisector(e0 , e00 ) do
e0 , e00 ← e00 , next(e00 )
end while
e ← twin(e0 )
e0 ← next(e)
else
e0 ← next(e0 )
end if
end while
return face(e)

loop invariant ensures that the loop will terminate iff the query point is inside
the current face. In line 6 we check if the query point lies strictly to the right
of the perimeter edge. If this is the case we know, due to the face convexity,
that the query point is outside the face. From that moment on the only goal
is to find an edge that allows us to walk to a neighboring face. Therefore, in
line 7 we introduce a second perimeter edge variable that will be used to iterate
over pairs (e0 , e00 ) of consecutive edges rather than singletons. In line 8 we state
the negation of the condition that we are looking for: either the internal angle
of the face is acute2 or, otherwise, the query point falls strictly to the right of
the approximate bisector. In line 9 we shift the pair of edges along the face
perimeter until said condition is reached. In line 11 we have exited the inner
loop so we know that e0 contains the edge that is “sufficiently closest” to the
query point for us to continue the walk to the neighboring face. Hence we reassign e and e0 and drop back into the outer loop as though we would have
restarted the algorithm on this new edge.
The latter can be seen as a form of tail-recursion. Indeed, a recursive formulation is trivially possible by transforming the entire procedure into a function
and substituting a recursive call at line 12. We have not opted for the recursive
formulation because it would be slightly less efficient. The efficiency loss would
be due to the bootstrapping check in lines 1-3 which, in the recursion, would be
superfluous.
2 Note that the obtuseness check in line 8 can potentially be memorized by the mesh at the
expense of one bit per half-edge.

9

6

Conclusion

We have shown how our celestial distance measure allows to design novel walking
strategies. In particular we improved the applicability of the popular visibility
walk by refining it into what we call the celestial walk. Our walk is almost
as simple to implement yet is guaranteed to terminate on arbitrary convex
subdivisions.
If we assume that all obtuseness checks are pre-computed or memoized (at
the expense of 1 or 2 bits per half-edge respectively) we can make the following
observations w.r.t. the expected complexity of our walk.
The worst case for our walk materializes when all angles in the mesh are
obtuse, i.e. in a regular hexagonal mesh, in this case we require 2 orientation
tests per visited half-edge (the same number as the straight walk).
For triangulations a simple counting argument proves that there is at most
a fraction of 13 of obtuse angles. This gives us a worst case rough estimate
of 43 orientation tests per visited half-edge (better than the straight walk and
only slightly worse than the visibility walk). In practice, we expect better
performance since the fraction of obtuse angles will be much less than 31 in
triangulations that exhibit some regularity.
For meshes that are guaranteed to be triangular, minor optimizations are
possible by unrolling the loops of the general algorithm. Furthermore, it is a
good heuristic to minimize the number of visited triangles by alternating clockwise and counterclockwise winding order in order to balance out the systematic
errors in approximating the angular bisectors.

Future Work
An obvious next question to consider is whether the walk can be generalized to
3 dimensions. We have some reason to believe this is the case.
In 2D as we progress according to the celestial walk we either move forward
in the direction of the query point or orient ourselves towards it. Similarly,
in 3 dimensions we expect to either progress in Euclidian sense or roll in the
direction of the query point. This intuition leads to developing a suitable nonincreasing distance measure quantifying this relationship. In particular to the
walking problem in tetrahedralization we find that a combination of Euclidian
distance as the first progression criterion and two successive angles are sufficient to guarantee the progress towards the tetrahedron containing the query
point [12].
Finally, we would also be interested in obtaining experimental results regarding the trade-off that exists between steering (i.e.: choosing a next edge
that is favorably oriented to the query point) and driving (i.e.: choosing a next
edge quickly based on easily computed criteria).

10

References
[1] J.-D. Boissonnat and F. P. Preparata. Robust plane sweep for intersecting segments. SIAM Journal on Computing, 29(5):1401–1421, 2000.
[2] P. Bose, A. Brodnik, S. Carlsson, E. D. Demaine, R. Fleischer, A. López-Ortiz,
P. Morin, and J. Ian Munro. Online routing in convex subdivisions. International
Journal of Computational Geometry & Applications, 12(04):283–295, 2002.
[3] P. Bose and P. Morin. Online routing in triangulations. SIAM journal on computing, 33(4):937–951, 2004.
[4] B. Chazelle and L. J. Guibas. Fractional cascading: I. a data structuring technique. Algorithmica, 1(1):133–162, 1986.
[5] M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf. Computational
Geometry: Algorithms and Applications. Springer-Verlag, Berlin, 1997.
[6] O. Devillers. The Delaunay Hierarchy. International Journal of Foundations of
Computer Science, 13:163–180, 2002.
[7] O. Devillers. Delaunay triangulations, theory vs practice. In EuroCG, 28th European Workshop on Computational Geometry, Assisi, Italy, 2012.
[8] O. Devillers and R. Hemsley. The worst visibility walk in a random Delaunay
√
triangulation is O( n) . Journal of Computational Geometry, 7(1):332–359, 2016.
[9] O. Devillers, S. Pion, and M. Teillaud. Walking in a Triangulation. International
Journal of Foundations of Computer Science, 13:181–199, 2002.
[10] L. Devroye, C. Lemaire, and J.-M. Moreau. Expected time analysis for delaunay
point location. Computational geometry, 29(2):61–89, 2004.
[11] D. Kirkpatrick. Optimal search in planar subdivisions. SIAM Journal on Computing, 12(1):28–35, 1983.
[12] W. Kuijper, V. Ermolaev, and H. Berntsen. Zig-zagging in a triangulation. arXiv
preprint arXiv:1705.03950, 2017.
[13] G. Liotta, F. P. Preparata, and R. Tamassia. Robust proximity queries: An
illustration of degree-driven algorithm design. SIAM Journal on Computing,
28(3):864–889, 1998.
[14] E. P. Mücke, I. Saias, and B. Zhu. Fast randomized point location without preprocessing in two-and three-dimensional Delaunay triangulations. Computational
Geometry, 12(1-2):63–83, 1999.
[15] F. P. Preparata. Planar point location revisited. International Journal of Foundations of Computer Science, 1(01):71–86, 1990.

11

