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).