Voronoi Diagram and Delaunay Triangulation ========================================== Overview: --------- Define Voronoi Diagram and Delaunay Triangulation Explain Duality Motivate their usefulness Examine numerous properties of both For example: Connection between minimum spanning trees and Delaunay triangulation Devise and analyze algorithms for constructing them Connection between convex hulls and Delaunay triangulation Definition: ----------- Voronoi Diagram: Believed to be first formally introduced by the Russian mathematician Georges Voronoi in 1907. Was discussed 50 years earlier by Dirichlet; sometimes the Voronoi digram is called Dirichlet tessellation. There are notes indicating that Rene Descartes has used Voronoi diagrams as early as 1644. Given n sites in the plane, the Voronoi diagram divides the plane into regions associated with each site, such that all points in a region are closest to the point associated with that region. A Voronoi diagram consists of cells associated with a single site, edges, equidistant to two sites, and vertices, equidistant to three sites. Since we assume general position, no four sites can be equidistant to a point in the plane, i.e., no four sites are co-circular. More formally, the Voronoi region of a site s_i is defined: V(s_i) = {x : |x - s_i| <= |x - s_j| for all j =/= i} Delaunay Triangulation: Boris Nikolaevich Delaunay (1890-1980) also was a Russian mathematician and worked on mathematical crystallography. In 1934 he proved that dual graph of the Voronoi diagram produces a planar triangulation of the Voronoi sites. Best defined as the dual of the Voronoi diagram (see below). It has many very interesting properties: * maximizes minimum angle in each triangle * minimizes maximum radius of circumcircle and inscribed circle * minimizes sum of inscribed radii * many more See Voronoi/Delaunay applet: http://www.cs.cornell.edu/Info/People/chew/Delaunay.html Duality: -------- Two Voronoi regions associated with Voronoi sites a,b share an edge e <===> The sites a,b are connected by a straight edge in the Delaunay triangulation. Applications: ------------- Just a few applications of the Voronoi diagram: Anthropology and Archaeology -- Identify the parts of a region under the influence of different neolithic clans, chiefdoms, ceremonial centers, or hill forts. Astronomy -- Identify clusters of stars and clusters of galaxies (Here we saw what may be the earliest picture of a Voronoi diagram, drawn by Descartes in 1644, where the regions described the regions of gravitational influence of the sun and other stars.) Biology, Ecology, Forestry -- Model and analyze plant competition ("Area potentially available to a tree", "Plant polygons") Cartography -- Piece together satellite photographs into large "mosaic" maps Crystallography and Chemistry -- Study chemical properties of metallic sodium ("Wigner-Seitz regions"); Modeling alloy structures as sphere packings ("Domain of an atom") Finite Element Analysis -- Generating finite element meshes which avoid small angles Geography -- Analyzing patterns of urban settlements Geology -- Estimation of ore reserves in a deposit using information obtained from bore holes; modeling crack patterns in basalt due to contraction on cooling Geometric Modeling -- Finding "good" triangulations of 3D surfaces Marketing -- Model market of US metropolitan areas; market area extending down to individual retail stores Mathematics -- Study of positive definite quadratic forms ("Dirichlet tessellation", "Voronoi diagram") Metallurgy -- Modeling "grain growth" in metal films Meteorology -- Estimate regional rainfall averages, given data at discrete rain gauges ("Thiessen polygons") Pattern Recognition -- Find simple descriptors for shapes that extract 1D characterizations from 2D shapes ("Medial axis" or "skeleton" of a contour) Physiology -- Analysis of capillary distribution in cross-sections of muscle tissue to compute oxygen transport ("Capillary domains") Robotics -- Path planning in the presence of obstacles Statistics and Data Analysis -- Analyze statistical clustering ("Natural neighbors" interpolation) Zoology -- Model and analyze the territories of animals Applications for Delaunay: Finite element analysis Function interpolation Properties of the Voronoi diagram and the Delaunay Triangulation: ----------------------------------------------------------------- * V01: All Voronoi regions are convex polygons. The Voronoi region is the intersection of halfplanes created by all the other sites, hence it is a convex polygon. * * V02: A site s has an unbounded regions if and only if it s lies on the convex hull of the sites. <= If s is a site on the convex hull, then we take a supporting line l through x such that all the sites lie on the same side of l. Every point on a ray normal to l and through s away from the other points is closest to s, thus part of its Voronoi regions. => If s has an unbounded region, since the region is convex, there must be an infinite ray from s. A line l normal to this ray through s, we see that all other sites must lie on one side of l, therefore s is an extreme point. * * V03: For n>=3 sites, the Voronoi diagram has at most 2n-5 vertices and 3n - 6 edges. v <= 2n-5 e <= 3n-6 Euler's formula states that every planar connected graph fulfills the following: v - e + f = 2 We introduce a vertex at infinity and connect all half-lines from open Voronoi regions to that point to obtain a connected graph. We get: v+1 - e + f = 2 Every vertex has at least degree three, we can bound: 2e >= 3(v+1) So: (v+1) - e + n = 2 2(v+1) - 2e + 2n = 4 2(v+1) - 3(v+1) + 2n >= 4 2n - (v+1) >= 4 2n -5 >= v (v+1) - e + n = 2 3(v+1) -3e + 3n = 6 2e - 3e + 3n >= 6 3n - e >= 6 3n - 6 >= e * The complexity of the Voronoi diagram in 2D is linear; in nD it is exponential. * V04: A point q is a vertex of the Voronoi diagram if and only if its largest empty circle C(q) contains exactly three sites on its boundary (again assuming general position of the sites). <= If a q and C(q) with three sites on its boundary exists then q is clearly equidistant to all three sites. Then it must be a vertex of the Voronoi diagram. => A vertex q of the Voronoi diagram must be equidistant to at least three sites, corresponding to the three Voronoi regions adjacent to q. If there were any other sites closer than the three considered, the Voronoi regions could not be adjacent to q. * * V05: The bisector of sites s_i and s_j defines a Voronoi edge if and only if there exists a point q on the bisector such that C(q) only contains s_i and s_j on its boundary. <= Suppose there is a point q on the bisector such that C(q) only contains sites s_i and s_j that give rise to the bisector. Since all other sites are further away, there has to be an edge between the Voronoi regions of s_i and s_j that q lies on. => Let q be a point on the Voronoi edge of s_i and s_j. By definition, s_i and s_j must be equidistant to q and no other site can be close. Therefore, the maximum circle around q only contains s_i and s_j on its boundary. * * D01: An edge ab is part of the Delaunay triangulation if and only if there is a circle (called witness to the Delaunay-hood of the edge) through a and b that contains no other site. => Let ab be an edge in the Delaunay triangulation. Sites a and b - by definition - must share a Voronoi edge. Take a point p in the interior of this edge, then a and b are at the same distance r from p, and all other sites are further away. Hence, the circle of radius r about r is the witness. <= The existence of a circle through a and b with all other sites strictly outside implies that the center of the circle is on a boundary edge between the Voronoi regions of a and b. * * D02: The Delaunay triangulation is planar. Assume that two edges ab and xy in the Delaunay triangulation cross. Then there have to be two circles with a,b and x,y on their boundary, respectively. The two corresponding circles would have to intersect in four points. Contradiction. * * D03: e = 3(n-1) -k f=2(n-1) -k (k = size of the convex hull) Again, assuming f triangles and using Euler's formula, after regarding the outer face as a "triangle", we get: n - e + (f+1) = 2 n - 3 + f = 1 Since every edge belongs to two faces we have 3f+k = 2e. We get: e = 3(n-1) -k f = 2(n-1) -k * By duality we obtain an exact bound for the Voronoi diagram, improving V03. * D04: No site lies inside the circumcircle of a triangle of the Delaunay triangulation. Consider the center of such a circle. It must be a vertex of the Voronoi diagram, because it is equidistant to three sites. According to V4 the associated circle must be empty. * * D05: We can check the Delaunay-hood of a triangulation in O(n^2). Sure, using D04, we can just go through all triangles and for each triangle check if each site is outside its circumcircle. Since we have O(n) triangles and O(n) sites, this yields an O(n^2) algorithm. * * D06: (Theorem) Given a triangulation of n sites such that for every pair of adjacent triangles abc and bcd, a is not in the circumcircle of bcd, then that triangulation is the Delaunay triangulation. See slides for D06. Note that this is stronger then D04. Now we claim that we only have to check the sites of adjacent triangles, rather than all other sites. Delaunay proved this astonishing theorem. We define the power of a point x wrt/ a circle C as follows. Let l be any line through x intersecting C in points a and b. The power of x wrt/C is ax * bx. This is independent of what line you chose. Signs are chosen such that power(x,C) > 0 (=0, <0) iff x is outside (on, inside) C. Given a triangulation, we pick a triangle abc and a site x across from bc. If xbc is part of the triangulation, we are done, because this is the case in which we will check anyway. Otherwise, we have a point x that is outside of the triangle abc. Pick a point d inside abc so that xd crosses bc (that choice is arbitrary). As we move along xd we must pass through triangles. Let's say the first one is vwx, containing x as a vertex. We obviously end in abc. Let's take a look at the circumcircles of those triangles. The first one goes through vwx and the power(x,C(vwx))=0, because x is on vwx. We will argue that the power of x with respect to any subsequent triangle's circumcircle can only increase, as we move away from x. Therefore, power(x,C(abc))>0 and hence x is outside of the circumcircle of abc, we're done. Now, is it true that the power has to increase? Consider two triangles lmn and klm, and say that moving from x to d we will move through lmn first. The circumcircles of both triangles are members of the family of circles through l and m. Imagine pushing C(lmn) through the points l and m to obtain C(klm). Both intersection points of a line l through x and the circle will move further away and hence the power will increase. * * D07: Now we can check Delaunay-hood in O(n). The decision problem for Delaunay is O(n). We use D06 to argue that we only have to check adjacent triangles. * * D08: (Theorem) The minimum spanning tree of a set of points (this holds in any dimensions) is a subgraph of the Delaunay triangulation. See slides for D08. Minimum spanning tree can be computed with Kruskal's algorithm in O(m log m), where m is the number of edges in the graph. Applying it to n sites in the plane, we'd have n^2=m edges and would get a running time of O(n^2 log n). We can compute the Delaunay triangulation in O(n log n) and then use the O(n) edges of this triangulation for Kruskal's algorithm to get the result in O(n log n), factor of n faster than by just using Kurskal's algorithm. But we first have to show that the theorem is true. So let's prove the theorem D08. Assume T is a minimum spanning tree and assume that edge ab is part of that tree but not an edge of the Delaunay triangulation. Therefore, no empty circle through ab exists. Let's pick the on with diameter ab. There must be a cite c inside this circle. Let's remove ab. It is easy to see that adding either ac or ab would reconnect the spanning tree with a smaller weight that it had with ab. Contradiction. * Fortune's Sweepline Method to Compute the Voronoi Diagram --------------------------------------------------------- Naive method of computing each individual Voronoi region by intersecting halfplanes can be done in O(n^2 log n). (Intersection of n halfplanes is dual to computation of convex hull.) Incremental algorithm by Green & Sibson 1977 uses O(n) per insertion for a total of O(n^2). Shamos & Hoey 1975 presented a divide-and-conquer algorithm, which is hard to implement, but runs in O(n log n) time. In 1985 Fortune introduced a sweepline algorithm that compute the Voronoi diagram in O(n log n). But how can a sweepline algorithm work? The shape of Voronoi region to be constructed by the sweep can depend on vertices that are beyond the sweepline! Fortune solved this with a very clever trick. Cones, representing time. Sweepline is a sweep plane at 45 degree. View from z=-infinity. Intersects cones in parabolas, called parabolic front. Where parabola meet, must be Voronoi edge. [lots of stuff missing here] Fortune avoided the problem of the sweepline by not constructing the Voronoi diagram to the left of the line, but underneath the plane, so it lags behind the line. Sweeping the Cones (see slide 21): Lemma 1: The two triangles are identical: They share an edge and have identical angles. See slide to convince yourself. Parabolic Front (see slide 22): The intersection of the sweep plane with the cones. The parabolic front is intersected by any horizontal line exactly once. It is a single chain of arcs from -infinity to + infinity. To the left of it, we see the Voronoi diagram. Evolution of the parabolic front (see slide 23): We need to show that the breakpoints of the parabolic front trace *all* the points on the Voronoi edges. Convince yourself also that the other direction (that no other points can be breakpoints than those on Voronoi edges) is obvious by construction. Lemma 2: Every point of every edge of the Voronoi digram is a breakpoint of the parabolic front at some time during the sweep. Let e be a Voronoi edge and u be an internal point on e. We know that pu = qu = distance from u to the sweepline, let's call it r. We know that no other site can be in the circle C with radius r around u (V05, see above). Therefore, u is on the intersections of the parabola arising from the sweep plane cutting the cones of p and q. But u could be cut off by some other cone, i.e., parabola, between u and the sweepline (u would be "behind" this parabola). Assume a site v is on such a parabola, just to the right of u. The site giving rise to the cone that gives rise to the parabola v is on, must be on the circle with center v that touches the sweepline. But this circle is inside C and we know there cannot be another site. * Okay, now we know that if we keep track of the parabolic front, we'll trace the Voronoi diagram. But how to do that? As always with sweepline techniques, we have to identify discrete events that we can process in order of their occurrence. There are two types of events happening on the parabolic front: A new parabola is introduced, because of a new site that is being swept (we'll call this a site event), or a parabola disappears because the two parabola around it "eat" it (we'll call this a circle event for reasons that will become clear in Lemma 4). Lemma 3: (Site Event, see slide 24) The only way for an arc to enter the parabolic front is through a site event. What could go wrong? How could another parabolic segment of the parabolic front be introduced if not by another site? Some parabola alpha could "overtake" another one (beta) and break through it, splitting it up into two pieces and inserting itself in between. However, if alpha is just about to overtake beta, then it would have to have higher curvature (after all it is supposed to "pierce" beta) but at the same time be farther away from the sweepline (it has not broken through yet!). This is contradictory, and therefore no alpha can overtake beta by cutting in between. But is it possible for alpha to overtake the intersection of two parabolas beta and gamma? To do that, it would have to move faster than the parabolas beta and gamma at their intersection, let's call the point of intersection u. But it turns out that the "sides" of the parabola (where the absolute slope is smaller) advance faster than the middle (you can do some algebra to prove it or look at the animation at the link provided below to convince yourself; this must be the case because ultimately the parabola degenerates into a line, that is its sides "catch up"). Lemma 4: (Circle Event, see slide 25) The only way for an arc to leave the parabolic front is through a circle event. A circle event occurs when the sweep line becomes tangent to an empty circle around a Voronoi vertex with three sites on its boundary. Given the property V04, and looking at the pictures on the slide, this is obvious. We now need to discuss a data structure to maintain the parabolic front and a priority queue of future events. These data structures will support insertion in O(log n) so that the n events can be handled in O(n log n). For more details see your textbook or section A4 of a paper by Guibas and Stolfi, called "Ruler, Compass, and Computer: The Design and Analysis of Geometric Algorithms," posted to the class web site. Data structures: Event queue contains all site events after sorting in O(n log n), and all the circle events for three consecutive arcs on the parabolic front. Note that not necessarily all consecutive arcs define a circle event. As new arc adjacencies imply new circle events they are inserted into the queue. The parabolic front is stored as a balanced tree. Leaves correspond to parabolic arcs, nodes to break points. The leaves contain pointers to the corresponding site, and the break points to the Voronoi edge they create. Event scheduling (see slide 26): All site events are in the queue in sorted order. When a new arc beta is added, the corresponding circle events (two, one with beta as the first, one with beta as the last element). When new adjacencies arise due to a deleted arc, we need to check if the circle event occurs to the left of the sweep line or to the right, in which case it was already processed and does not need to be inserted into the event queue. There is a detail about inserted events being invalidated by the introduction of new events (see the two left figures on the slide). We won't go into more detail here, but this is the result: Theorem: Fortune's algorithm computes the Voronoi diagram of n sites in O(n log n) time and O(n) space. See animation of Fortune's algorithm at: http://www.diku.dk/hjemmesider/studerende/duff/Fortune/ Randomized Incremental Method for Constructing the Delaunay triangulation ------------------------------------------------------------------------ We first describe a very natural incremental method that is not randomized. Its runtime will be worst case O(n^2). The we show how by introducing randomization -- a new tool that we have so far not used in this class -- we can modify the algorithm to obtain an expected running time of O(n log n). Of course, the difference in running time between the two methods will be mainly in the analysis. But it is possible to construct sequences of sites that when inserted in order will cause the algorithm to be O(n^2). The advantage of the randomized analysis is that we can consider the expected running time. Incremental Method: ------------------- Imagine "huge" triangle that contains all sites s_1, ... , s_n that we want to insert. For every site s to be inserted: - locate the triangle xyz with circumcircle C that s falls into (we can do this in O(log n) and will describe it later) - add edges sx sy sz these edges are Delaunay: the circle at through s, tangent to C at x is a Delaunay witness of sx - edges xy yz zx are suspect and need to be checked if they fail the Delaunay test the edge is swapped and two new suspect edges are created - spread "wave of suspicion" until no more edges are suspicious or the convex hull of the points is reached Lemma: When a new site s is inserted, no edge is tested more than once. Proof: The wave of suspicion moves away from the new site. Edges on the convex hull are always Delaunay. Running time: O(n) sites are added and for each insertion O(n) edges need to be checked, therefore the running time is O(n^2). Randomized Method: ------------------ Same as above, except that we add the sites in random order. Theorem: The number of Delaunay triangles arising during the execution of the algorithm is O(n). Proof: Provided in class. Rest of this after randomized point location. Randomized Point Location: -------------------------- See handout. Connection between Delaunay triangulation and convex hulls ---------------------------------------------------------- See handout.