10. Delaunay Triangulation
Motivation. Let us consider a function f : R2
→ R. Its values are known only in the
final set of points P. One way to approximate such a function is to divide the convex
hull of the set P into triangles with vertices at the points of the set P and replace the
function f by a function which is linear on these triangles. This means the following:
Each point x of a triangle pi, pj, pk is a convex combination of the vertices
x = t1p1 + t2p2 + t3p3, where t1 + t2 + t3 = 1, t1 ≥ 0, t2 ≥ 0, t3 ≥ 0,
and hence we put
f(x) = t1f(p1) + t2f(p2) + t3f(p3).
The graph of the approximation of the function f is thus formed by triangles in R2
×R
with vertices (pi, f(pi), (pj, f(pj) and (pk, f(pk)).
FIGURE 10.1 Linear approximation of the function f.
This approximation depends substantially on how we do the triangulation of the
convex hull of our set. It seems that a good approximation should be one that has
no triangles with small angles. We will deal with such triangulations in the following
text.
Angle-optimal triangulation. First, we show that for a given set of points P in a
plane, all the triangulations of its convex hull have the same number of triangles. This
allows us to compare different triangulations lexicographically according to the sizes
of angles.
Theorem 10.1. Let P be a set of n points in the plane. Let the convex hull of P
have k edges. Then any triangulation of the convex hull of the set P has 2n − 2 − k
triangles and 3n − 3 − k edges.
Proof. Denote the number of triangles in the triangulation as m. Every triangle has
three edges. k edges are the edge of only one triangle, the remaining edges are the
edges of two triangles. So the number of all edges is equal to
h =
3m + k
2
.
The Euler Formula says that
n − h + (m + 1) = 2.
Substituting h from the first formula we derive that
m = 2n − 2 − k.
From here substituting into the first formula we get
h = 3n − 3 − k.
1
2
Let T be a triangulation of the convex hull of a set P with m = 2n−2−k triangles.
The triangles of this triangulation have 3m angles. Let us order them into a sequence
α(T ) = (α1, α2, . . . , α3m), α1 ≤ α2 ≤ · · · ≤ α3m.
On these sequences we define the lexicographis order
α(T ) < α(T ),
if there is i ∈ {1, 2, . . . , 3m} such that
αj = αj for j < i and αi < αi.
A triangulation T is called angle-optimal if it is maximal in this ordering.
Legal edges and legal triangulations. For the purposes of our algorithm which
should create a suitable triangulation we will need the notion of legal edge a legal
triangulation. Before we start with a definition we repeat shortly high school geometry,
especially what we know about the angles on the arc of a circle.
Consider points p and q on the circle K with a center s. The points divide the
circle into two arcs. Choose a point a on the first of them and a point b on the
second. Denote the angles in the triangles paq and pbq at the vertices a and b α and
β, respectively. We call them arc angles of the chord pq. Arc angles have these two
importatnt properties:
• The size of angles α and β does not depend on the position of points a and b
on the given arcs.
• Their sum is α + β = 180o
.
From here we can easily derive that every triangle pcq where the point c lies in the
halfplane paq inside the circle K has the angle
γ > α
at the vertex c and that every triangle pdq where the point d lies in the halfplane paq
outside the circle K has the angle
δ < α
at the vertex d.
FIGURE 10.2 Arc angles
From the properties of arc angles it can be shown that we can curcimscribe a circle
to a quadrilateral if and only if the sum of opposite angles is equal to 180o
.
Now consider a triangulation T and in it an edge pq which has two adjacent triangles
pqr a pqt. Such an edge is called illegal edge of the triangulation T if the point t lies
inside the circle circumscribed to the triangle pqr. It is equivalent to the fact that the
sum of angles at vertices r a t in the quadrilateral prqt is greater than 180o
and also
to the fact that the point r lies inside the circle circumcsribed to the triangle pqt.
FIGURE 10.3 Illegal edge pq.
The edges of the triangulation T which are not illegal are called legal and the
triangulation without illegal edges is called legal triangulation. Immediate consequence
of this definition is the following characterization of legal triangulations:
3
Lemma 10.2. A triangulation is legal if and only if for every edge pq with adjacent
triagles pqr and pqt the point t does not lie inside the circle circumscribed to the
triangle pqr.
The following assertion will play a crucial role in our algorithm.
Lemma 10.3. Let T be a triangulation with illegal edge pq. Denote the adjacent
triangles pqr a pqt. From the triangulation T let us make a new triangulation T in
such a way that we replace the edge pq by the new edge rt. This edge is legal in the
triangulation T and moreover in the arrangement defined in the previous section
α(T ) < α(T ).
The exchange of the edge pq for the edge rt is called flip.
FIGURE 10.4 Flip
Proof. Since the edge pq is illegal in the triangulation T , the sum of the angles at
vertices r and t in the quadrilateral prqt is bigger than 180o
and so the sum of the
angles at p and q in the same quadrilateral is less than 180o
. (The sum of all angles
in the quadrilateral is 360o
.) That is why the edge rt in the triangulation T is legal.
Denote the angles before and after the flip according to the following figure:
FIGURE 10.5 Notation.
Since all the other triangles of both triangulations are the same, for
α(T ) < α(T ),
it suffices to prove that
min{αi, i = 1, 2, . . . , 6} < min{βi, i = 1, 2, . . . , 6}.
We show subsequently that every βi is bigger than an angle αj. Obviously,
β1 > α1 a β4 > α5.
Further, β2 is an arc angle over the chord pr in the circle circumscribed to the triangle
prt. Since the point q lies outside this circle, we get
β2 > α5.
Similarly β6 is an arc angle for the chord pt in the circle circumsribed to the triangle
prt. Hence
β6 > α4.
Analogously, using the circle circumsribed to the triangle qrt we can prove that
β3 > α1 a β5 > α2.
Corollary 10.4. An angle-optimal triangulation is legal.
Proof. If an angle-optimal triangulation had an illegal edge, using the flip of this edge
we would have got a triangulation which is bigger in the defined arrangement which
would have been a contradiction with the definition of the angle-optimal triangulation.
4
Legal triangulations need not be angle-optimal as the following figure shows on the
example of two triangulations of the pentagon p1p2p3p4p5 to which we can circumscribe
a circle. Moreover, p1p2p3p5 is a square and p1p2p4 is a isosceles triangle. Both
triangulations T and T are legal but because α(T ) < α(T ), the triangulation T is
not angle-optimal.
FIGURE 10.6 The example of legal triangulation which is not angle-optimal.
Delaunay graph and triangulation. Besides angle-optimal triangulation and legal
triangulation we introduce another notion of Delaunay triangulation. First we define
what Delaunay graph is.
The Delaunay graph of the set P has P as the set of vertices. Two vertices pi and
pj are connected by an edge if there is a circle on which these two points are lying and
all the other points of the set P lie outside the circle. It means, that the centre of this
circle lies on the edge of the Voronoi diagram which separates the cell determined by
pi from the cell determined by pj. That is why Delaunay graph is dual to the Voronoi
diagram.
FIGURE 10.7 The duality between the Voronoi diagram (dashed) and the Delaunay
graph (red).
Lemma 10.5. The Delaunay graph is a planar graph.
Proof. Suppose that pipj and pkpl are two edges of the Delaunay graph which intersect
in the inner point. Let the quadrilateral pipkpjpl has the sum of angles at the vertices
pk a pl at least 180o
. Consider a circle which has the chord pipj and pk lies outside it.
Then pl does not lie outside the circle. (Otherwise the sum of the angles at pk and pl
would be smaller than 180o
.) However this is a contradiction with pipj being an edge
of the Delaunay graph.
FIGURE 10.8 Illustration of the proof of Lemma 10.5
Bounded faces of the Delaunay graph are polygons. Let points p1, p2, . . . , pk form
the vertices of one of these polygons. This polygon determines a vertex v in the dual
graph, i.e. in the Voronoi diagram. Hence the points p1, p2, . . . , pk have the same
distance to the point v and lie on one circle with the center v.
FIGURE 10.9 The quadrilateral p1p2p3p4 in the Delaunay graph (black) and the
Voronoi diagram (red).
Delaunay trangulation is defined as an arbitrary triangulation of these polygons in
the Delaunay graph. The Delaunay triangulation need not be determined uniquely.
Nevertheless, the centres of the circles circumscribed to the triangles in the Delaunay
triangulation form the set of the vertices of the Voronoi diagram.
FIGURE 10.10 Delaunay graph (black) and Delaunay triangulation (black and red).
If no four points of the set P lie on one circle then the Delaunay graph is already
the Delaunay triangulation. In this case the Delaunay triangulation is determined
uniquely. We say that in this case the points of the set P are in general position.
5
Every Delaunay triangulation can be characterized by the following property.
Lemma 10.6. A triangulation of the convex hull of the set P is Delaunay if and only
if for every its edge pipj holds: the circle circumscribed to an adjacent triangle pipjpk
does not contain any other point of P in its interior.
Proof. To the proof we use the duality between the Delaunay graph and the Voronoi
diagram.
Let a triangulation is Delaunay. Then the triangle pipjpk has come into being by
dividing a polygon of the Delaunay graph into triangles. The centre of the circle
circumscribed to the triangle pipjpk is also the centre of this polygon and the vertex
in the Voronoi diagram. That is why in the interior of this circle cannot lie any other
point of the set P.
On the contrary, suppose that all edges of a triangulation have the property described
in the statement. Consider a triangle pipjpk and the circle circumscribed to
it. No remaining point of P lies in the interior of this circle. Some of them can lie on
the circle. So the triangle pipjpk is a part of polygon having the same circumscribed
circle. It is easy to see that to every edge of this polygon there is a circle having the
edge as a chord and having all the other points of P outside. So the edges of the
polygon are the edges of the Delaunay graph and our triangulation come into being
by the triangulation of the Delaunay graph, so it is Delaunay.
Using this lemma and the characterization of legal triangulations by Lemma 10.2,
we get immediately that every Delaunay triangulation is legal. To prove that every
legal triangulation is Delaunay is a little bit more complicated.
Theorem 10.7. The set of legal triangulations is equal to the set of Delaunay trian-
gulations.
Proof. Let us prove by contradiction that every legal triangulation is Delaunay. Let T
be a legal triangulation which is not Delauney. In this triangulation there is a triangle
pipjpk and a point pl such that pl lies in the interior of the circle circumsribed to pipjpk
and the angle ∠piplpj is maximal for all such triangles and points. Let pipjpm be an
adjacent triangle to the triangle pipjpk in the triangulation T . Since this triangulation
is legal the point pm does not lie in the interior of the circle K circumscribed to the
triangle pipjpk. The point pl cannot lie inside the triangle pipjpm. Let pl lies on the
opposite side of the line pjpm than the point pi. From the following figure we can see
that for the sizes of agles we have
∠pmplpj > ∠piplpj
which contradicts the choice of the triangle pipjpk and the point pl.
FIGURE 10.11 To the proof of Theorem 10.7
From here and from Corollary 10.4 we get:
Corollary 10.8. (1) Every angle-optimal triangulation is Delaunay.
(2) If the points of a set P are in general position, there is just one Delaunay triangulation,
just one legal triangulation and just one angle-optimal triangulation
and these triangulations are the same.
6
Basic idea of randomized incremental algorithm. One possibility how to construct
algorithmically Delaunay triangulation is to use algorithm for the Voronoi diagram,
to find the dual graph to it and divide its polygons into triangles.
We will present a differnt approach based on the fact that every legal triangulation
is Delaunay and which insists in removing gradually illegal edges.
A naive version of such an algorithm is following: Make a triangulation of the convex
hull of the set P and go through its edges. When you find an illegal edge, replace it
with a legal edge by a flip and proceed by the same way further. This process is finite
because after every flip we get a triangulation which is bigger than the previous one in
described lexicographic order. The number of triangulations is finite, so this process
must sometimes end.
A sofisticated version of this approach is the following randomized incremental algorithm.
Its basic idea can be captured in several points:
(1) We denote the point from the set P maximal in lexicographis order first with
respect y-coordinate and then x-coordinate as p0. We add points p−1 and p−2
to the points of the set P so that all other points of P lie inside the triangle
p0p−1p−2.
(2) We order these points randomly into a sequence p1, p2, . . . , pn. Gradually for
r = 1, 2, . . . , n we add the point pr to a legal triangulation Tr−1 of the triangle
p0p−1p−2 where the vertices of single triangles are in the set
p−2, p−1, p0, p1, . . . , pr−1.
(3) Using a simultaneously constructed search structure we find in which triangle
pipjpk of the triangulation Tr−1 the point pr lies and if it lies in the interior
or on the edge. Then we connect the point pr by edges with vertices of this
triangle (if it lies in the interior) or with vertices of adjacent triangles (if it lies
on an edge). This newly created triangulation of the triangle p0p−1p−2 with
vertices in the set
{p−2, p−1, p0, p1, . . . , pr−1, pr}
can be ”legalized” using flips which exchange illegal edges for legal ones. After
removing all illegal edges we get a triangulation Tr.
(4) At the end we remove the points p−1 and p−2 and the edges going from these
points from the legal triangulation Tn, Due to a suitable choice of the points p−1
and p−2 the remaining triangulation is the legal and also Delaunay triangulation
of the convex hull of the set P = {p0, p1, p2, . . . , pn}.
FIGURE 10.12 The set P in the triangle p0p−1p−2.
Now we describe the single steps in more details. We start with (2) and (3).
Legalization. Let us consider a legal triangulation Tr−1 described in the step (2)
above. If we find using a search structure that the point pr lies in the interior of a
triangle pipjpk, we connect it with single vertices by edges.
FIGURE 10.13 The point pr in the interior of a triangle pipjpk.
7
This will create a new triangulation. What is important is that the newly added
edges based on pr are legal. This can be proved as follows: Let K be a circle circumscribed
to the triangle pipjpk. Since the triangulation Tr−1 is legal and hence also
Delaunay, there is no point from the set Pr−1 = {p−2, p−1, p0, . . . , pr−1} in the interior
of K. Let L be a circle homothetic to K which has the cord prpi. See Figure 10.14.
In the interior of this circle there are no points from Pr−1. So the edge prpi is an edge
of the Delaunay graph for the set Pr and this has to be legal. Analogously, one can
prove that also the edges prpj and prpk are legal.
FIGURE 10.14 To the proof that prpi is legal.
It can happen that some of the edges pipj, pipk a pjpk become illegal in the new
triangulation. Using a flip we can replace them by edges going from the point pr. For
example let pipj be illegal. Then it is the edge not only of the triangle pipjpr but
elso of an adjacent triangle pipjpl. We replace it by the edge prpl which is already
legal. Now some of the edges pipl and pjpl can become illegal. That is why we have
to replace it by a flip with an edge based in the point pr.
FIGURE 10.15 The edge pipj is illegal.
We proceed in this way as long as there is no illegal edge. (After every flip the
new triangulation is bigger in the lexicographic order than the previous one.) The
procedure is captured in the following animation and pseudocode. The input is the
point pr, the edge pipj and a legal triangulation Tr−1.
ANIMATION
PSEUDOCODE Legalize(pr, pipj, Tr−1) from pseudo.pdf, page 37, but use the form
from the Czech version of e-learning.
If we find that the point pr lies on an edge pipj, then we replace the edge pipj by the
edges prpi and prpj and we create new edges connecting the point pr with the vertices
pk and pl of adjacent triangles to the edge pipj. Similarly as above one can show that
the new edges prpi, prpj, prpk a prpl are legal. However the edges pipk, pjpk, pipl a pjpl
can become illegal. This can be corrected by making gradually flips as in the above
case.
FIGURE 10.16 The point pr on the edge pipj.
Search structure. In the course of the algorithm we have a relevant search structure
for every triangulation. This is an oriented graph the leaves of which are the
triangles of the actual triangulation and inner nodes are the triangles of the previous
triangulations.
At the beginning the triangulation T0 has the only triangle p0p−1p−2 and the corresponding
search structure D0 has only one node given by this triangle. For the
triangulation Tr−1 we have a search structure Dr−1. We change the triangulation Tr−1
gradually in the way described in the previous section. Simultaneously we change
the search structure. To every transition in the triangulation there is a corresponding
transition from one search structure to further search structure. This is illustrated in
the Figure 10.17.
8
FIGURE 10.17 The change of the search structure corresponding to newly added
edges inside a triangle.
FIGURE 10.18 The change of the serch structure corresponding to a flip.
Using this procedure we create a search structure Dr corresponding to the legal
triangulation Tr. The resulting search structure depends only on the random order of
the points p1, p2, . . . , pr.
Initial choice of points p0, p−1, p−2. Now let us explain the steps (1) and (4) of our
algorithm. We will refer to the geometric idea for better understanding. Important is
that the points p−1 a p−2 will be determined by their properties, we will not determined
their concrete coordinates. The decisions whether some edge is legal or illegal, if among
four points, which are considered, is some of the points p−1, p−2, will not be made by
computations but only by applying the rules which follow from the properties of p−1
a p−2.
As we have already said for the point p0 we choose that one from P which is maximal
in the lexicographic order first with respect to the coordinate y and then with respect
to the coordinate x.
We choose the point p−1 to the right of the set P and under it. (More exactly: its
y-coordinate is smaller than the minimum of y-coordinates of the points from P and
its x-coordinate is bigger than the maximum of x-coordinates of the points from P.)
The choice of p−1 is carried out so that
• it lies outside all circles circumscribed to all triples of points from P which do
not lie in one line,
• all points pi with i > 0 lie under the line p0p−1.
We choose the point p−2 to the left of the set P and over it so that
• it lies outside all circles circumscribed to all triples of points from P ∪ {p−1}
which do not lie in one line,
• all points pi with i > 0 lie under the line p0p−2,
• all points pi with i ≥ 0 lie over the line p−1p−2.
So all the points pi with i > 0 lie in the interior of the triangle p0p−1p−2 and every
quadrilateral p−2pip−1pj is not convex which implies that the point p−1 lies outside all
circles circumscribed to the triangles p−2pipj.
FIGURE 10.19 The set P in the triangle p0p−1p−2.
Under these assumptions every Delaunay (legal) triangulation of P ∪ {p−1, p−2}
consists of
(1) the edges going from the point p−1 to the vertices of the right hand side of the
convex hull of the set P,
(2) the edges going from the point p−2 to the vertices of the left hand side of the
convex hull of the set P,
(3) a Delaunay triangulation of the set P.
Consequetly we can obtain a Delaunay triangulation of the set P in the way descrobed
in the step (4).
The rules how to work with our algorithm are based on the following lemma:
9
Lemma 10.9. If the points p0, p−1 a p−2 are chosen as described above the following
assertions are equivalent:
• pj lies to the left of the oriented line pip−1.
• pj lies to the left of the oriented line p−2pi.
• pj > pi in the lexicographic order first with respect to y and then with respect
to x.
Instead of doing a proof we demonstrate the situation by a picture.
FIGURE 10.20 To Lemma 10.9. The points p−1 and p−2 are chosen in such a way
that in the coloured angles there is no other point of the set P.
When finding weather an edge is legal or illegal we use the following rules:
• All edges of the triangle p0p−1p−2 are legal.
• The edge p−2pj with adjacent vertices p−1 and pk is legal.
FIGURE 10.21 Illustration of the previous statement.
• The edge p−1pj with adjacent vertices p−2 and pk is legal.
• Let the edge pipj have adjacent vertices pk and pl, let the quadrilateral pipkpjpl
be convex and let there be just one negative number among the numbers
i, j, k, l. Then the edge pipj is legal if and only if
min(k.l) < min(i, j).
FIGURE 10.22 min(−2, l) < min(i, j), the edge pipj is legal.
FIGURE 10.23 min(k, l) > min(−1, j), the edge p−1pj is illegal.
Pseudocode of the algorithm and its expected time complexity. A brief pseudocode
of the algorithm described in details in the previous sections can look like this:
PSEODOCODE DelaunayTriangulation(P) from pseudo.pdf, page 36
Without proof we give the following statement about expected time complexity of
the algorithm. See [Berg et al] for the proof.
Theorem 10.10. Expected time complexity of the randomized incremental algorithm
for the Delaunay triangulation of the set with n + 1 elements is O(n log n).