Lecture 17:  Voronoi Diagrams and Delauney Triangulations.

In the last lecture, we covered Fortune's algorithm for creating a
voronoi diagram.  The following applet given a demonstration of this
algorithm at work:

voronoi sweep line

The dual of the Voronoi Diagram is called the Delaunay
Triangulation.  In today's class we will discuss a number of
properties of the Delaunay traingulation.  To give you an intuitive
idea of what the relationship is, you can look at the following
applet:

Voronoi and Delaunay

So let us consider more formally the Delaunay triangulation.

Delaunay Triangulations: Last time we gave an algorithm for computing
Voronoi diagrams. Today we consider the related structure, called a
Delaunay triangulation (DT). Since the Voronoi diagram is a planar
graph, we may naturally ask what is the corresponding dual graph. The
vertices for this dual graph can be taken to be the sites
themselves. Since (assuming general position) the vertices of the
Voronoi diagram are of degree three, it follows that the faces of the
dual graph (excluding the exterior face) will be triangles. The
resulting dual graph is a triangulation of the sites, the Delaunay
triangulation.


Figure 60: Delaunay triangulation.

Delaunay triangulations have a number of interesting properties, that
are consequences of the structure of the Voronoi diagram.

Convex hull: The exterior face of the Delaunay triangulation is the
convex hull of the point set.

Circumcircle property: The circumcircle of any triangle in the
Delaunay triangulation is empty (contains no sites of P ).

Empty circle property: Two sites pi and pj are
connected by an edge in the Delaunay triangulation, if and only if
there is an empty circle passing through pi and
pj . (One direction of the proof is trivial from the
circumcircle property. In general, if there is an empty circumcircle
passing through pi and pj , then the center c of
this circle is a point on the edge of the Voronoi diagram between
pi and pj , because c is equidistant from each
of these sites and there is no closer site.)

Closest pair property: The closest pair of sites in P are neighbors in
the Delaunay triangulation. (The circle having these two sites as its
diameter cannot contain any other sites, and so is an empty circle.)

If the sites are not in general position, in the sense that four or
more are cocircular, then the Delaunay triangulation may not be a
triangulation at all, but just a planar graph (since the Voronoi
vertex that is incident to four or more Voronoi cells will induce a
face whose degree is equal to the number of such cells). In this case
the more appropriate term would be Delaunay graph. However, it is
common to either assume the sites are in general position (or to
enforce it through some sort of symbolic perturbation) or else to
simply triangulate the faces of degree four or more in any arbitrary
way. Henceforth we will assume that sites are in general position, so
we do not have to deal with these messy situations.


Given a point set P with n sites where there are h sites on the convex
hull, it is not hard to prove by Euler's formula that the Delaunay
triangulation has 2n-2-h triangles, and 3n - 3 - h edges. (The ability
to predict the number of triangles from n and h only works in the
plane. In 3 space, the number of triangles in the Delaunay triangles
can range from O(n) up to O(n2 ). In dimension n, the number of
triangles can range as high as O(nfloor(d/2) ). This is not easy to see
and requires some pretty deep understanding of the combinatorics of
convex polytopes.)

Minimum Spanning Tree: The Delaunay triangulation possesses some
interesting properties that are not directly related to the Voronoi
diagram structure. One of these is its relation to the minimum
spanning tree. Given a set of n points in the plane, we can think of
the points as defining a Euclidean graph whose edges are all (n choose
2) (undirected) pairs of distinct points, and edge (pi,
pj) has weight equal to the Euclidean distance from
pi to pj. A minimum spanning tree is a set of n
- 1 edges that connect the points (into a free tree) such that the
total weight of edges is minimized. We could compute the MST using
Kruskal's algorithm. Recall that Kruskal's algorithm works by first
sorting the edges and inserting them one by one. We could first
compute the Euclidean graph, and then pass the result on to Kruskal's
algorithm, for a total running time of O(n2 log n).

However there is a much faster method based on Delaunay
triangulations. First compute the Delaunay triangulation of the point
set. We will see later that it can be done in O(n log n) time. Then
compute the MST of the Delaunay triangulation by Kruskal's algorithm
and return the result. This leads to a total running time of O(n log
n). The reason that this works is given in the following theorem.

Theorem: The minimum spanning tree of a set of points P (in any
dimension) is a subgraph of the Delaunay triangulation.

Proof: Let T be the MST for P , let w(T ) denote the total weight of T. 
Let a and b be any two sites such that ab is an edge of T . Suppose
to the contrary that ab is not an edge in the Delaunay
triangulation. This implies that there is no empty circle passing
through a and b, and in particular, the circle whose diameter is the
segment ab contains a site, call it c. (See the figure below.)

The removal of ab from the MST splits the tree into two
subtrees. Assume without loss of generality that c lies in the same
subtree as a. Now, remove the edge ab from the MST and add the edge bc
in its place. The result will be a spanning tree T' whose weight is

w(T') = w(T) + |bc| - |ab| > w(T ):

The last inequality follows because ab is the diameter of the circle,
implying that |bc| < |ab|.  This contradicts the hypothesis that T is
the MST, completing the proof.


Figure 61: The Delaunay triangulation and MST.

By the way, this suggests another interesting question. Among all
triangulations, we might ask, does the Delaunay triangulation minimize
the total edge length? The answer is no (and there is a simple four
point counterexample.  The parenthetical comments typically indicate
something that you should work out.  Got it?!).  However, this claim
was made in a famous paper on Delaunay triangulations, and you may
still here it quoted from time to time. The triangulation that
minimizes total edge weight is called the minimum weight
triangulation. To date, no polynomial time algorithm is known for
computing it, and the problem is not known to be NP complete.

3D Convex Hulls and Delaunay Triangulations There is a
fascinating relationship between delaunay triangulations and the
convex hull of a particular set of 3D points.


At first, Delaunay triangulations and convex hulls appear to be quite
different structures, one is based on metric properties (distances)
and the other on affine properties (collinearity, coplanarity). Today
we show that it is possible to con vert the problem of computing a
Delaunay triangulation in dimension d to that of computing a convex
hull in dimension d + 1. Thus, there is a remarkable relationship
between these two structures.

We will demonstrate the connection in dimension 2 (by computing a
convex hull in dimension 3). Some of this may be hard to visualize.
The connection between the two structures is the paraboloid z =
x2 + y2. Observe that this equation defines a
surface whose vertical cross sections (constant x or constant y) are
parabolas, and whose horizontal cross sections (constant z) are
circles. For each point in the plane, (x, y), the vertical projection
of this point onto this paraboloid is (x,y, x2
+y2) in 3 space.

Given a set of points S in the plane, let S0 denote the
projection of every point in S onto this paraboloid.  Consider the
lower convex hull of S0. This is the portion of the convex
hull of S0 which is visible to a viewer standing at z =
-1. We claim that if we take the lower convex hull of S0,
and project it back onto the plane, then we get the Delaunay
triangulation of S. In particular, let (p, q, r) be elements of S, and
let p0 , q0 , r0 denote the
projections of these points onto the paraboloid.  Then p0
q0 r0 define a face of the lower convex hull of
S0 if and only if pqr is a triangle of the Delaunay
triangulation of S. The process is illustrated in the following
figure.


Figure 93: Delaunay triangulations and convex hull.

The question is, why does this work? To see why, we need to establish
the connection between the triangles of the Delaunay triangulation and
the faces of the convex hull of transformed points. In particular,
recall that

Delaunay condition: Three points p, q, r, in S form a Delaunay
triangle if and only if the circumcircle of these points contains no
other point of S.

Convex hull condition: Three points p0,
q0, r0 in S0 form a face of the
convex hull of S0 if and only if the plane passing through
p0, q0, and r0 has all the points of
S0 lying to one side.

Clearly, the connection we need to establish is between the emptiness
of circumcircles in the plane and the emptiness of halfspaces in 3
space. We will prove the following claim.

Lemma: Consider 4 distinct points p, q, r, s in the plane, and
let p0, q0, r0 s0 be their
respective projections onto the paraboloid, z = x2 +
y2 . The point s lies within the circumcircle of p, q, r if
and only if s0 lies on the lower side of the plane passing
through p0, q0, r0.

To prove the lemma, first consider an arbitrary (nonvertical) plane in
3 space, which we assume is tangent to the paraboloid above some point
(a, b) in the plane. To determine the equation of this tangent plane,
we take derivatives of the equation z = x2 + y2
with respect to x and y, giving:

dz/dx = 2x

dz/dy = 2y

At the point (a, b, a2 +b2 ) these evaluate to
2a and 2b. It follows that the plane passing through these point has
the form

z = 2ax + 2by + k:

To solve for k we know that the plane passes through 
(a, b, a2 + b2 ) 
so we solve giving 

a2+ b2 = 2 a2 + 2 b2 + k;

Implying that k  = - (a2 + b2  ). Thus the plane
equation is z = 2ax + 2by - (a2 + b2 ):

If we shift the plane upwards by some positive amount r2 we
get the plane

z = 2ax + 2by - (a2 + b2 ) + r2 :

How does this plane intersect the paraboloid? Since the paraboloid is
defined by z = x2 + y2

we can eliminate z giving

x2 + y2 = 2ax + 2by - (a2 + b2 ) + r2.

which after some simple rearrangements is equal to

(x - a)2 + (y - b)2 = r2

This is just a circle. Thus, we have shown that the intersection of a
plane with the paraboloid produces a space curve (which turns out to
be an ellipse), which when projected back onto the (x, y) coordinate
plane is a circle centered at (a, b).

Thus, we conclude that the intersection of an arbitrary lower
halfspace with the paraboloid, when projected onto the (x, y) plane is
the interior of a circle. Going back to the lemma, when we project the
points p; q; r onto the paraboloid, the projected points p0
, q0 and r0 define a plane. Since p0,
q0, and r0, lie at the intersection of the plane
and paraboloid, the original points p, q, r lie on the projected
circle. Thus this circle is the (unique) circumcircle passing through
these p, q, and r. Thus, the point s lies within this circumcircle, if
and only if its projection s0 onto the paraboloid lies
within the lower halfspace of the plane passing through


Figure 94: Planes and circles.

An applet that show simultaneously the delaunay triangulation and the
3D convex hull can be found 
here
It may be necessary to try a
collection of patterns or sets of points (click on the left coordinate
system near the bottom of the page).