AbstractThis note proposes a simple rule to determine a unique triangulation among all Delaunay triangulations of a planar point set, based on two preferred directions. We show that the triangulation can be generated by extending Lawson's edge-swapping algorithm and that point deletion is a local procedure. The rule can be implemented exactly when the points have integer coordinates and can be used to improve image compression methods

Dyken, Christopher

Floater, Michael S.

Elsevier - Publisher Connector 

Preferred directions for resolving the non-uniqueness of Delaunay triangulations 

Computational Geometry 34 (2006) 96–101
www.elsevier.com/locate/comgeo
Preferred directions for resolving the non-uniqueness
of Delaunay triangulations
Christopher Dyken ∗, Michael S. Floater
Centre of Mathematics for Applications, Department of Informatics, University of Oslo, Postbox 1053, Blindern, 0316 Oslo, Norway
Received 3 September 2004; received in revised form 6 February 2005; accepted 31 August 2005
Available online 25 October 2005
Communicated by J.D. Boissonnat
Abstract
This note proposes a simple rule to determine a unique triangulation among all Delaunay triangulations of a planar point set,
based on two preferred directions. We show that the triangulation can be generated by extending Lawson’s edge-swapping algo-
rithm and that point deletion is a local procedure. The rule can be implemented exactly when the points have integer coordinates
and can be used to improve image compression methods.
 2005 Elsevier B.V. All rights reserved.
Keywords: Delaunay triangulation; Uniqueness; Edge-swapping; Integer arithmetic
1. Introduction
Delaunay triangulations [3] play an important role in computational geometry [9]. A recent application that has
emerged is their use in compressing digital images [4]. Such images are represented by rectangular arrays of grey scale
values or colour values and one approach to compression is to start by representing them as piecewise linear functions
over regular triangulations and then to approximate these functions by piecewise linear functions over triangulations
of subsets of the points. If one could agree on a unique method of triangulating the points, one would obtain higher
compression rates because the sender would only need to encode the points and the height values, not the connectivity
of the triangulation: the receiver would be able to reproduce the triangulation exactly.
One advantage of Delaunay triangulations is that they are almost unique. In fact they are unique for point sets
containing no sets of four co-circular points. However, in the case that a set of planar points is a subset of a rectangular
array of points, there will typically be many co-circular points and therefore a large number of Delaunay triangulations.
An obvious and simple approach to the non-uniqueness problem is to perturb the points randomly before trian-
gulating but there is no guarantee that the perturbed points will have a unique Delaunay triangulation. An alternative
approach could be a symbolic perturbation method as discussed in [1,5,6,8,10]. Such a method should lead to a unique
Delaunay triangulation of the ‘perturbed’ points, but may not be a valid triangulation of the original ones.
* Corresponding author.
E-mail addresses: dyken@cma.uio.no (C. Dyken), michaelf@ifi.uio.no (M.S. Floater).0925-7721/$ – see front matter  2005 Elsevier B.V. All rights reserved.
doi:10.1016/j.comgeo.2005.08.004
C. Dyken, M.S. Floater / Computational Geometry 34 (2006) 96–101 97Fig. 1. A strictly convex quadrilateral.
The purpose of this note is to point out that non-uniqueness can be resolved without perturbing the points. We show
that a simple rule based on two preferred directions can be used to determine a unique member of all the Delaunay
triangulations of a set of points in the plane. We show that the rule can simply be incorporated into Lawson’s swapping
algorithm and that point deletion is a local procedure. We further show that, importantly, the rule can be computed
exactly in integer arithmetic. The preferred direction method could immediately be applied to the compression algo-
rithms described above.
2. Triangulating quadrilaterals
Let Q be a quadrilateral in the plane, with ordered vertices v1, w1, v2, w2 as in Fig. 1. If Q is strictly convex,
i.e., convex and such that no three vertices are collinear, then there are two ways to triangulate it, either by placing
a diagonal edge between v1 and v2, as in Fig. 1, or between w1 and w2. We want to propose a simple rule which
determines the diagonal of any such quadrilateral uniquely. A natural rule seems to be to use a preferred direction: we
choose that diagonal which makes the smallest angle with some arbitrary, fixed straight line. However, the two angles
could be equal, and in order to distinguish the two diagonals in this case, we will use a second preferred direction.
Thus we choose any two non-zero vectors d1 and d2 which are neither parallel nor orthogonal to each other. One such
choice would be d1 = (1,0) and d2 = (1,1). For any line (or line segment) , let αi(), where 0 αi() π/2, be the
angle between  and the (undirected) vector di , i = 1,2. Then we define the score of  as the ordered pair of angles
score() := (α1(),α2()
)
.
We compare scores lexicographically. Thus for two arbitrary lines  and m, we say that score() < score(m) if either
α1() < α1(m) or α1() = α1(m) and α2() < α2(m).
Lemma 2.1. Two lines  and m have the same score if and only if they are parallel.
Proof. If  and m are parallel, then clearly αi() = αi(m) for both i = 1,2 and so  and m have the same score.
Conversely, suppose  and m are not parallel but that they have the same score. Then at the point p of intersection
between  and m, the two lines through p in the directions d1 and d2 must bisect the lines  and m. But this can only
occur if d1 is either parallel or orthogonal to d2, which is a contradiction. 
Our preferred direction rule for the quadrilateral Q in Fig. 1 simply chooses that diagonal, [v1, v2] or [w1,w2],
with the lowest score. Since the two diagonals are never parallel, they always have distinct scores by Lemma 2.1.
3. Triangulating convex polygons
We next use the preferred direction rule to triangulate uniquely any strictly convex polygon P , illustrated in Fig. 2.
Consider EI (P ), the set of all interior edges of P , i.e., all line segments [v1, v2] connecting non-neighbouring pairs
of vertices v1 and v2 of P . We start by ranking all the edges of EI (P ) according to their score. We then employ an
insertion algorithm, inserting edges of EI (P ) into P in order of their ranking. In the first step, we insert all edges of
EI (P ) (one or more) which share the lowest score. Note that if there are more than one of these, they must all be
parallel by Lemma 2.1, and so they do not cross each other. In the general step, we insert all edges of EI (P ) with the
98 C. Dyken, M.S. Floater / Computational Geometry 34 (2006) 96–101Fig. 2. A strictly convex polygon and triangulation.
Fig. 3. The edges e and e′.
current lowest score which do not cross edges previously inserted. We continue until we have triangulated P , denoting
the triangulation by Td1,d2(P ).
We will show that Td1,d2(P ) has the useful property that it can be generated by a local optimisation procedure based
on edge swapping. Suppose T (P ) is any triangulation of P . Every interior edge e of T (P ) is the diagonal of a strictly
convex quadrilateral Q, and can therefore be swapped by the other diagonal e′ of Q to form a new triangulation. We
say that e is locally optimal if its score is lower than that of e′. Otherwise, we optimize Q by swapping e with e′. We
keep applying this local optimisation procedure until no more swaps can be performed, in which case we say that the
triangulation is locally optimal.
Lemma 3.1. Let T (P ) be any triangulation of P . If there is an edge e ∈ EI (P ) which crosses an edge of T (P ) with
a higher score, then T (P ) is not locally optimal.
Proof. Let e′ = [v1, v2] be an edge of T (P ) with the highest score among all edges of T (P ) which cross e. We will
show that e′ is not locally optimal. Let p be the intersection of e and e′. Without loss of generality, and by translating
and rotating P and d1 and d2 about p, we may assume that p = (0,0) and d1 = (1,0). Then since e′ has a higher
score than e, e′ cannot lie along the x-axis, and so we can assume that y1 < 0 < y2 (where vi = (xi, yi)). Further, by
reflecting all points, edges, and d1 and d2 about the y-axis if necessary, we may assume that x1  0 x2; see Fig. 3.
Next let ′ be the infinite straight line passing through e′ and let ′′ be its reflection about the y-axis, and let A1
and A2 be the two open semi-infinite cones bounded by ′ and ′′, with A1 in the positive y half-plane and A2 in the
negative y half-plane; see Fig. 4. Clearly because α1(e)  α1(′) = α1(′′), e does not intersect A1 ∪ A2. Thus e is
contained in the union of the two closed regions B1 and B2 shown in Fig. 4. Note, moreover, that e is not contained
in the line ′. It may however be contained in the line ′′.
Let Q be the convex quadrilateral of T (P ) with e′ as its diagonal. Let the two vertices in Q other than v1 and v2
be w1 and w2, with w1 lying on the side of ′ containing the positive x-axis and w2 lying on the side containing the
negative x-axis; see Fig. 5 for an example. Next we show that w2 ∈ B1. Indeed if w2 ∈ A1 then by the convexity of P ,
the edge [v1,w2] of T (P ) intersects e and [v1,w2] has a higher score than e′ which contradicts the definition of e′.
Furthermore, if w2 ∈ ′′ and e ⊂ ′′, then again the edge [v1,w2] would intersect e and have a higher score than e′.
Therefore, if w2 ∈ ′′ then also e ⊂ ′′. In this case, since no three vertices of P are collinear, we conclude that w2 is
an end point of e.
C. Dyken, M.S. Floater / Computational Geometry 34 (2006) 96–101 99Fig. 4. Regions of vertices.
Fig. 5. A possible quadrilateral.
Similarly, w1 ∈ B2, because if w1 ∈ A2 then the edge [v2,w1] would intersect e and have a higher score than e′
which contradicts the definition of e′. A similar argument shows that if w1 ∈ ′′ then e ⊂ ′′ and w1 must be an end
point of e.
To complete the proof, we will show that e′ is not locally optimal by showing that [w1,w2] has a lower score
than e′. First observe from Fig. 5 that α1([w1,w2])  α1(e′). Moreover, if either w2 is in the interior of B1 or w1
is in the interior of B2, or both, then the inequality is strict and so score([w1,w2]) < score(e′). The only remaining
possibility is that both w1 and w2 lie on the line ′′, in which case α1([w1,w2]) = α1(e′). We have shown, however,
that in this case both w1 and w2 are the end points of the edge e, so that e = [w1,w2]. Since e has a lower score than
e′ by assumption, e′ is therefore not locally optimal. 
An immediate consequence of Lemma 3.1 is
Theorem 3.2. If T (P ) is locally optimal then T (P ) = Td1,d2(P ).
Proof. Lemma 3.1 implies that if e is an interior edge of T (P ) then the only edges of EI (P ) which cross it have
a higher score than e. This means that when the score of e is reached in the insertion algorithm, e will be inserted
into P . Thus e is in Td1,d2(P ). 
Lemma 3.1 and Theorem 3.2 together clearly imply
Theorem 3.3. A line segment e ∈ EI (P ) is an edge of Td1,d2(P ) if and only if it is not crossed by an edge in EI (P )
with a lower score.
4. Delaunay triangulations
Finally we use the preferred direction rule to determine a unique triangulation among all possible Delaunay trian-
gulations of a given set of points in the plane.
100 C. Dyken, M.S. Floater / Computational Geometry 34 (2006) 96–101Any set of planar points V which are not all collinear admits a unique Delaunay pretriangulation, which is a tiling
of the points, whose boundary is the convex hull of V [11]. Two points form an edge in the tiling if and only if
they are strong neighbours in the Voronoi diagram of V [9,11]. Two points are Voronoi neighbours if their Voronoi
tiles intersect. Two such tiles intersect either in a line segment or a point. If the intersection is a line segment the
two points are strong neighbours, and they are weak neighbours otherwise. The vertices of each tile in the Delaunay
pretriangulation lie on a circle. Tiles with three vertices are triangles. By triangulating any tile with four or more
vertices arbitrarily, we convert the Delaunay pretriangulation into a Delaunay triangulation.
Since each tile P of the Delaunay pretriangulation is strictly convex, we can simply triangulate it using Td1,d2(P )
and in this way we determine a unique Delaunay triangulation of V , which we will denote by Td1,d2(V ).
Similarly to the case of convex polygons, we next show that given an arbitrary triangulation T (V ), we always
reach Td1,d2(V ) by edge swapping. We simply augment Lawson’s local optimisation procedure for Delaunay triangu-
lations [7] with the preferred direction rule as follows. Suppose e = [v1, v2] is an interior edge of some triangulation
T (V ). If Q, the quadrilateral having e as its diagonal, is not strictly convex we say that e is locally optimal. Other-
wise, if e′ = [w1,w2] denotes the opposite diagonal of Q, we say that e is locally optimal if w2 lies strictly outside
the circumcircle C through v1, v2, w1. If w2 lies strictly inside C, we swap e with e′. Otherwise v1, v2,w1,w2 are co-
circular and we use the preferred direction rule as a tie-breaker: e is locally optimal if score(e) < score(e′). Otherwise
we swap e with e′. We say that T (V ) is locally optimal if all its interior edges are locally optimal.
Theorem 4.1. If T (V ) is locally optimal then T (V ) = Td1,d2(V ).
Proof. If T (V ) is locally optimal then it is a Delaunay triangulation [7,11]. Thus every edge of T (V ) which is in the
Delaunay pretriangulation of V is also in Td1,d2(V ). Every remaining edge e of T (V ) is an interior edge of some tile
P of the pretriangulation. But since e is then in Td1,d2(P ) it must also be in Td1,d2(V ). 
Note further that an important operation on Delaunay triangulations is point deletion. It is well known that when
an interior vertex of a Delaunay triangulation is deleted, a Delaunay triangulation of the remaining points can be
constructed by simply retriangulating the hole created by the removal of v, and thus updating the triangulation is a
local operation, and can be implemented efficiently. We now show that the triangulation Td1,d2(V ) has an analogous
property.
Theorem 4.2. Let v ∈ V be an interior vertex of the triangulation Td1,d2(V ). Then every edge e of Td1,d2(V ) which is
not incident on v also belongs to the triangulation Td1,d2(V \ v).
Proof. Suppose first that e is an edge of the Delaunay pretriangulation of V . Then its end points are strong neighbours
in the Voronoi diagram of V , and are therefore also strong neighbours of the Voronoi diagram of V \ v, and so e is
also an edge of the Delaunay pretriangulation of V \ v, and therefore an edge of Td1,d2(V \ v).
The remaining possibility is that e is an interior edge of some strictly convex tile P of the Delaunay pretriangulation
of V , and all vertices of P lie on a circle. If v is not a vertex of P , then P will also be a tile in the Delaunay
pretriangulation of V \ v, and so e, being an interior edge of Td1,d2(P ), will be contained in Td1,d2(V \ v). Otherwise
v is a vertex of P . Then the polygon P ′, formed by removing v from P , is a tile in the Delaunay pretriangulation of
V \ v. Now if v and the two end points of e are consecutive vertices of P , then e will be an edge of P ′, and therefore
of Td1,d2(V \ v). Otherwise, since EI (P ′) ⊂ EI (P ), we see by Theorem 3.3 that since e is not crossed by any edge in
EI (P ) with a lower score, it is not crossed by any edge of EI (P ′) with a lower score and is therefore an interior edge
of Td1,d2(P ′), and hence also an edge of Td1,d2(V \ v). 
5. Numerical implementation
Using the implementation of the circumcircle test proposed by Cline and Renka [2], one can construct a Delaunay
triangulation of a set V of points with integer coordinates using integer arithmetic. So in order to construct the trian-
gulation Td1,d2(V ) in this case, we only need to show that the preferred direction rule can be implemented exactly. To
see that this is possible, consider the angle test
αi
([v1, v2]
)
< αi
([w1,w2]
)
, i = 1,2.
C. Dyken, M.S. Floater / Computational Geometry 34 (2006) 96–101 101Fig. 6. Triangulations with d1 = (1,0), d2 = (1,1) and d1 = (1,−1), d2 = (0,1).
One way to convert this to an integer comparison is to observe that it is equivalent to
cosαi
([v1, v2]
)
> cosαi
([w1,w2]
)
which, using scalar products, is equivalent to
∣∣∣∣
(v2 − v1)
|v2 − v1| ·
di
|di |
∣∣∣∣ >
∣∣∣∣
(w2 − w1)
|w2 − w1| ·
di
|di |
∣∣∣∣.
By squaring and removing the common denominator we get the equivalent test
|w2 − w1|2
(
(v2 − v1) · di
)2 − |v2 − v1|2
(
(w2 − w1) · di
)2
> 0.
Since the left hand side is a polynomial in the coordinates of the points v1, w1, v2, w2, and the vector di , we see that
provided d1 and d2 also have integer coordinates, the left hand side is also an integer. Clearly, all angle tests involved
in comparing two scores involve testing the sign (positive, zero, or negative) of such an integer.
Fig. 6 shows the Delaunay triangulation Td1,d2(V ) for two different choices of d1 and d2, where V is a subset of
points on a square grid. The triangulations in the figure were found by recursive point insertion, applying Lawson’s
swapping algorithm augmented with the preferred direction rule, and using integer arithmetic. One could, however,
use any method to find one of the possible Delaunay triangulations and then apply edge swapping with the preferred
direction rule to reach the triangulation Td1,d2(V ). Since the edges of the Delaunay pretriangulation will never be
swapped, the only swapping required will occur inside the tiles of the pretriangulation, and so the number of swaps
will generally be low. The Delaunay triangulations computed in [4] were found by recursive point removal and this
algorithm could be augmented with the preferred direction rule and the updates would be local due to Theorem 4.2.
References
[1] P. Alliez, O. Devillers, J. Snoeyink, Removing degeneracies by perturbing the problem or the world, Reliable Computing 6 (2000) 61–79.
[2] A.K. Cline, R.J. Renka, A storage efficient method for construction of a Thiessen triangulation, Rocky Mountain J. Math. 14 (1984) 119–140.
[3] B. Delaunay, Sur la sphère vide, Bull. Acad. Sci. USSR (VII), Classe Sci. Mat. Nat. (1934) 793–800.
[4] L. Demaret, N. Dyn, M.S. Floater, A. Iske, Adaptive thinning for terrain modelling and image compression, in: N.A. Dodgson, M.S. Floater,
M.A. Sabin (Eds.), Advances in Multiresolution for Geometric Modelling, Springer, Heidelberg, 2004, pp. 319–338.
[5] O. Devillers, M. Teillaud, Perturbations and vertex removal in a 3D Delaunay triangulation, in: SODA 2003, pp. 313–319.
[6] H. Edelsbrunner, E.P. Mücke, Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms, ACM Trans.
Graphics 9 (1990) 67–104.
[7] C.F. Lawson, Software for C1 surface interpolation, in: J. Rice (Ed.), Mathematical Software III, Academic Press, 1977, pp. 161–194.
[8] E.P. Mücke, A robust implementation for three-dimensional Delaunay triangulations, Int. J. Comput. Geom. Appl. 8 (1998) 255–276.
[9] J. O’Rourke, Computational Geometry in C, second ed., Cambridge University Press, Cambridge, 1998.
[10] R. Seidel, The nature and meaning of perturbations in geometric computing, Discrete Comput. Geom. 19 (1998) 1–17.
[11] R. Sibson, Locally equiangular triangulations, Comput. J. 21 (1978) 243–245.


file:///data/core-remote/dit/data/elsevier/pdf/881/aHR0cDovL2FwaS5lbHNldmllci5jb20vY29udGVudC9hcnRpY2xlL3BpaS9zMDkyNTc3MjEwNTAwMDgxNw==.pdf

Preferred directions for resolving the non-uniqueness of Delaunay triangulations

Abstract

Similar works

Full text

Available Versions

Elsevier - Publisher Connector

Elsevier - Publisher Connector