Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00003998/00001
## Material Information- Title:
- Voronoi diagrams properties, algorithms, and applications
- Creator:
- Romano, Don
- Publication Date:
- 1997
- Language:
- English
- Physical Description:
- v, 87 leaves : illustrations ; 29 cm
## Subjects- Subjects / Keywords:
- Voronoi polygons ( lcsh )
Voronoi polygons ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaves 81-87).
- General Note:
- Submitted in partial fulfillment of the requirements for the degree, Master of Science, Computer Science.
- General Note:
- Department of Computer Science and Engineering
- Statement of Responsibility:
- by Don Romano.
## Record Information- Source Institution:
- |University of Colorado Denver
- Holding Location:
- |Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 37832088 ( OCLC )
ocm37832088 - Classification:
- LD1190.E52 1997m .R66 ( lcc )
## Auraria Membership |

Full Text |

VORONOI DIAGRAMS:
PROPERTIES, ALGORITHMS, AND APPLICATIONS by Don Romano B.S., University of Colorado, 1993 A thesis submitted to the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Computer Science 1997 This thesis for the Master of Science degree by Don Romano has been approved by aLl2j227 Date Romano, Don (M.S., Computer Science) Voronoi Diagrams: Properties, Algorithms, and Applications Thesis directed by Assistant Professor Jay Rothman ABSTRACT This thesis presents an introduction to the concept of the Voronoi diagram, a fundamental structure in computational geometry. A number of mathematical properties of the Voronoi diagram and related geometric structures such as the Delaunay tessellation are examined. Particular emphasis is placed upon several efficient methods for the construction of the Voronoi diagram. Computer implementation details, algorithm performance, and a data structure used for the representation of Voronoi diagrams are also described. A survey of applications of the Voronoi diagram is provided to illustrate the usefulness of this important geometric structure. This abstract accurately represents the content of the candidates thesis I recommend its publication. Signed Jay Rothman m CONTENTS Chapter 1. Introduction.............................................................1 1.1 The Nearest Neighbor Problem..........................................1 1.2 The Voronoi Diagram...................................................2 1.3 The Delaunay Tessellation.............................................3 1.4 Formal Definitions....................................................4 2. Properties of Voronoi Diagrams and Delaunay Tessellations...............8 2.1 Properties of the Voronoi Diagram.....................................8 2.2 Properties of the Delaunay Tessellation..............................14 2.3 Graph-theoretic Properties of 'V(P) and T)(P)........................17 2.3.1 Abstract and Geometric Graphs........................................18 2.3.2 Eulers Formula......................................................19 2.3.3 Voronoi Graphs.......................................................20 2.3.4 Delaunay Graphs......................................................23 2.3.5 Subgraphs of the Delaunay Triangulation..............................23 3. Representing Voronoi Diagrams in a Data Structure.......................27 3.1 Winged-edge Data Structure...........................................27 IV 3.2 Geometric Information Retrieval 32 3.3 Complexity of Constructing 'V(P)................................... 33 3.4 Other Data Structures Used To Represent 'V(P).......................35 4. Algorithms For Constructing Voronoi Diagrams...........................36 4.1 The Polygon Construction Algorithm..................................38 4.2 The Incremental Method..............................................40 4.3 The Divide-and-Conquer Method.......................................47 4.4 The Plane Sweep Method..............................................56 4.5 Comparative Performance of Construction Methods.....................62 4.6 Numerical Accuracy and Implementation Issues........................63 4.7 Constructing the Voronoi Diagram in Higher Dimensions...............65 4.8 Algorithms for Constructing the Delaunay Tessellation...............67 5. Applications of Voronoi Diagrams.......................................69 5.1 Applications in Mathematics.........................................69 5.2 Applications in the Natural Sciences................................75 5.3 Applications in Computer Science....................................77 References.................................................................81 v 1. Introduction A number of problems arise in different fields of study which require solutions involving proximity relationships between distinct points and objects in n-dimensional space. More specifically, the frequent need to solve what is known as the nearest neighbor problem, has motivated considerable research into Voronoi diagrams. Being very well-suited to provide solutions to these problems, the Voronoi diagram along with its generalizations and related structures have become a focus of computational geometry. This thesis presents a survey of this rich and fascinating subject. To lay a groundwork for the sections that follow, this chapter will describe the basic nature of the Voronoi diagram and its geometric dual, the Delaunay tessellation. After an introductory description of these diagrams, the formal mathematical definitions are presented. Voronoi diagrams and their duals possess a number of unique and useful propertiescharacteristics which are often exploited by specific applications. Chapter 2 presents a number of important properties associated with these remarkable structures. Several properties of the graph- theoretic representations are also given. Voronoi diagrams hold much information about the space in which they exist. This information must be efficiently stored and easily retrieved for the diagrams to be useful in practical applications. In Chapter 3, we define a common data structure for representing a Voronoi diagram. Several algorithms for the construction of Voronoi diagrams, along with analysis of their performance, are given in Chapter 4. Finally, and for a better understanding of the usefulness of these structures, Chapter 5 presents a survey of the principal applications. It is assumed that the reader is familiar with the fundamental concepts of plane geometry, graph theory, and data structures. Common set notation and the nomenclature used to represent algorithm performance (i.e., big-0 notation) are also used without preliminary explanation. 1.1 The Nearest Neighbor Problem Nearest neighbor problems occur quite often and are truly interdisciplinary. Problems as diverse as locating the placement of public facilities in a city, finding the regions of influence for certain species of animals, studying the growth of 1 crystals, and determining the distribution of galaxies in the universe are but a few of the myriad instances of this wide-ranging class of applications. During the last two decades, the evolution of high-speed computers and the development of efficient methods of constructing Voronoi diagrams have broadened the use of nearest neighbor analysis. In Chapter 5, we provide an overview of nearest neighbor and related problems from computer science, mathematics, and the natural sciences. Other disciplines including, but not limited to, physics, meteorology, regional planning, geography, cartography, and operations research have all used nearest neighbor analysis. 1.2 The Voronoi Diagram In general, given a finite set of distinct, fixed points in the real plane (or any continuous space), the solution to the nearest neighbor problem answers the question which site is closest to a chosen point p in the space?. This problem is solved by computing the Voronoi diagram for the original set of points. The resulting geometric structure partitions the space into a set of regions, each of which is associated with one of the Voronoi points. An intuitive way of looking at this is to see that for each point in the Voronoi set, there exists a region of the space closest to it. All locations within a particular region are closer to the associated Voronoi point defining the region than to any other Voronoi point in the space. Figure 1.1 illustrates a simple planar example of a Voronoi diagram for 10 points. Figure 1.1 A Voronoi diagram for 10 points in the plane. 2 1.3 The Delaunay Tessellation The structure most closely related to the Voronoi diagram is the Delaunay tessellation. Every Voronoi diagram has a Delaunay tessellation for its geometric (and graph-theoretic) dual. These tessellations are simply triangulations in the plane, tetrahedral tessellations in R3, and more generally, a tessellation of simplices in Rn. For purposes of clarity, we will confine our discussion (in this introduction) to diagrams and tessellations in the plane, while keeping in mind the immediate extensions to higher-dimensional spaces. Delaunay tessellations can be easily constructed from the corresponding Voronoi diagrams by joining those points whose regions share a line (or face, or simplex). Figure 1.2 shows the Delaunay triangulation which is the dual of the Voronoi diagram illustrated in Figure 1.1 above. It is interesting to note that a Delaunay tessellation contains the same information as its Voronoi dual, but the information is represented differently. The vertices of the Voronoi diagram always have a one-to-one correspondence to the triangles of the Delaunay tessellation. Similarly, Voronoi regions correspond to Voronoi points which are the vertices of the Delaunay triangles. What follows from this principle of duality is a one-to-one relationship between the edges of both types of 3 structures. The corresponding edges of a dual pair of diagrams do not necessarily intersect, but are always orthogonal. Moreover, the edges of the Delaunay tessellation are bisected by the edges (or extensions thereof) of the Voronoi diagram. Figure 1.3 superimposes both structures to illustrate their geometric duality. Figure 1.3 The dual relationship of the Voronoi and Delaunay structures. 1.4 Formal Definitions With an intuitive description of Voronoi diagrams and Delaunay tessellations given, we can now state the formal mathematical definitions. Voronoi Diagram Let P = {p, ,p2,..., p} denote a set of two or more distinct points in the Euclidean plane. These points are the sites of the Voronoi diagram. Now partition the plane so that every point in the plane is assigned to its nearest site with respect to the Euclidean metric. The result tessellates the plane into regions, each of which is associated with a site. The set of points assigned to site pt form the Voronoi region V(pJ. The points of F(p,) are closer to pt than any other site in P. 4 Stated formally we have V(pi) = | ||p, x|| s \\pj x|| for all j*i} (1.1) representing the Voronoi polygon (or region) associated with site pt. The set of all polygons given by y(P) = {V(PV{p2),F(pJ} (1.2) is the Voronoi diagram generated by P. Note that Voronoi polygons are defined as a closed sets (we used <, not <, in equation 1.1). As a closed set, a Voronoi polygon contains its boundary which consists of Voronoi edges. The end point of a Voronoi edge is a Voronoi vertex. Voronoi vertices are always the intersection of three or more Voronoi edges. If a Voronoi diagram contains at least one vertex at which four or more edges meet, it is said that 'V{F) is degenerate, otherwise ^V(P) is non-degenerate. Degenerate Voronoi diagrams often occur when the set of generating sites are spaced or arranged in a regular fashion. Figure 1.4 shows an example of a degenerate Voronoi diagram containing a vertex with four joining edges. Figure 1.4 A Voronoi diagram that is degenerate. 5 As we explore various properties of Voronoi diagrams, it will be convenient to assume that they are non-degenerate. This is the usual case as most of the natural occurrences of the nearest neighbor problem have non-regular arrays of points as their basis. In many applications, the degenerate cases also involve unnecessarily complex treatments which are not essential to the concepts at hand. Throughout this paper, we will assume that every Voronoi vertex in a Voronoi diagram has exactly three Voronoi edges, thereby guaranteeing non-degeneracy. Before stating our definition of the Delaunay tessellation (specifically in this case, a Delaunay triangulation), two further assumptions required. For a set of points P = {p,, p2,..., /?}c m we may assume that the points Pi ,p2, , p are not collinear. Additionally, we assume that 3 < n, as without at least three points, we cannot obtain a triangulation. Delaunay Triangulation Let 'V (P) have a generating set of distinct sites P = {p, ,p2, ...,p}, 3 < n and not all P, collinear. Construct line segments between all pairs of generators whose Voronoi polygons share a common edge. The resulting structure is a tessellation. If this tessellation is made up entirely of triangles, we have a Delaunay triangulation. The set given by V(P) = {Tj, T2,..., Tn}, 1 If some regions of D(P) are polygons with greater than three sides, a Delaunay pre-triangulation existsa configuration that results when the corresponding Voronoi diagram is degenerate. If T)(P) is not a pre-triangulation, we sometimes refer to it as a proper triangulation. By partitioning the non-triangular polygons into triangles such that no line segments from vertex to vertex intersect with a certain condition imposed (i.e., the min-max angle criterion described in Chapter 2), the pre-triangulation can be transformed into a Delaunay triangulation. As in the case of Voronoi polygons, Delaunay triangles are defined as closed sets which contain the boundary of edges. These edges are known as Delaunay edges and specifically, if a Delaunay edge is shared by two Delaunay triangles it is said to be an internal Delaunay edge; otherwise it is said to be an external Delaunay edge. 6 Note that a Delaunay triangulation must consist of at least one triangle, a condition which is given by the assumption of non-collinearity and the requirement that the structure has at least three points. It is also worth noting that the boundary of the convex hull of the generating set P is equivalent to the set of external Delaunay edges of the Delaunay triangulation of 'V (P). 7 2. Properties of Voronoi Diagrams and Delaunay Tessellations Voronoi diagrams and Delaunay triangulations enjoy many useful geometric and graph-theoretic properties. These properties are essentially exploited by a number of problems that rely on these complex structures for a solution. A specific application may utilize one or more properties which happen to be well-suited to representing and providing an answer the problem at hand. Proficiency in applying Voronoi and Delaunay structures to real-world problems requires not only a thorough understanding of this catalog of properties, but also an awareness of the relationships between a Voronoi diagram and its corresponding Delaunay triangulation. The list of properties that follow is not exhaustive, but is representative of some of the more common, and useful ones. For a thorough description of Voronoi properties with derivations, the reader is referred to [65] from which many of the following properties were gathered. 2.1 Properties of the Voronoi Diagram The following properties deal mainly with Voronoi Diagrams in the planar case, although some properties can be extended to include diagrams in Rn. It should be noted also that some of the properties may not hold for generalized Voronoi diagrams (see Chapter 5). We begin our discussion of properties with an observation about Voronoi regions and how they collectively tessellate the real plane. VI: Each Voronoi region V(p,) is a non-empty convex polygon. Let P = {p, ,p2 ...,/?}c M2, 2 < n be a set of distinct points. Then the set V(p,) defined by V(p,) = {xeR2 | ||pt x|l < ||pj x|| for all j*i, j e I} (2.1) is a non-empty convex polygon. V2: The Voronoi diagram 'V (P) is a unique tessellation for the set P of Voronoi points. 8 Our Voronoi diagram defined as V(i>) {V(p,), V(p2),F(p)} (2.2) satisfies UieF-(A) = R2, (2.3) [F(p,)W(p,)] n [F(p,)W(p)] = e>, ]*U hi el, (2.4) which is to say that the set of Voronoi polygons uniquely tessellate the real plane. This property is also valid if polygons are replaced by polyhedons (or polytopes) and R2 is replaced by Rm Since finite Voronoi polytopes are a convex sets, they are sometimes referred to as Voronoi polytopes. It is evident from Figure 1.1 that some regions of V(P) are bounded while others are infinite. The following property formally identifies the distinction. VS: A Voronoi polygon V(p,) is infinite (unbounded) if and only if pt is on the boundary of the convex hull of the point set P. That is, pt e 5CH(P). For a given Voronoi diagram V(P), it is easy to construct the convex hull CH(P) of P by simply joining the Voronoi points whose Voronoi polygons are unbounded. Property V3 also holds for a Voronoi polytope in Mm . V4: Given a Voronoi diagram constructed from a set of distinct generators P= {p, ,p2, R2, 2 < n the following related properties hold: (1) The Voronoi edges of "V(P) are infinite lines if and only if P is collinear set. (2) A Voronoi edge, denoted by e(p{, pj), is a half line if and only if P is a non-collinear set and pt and p} are adjacent generators of the boundary of CH(P). (3) If Pi and pj generate a Voronoi edge e(p,, pj), then this edge is a finite line segment if and only if P is a non-collinear set and at least one ofPi, Pj is in the interior (i.e., not on the boundary) of CH(P). 9 Note that the Voronoi edges of 'V(P) are segments of bisectors, but bisectors do not always produce Voronoi edges. The following property gives the sufficient condition for the generation of a Voronoi edge. VS: The generator point closest to point pt generates a Voronoi edge of V(p,). Property V6 follows as an immediate consequence of property V5. V6: The generator point nearest to point pt is one of the generators whose Voronoi polygons share a Voronoi edge with polygon V(p^. Property V6 is very useful in solving several well-known proximity problems, two of which are now stated. Closest Pair Problem For a set of distinct points P, if the distance from p} e P to p} e P is the minimum distance between all pairs of points in P, then the pair (p,, pj} is the closest pair among all points in P. Given P, find the closest pair. All Nearest Neighbor Problem Given a set of distinct points P, find the nearest neighbor point for every Pi e P. In other words, find all nearest neighbors for a given set of points. For a complete description of the some of the more common proximity problems see [67] and [72]. The nearest generator point from another generator point can be solved by employing properties V5 and V6. Problems dealing with finding the closest generator point from an arbitrary point p is addressed by the next property. V7: Point pf is the nearest generator point from an arbitrary point p if and only if Voronoi region F(p,) contains p. This property, although basically the definition of a Voronoi polygon, is useful in solving the following problem. Nearest Search Problem If we are given a set of distinct points P, determine the nearest neighbor point in P from a given point p, with the condition that p is not necessarily contained in the set P. 10 A brute-force method of solving this problem is to find the minimum distance among all the possible n distances, {\\p pt\, i e /}. However, for a large number of probe points this computation becomes burdensome. Property V7 significantly reduces the computational load by first finding the Voronoi region in which the probe point is located. When this polygon V(p,) is found, property V7 immediately tells us that the associated generator point pt is the nearest generator point. The vertices of a Voronoi diagram also have a number of useful properties, some of which are now stated. It should first be noted that for a given set of points P, a circle C, centered at vertex q, is referred to as empty if the interior of C, contains no other points of P. V8: Let Q- {q,,q2,..., q) be the set of Voronoi vertices of a Voronoi diagram "VfP). For each vertex qt e Q, there exists a unique empty circle C, with its center at <7, which passes through three or more generators of y(P). Note that by our previously stated non-degeneracy assumption, the circumference of C, will contain exactly three generator points of 'V(P). If this assumption does not hold, lengthy, yet inconsequential details must be added to the derivations and proofs of the associated theorems. Therefore, throughout this paper we make the following assumption. .\on -cocircularity Assumption For a set of Voronoi points P = {pj ,p2,...,p}c R2, 4 ^ n, a circle C does not exist such that pu ,pi2,pik e P, k > 4, are on the boundary of C. It follows that every point in P \ {pu pa,..., p,k) are exterior to C. If circle C is replaced with a sphere (or hypersphere) and k > m + 2, then the non- cocircularity assumption holds in Rm as well. Figure 2.1 demonstrates the degeneracy of a Voronoi diagram by illustrating the cocircularity of four points. 11 Figure 2.1 (a) The non-cocircularity of four points; (b) four cocircular points The following property is an obvious consequence of the non-cocircularity assumption. V9: Every vertex q, of a Voronoi diagram V(P) is the common intersection of exactly three edges of "V(P). The next problem occurs quite often in a number geometric applications and is motivated by the property which follows it. Largest Empty Circle Problem Given P, a set of distinct points in R2, find the largest empty circle whose center in the interior of CH(P). VI0: Let Q = {q,, q2,..., q) be the set of Voronoi vertices of a Voronoi diagram ~V(P). A circle C: centered at vertex g, and passing through three or more generators of V(P) is the largest empty circle centered at <7,. An example of the largest empty circle problem for nine points is shown in Figure 2.2. The application of this problem is useful in solving collision avoidance problems and a variety of spatial-temporal processes [65]. 12 ORourke [67] gives an intuitive property dealing with circles in Voronoi diagrams that follows nicely from V9 and VI0 above. Vll: Let qt be a Voronoi vertex at the junction of V(p,), V(p2), and V(p3). Then <7, is the center of the circle C, determined by p, ,p2, and p3. To continue with this reasoning, ORourke [67] states another property concerning circles, but this time relating them to Delaunay triangles. V12: Let q, be a Voronoi vertex at the junction of V(p,), V(p2), and V(p3) with the resulting circle C, centered at qt. Then C, is the circumcircle of the Delaunay triangle corresponding to q,. This obvious result stems from a fundamental theorem of plane geometry: that the perpendicular bisectors of the sides of any triangle intersect at the center of the circle that circumscribes the triangle (i.e., the circumcenter). 13 A number of problems are concerned with determining if a given planar tessellation is in fact a Voronoi diagram. This is stated as the following problem. Generator Recognition Problem Given the set of Voronoi edges for a non-degenerate Voronoi diagram 'V(P), find the locations of generators P. A practical technique for locating the set of generators is given in [65]. This reference also provides a method of determining, more generally, whether any given planar tessellation is a Voronoi diagram. One final property of Voronoi diagrams (which is also a theorem) provides a transition into the next section. VI3: The straight-line dual of a Voronoi diagram "V(P) is a unique triangulation of P. Furthermore, it is a Delaunay triangulation. This property is somewhat obvious when one observes the superimposition of the structures (as in Figure 1.3). However, the proof of this theorem, which was first given by Delaunay [22] in 1934, is rather lengthy and omitted here. Preperata and Shamos [72] provide the details of the derivation. 2.2 Properties of the Delaunay Tessellation As in the case of Voronoi diagrams, many of these properties can be extended to dimensions higher than l2. We begin our discussion of Delaunay properties with several general attributes. Dl: The Delaunay region T,, defined above, is a unique non-empty polygon and the set T)(P) = {T,, T2,..., Tn) satisfies U(ai T, = CH(P) , (2.5) [(WI3]n[(9(15]-0, j*i, ijel,. (2.6) The implication here is that T>(P) is equivalent to the convex hull of the set P. Stated more precisely, the external edges of a Delaunay tessellation constitute the 14 boundary of the convex hull of the generating sites. Furthermore, D(P) is a unique tessellation over P. One may question under what circumstances is a Delaunay tessellation in R2 a triangulation. Recall that D(P) is not a proper triangulation if its dual 'V(P) is degenerate. ORourke [67] provides the following criterion which he claims is Delaunays Theorem. D2: T){P) is a triangulation if there are no four cocircular points in P. Consequently, all polygons of D(P) are triangles. The next property relates the point sets and the cardinality of the edge sets for any pair of dual structures. D3: Let: y(P) and D(P) be generated by P, Q be the set of Voronoi vertices, E be the set of Voronoi edges, Qd be the set of Delaunay vertices, Ed be the set of Delaunay edges, and Cd be the set of circumcenters of the Delaunay triangles. Then the following hold: (1) Qd = P- (2) ||Ed || > ||E ||, with equality if and only if "V(P) is non- degenerate. (3) Cd=Q. For Delaunay tessellations in the plane, we may question whether the triangles are empty. The following property confirms that they are. D4: For any Delaunay triangulation T>{P), the interior of every triangle in T>{P) contains no points of P. A natural extension to this result concerns the emptiness of the circumcircles about the triangles. By noting the relations in D3 and recalling property V8, we derive another attribute of Delaunay triangulations which can be of significant use in 15 solving the largest empty circle problem (see VI0 above). Figure 2.3 provides an example of this property for a Delaunay triangulation on 11 points. D5: For a Delaunay triangulation D(P) generated by the set P, all circumcircles of the resulting Delaunay triangles are empty circles. Figure 2.3 The circumcircles of all Delaunay triangles are empty circles. Our next property presents a criterion for the inclusion of an edge e in the set of edges Ed of TKJP). The derivation can be found in [67]. D6: For a given T>{P), if there is a circle Cd which passes through two points Pi and pj such that Cd contains no other p, then the line ip,, p]) is an edge of D(P). The converse is also true; for every edge e of T)(P), there exists an empty circle for which e is a chord. 16 The next property is rather intuitive, but also worthy of a formal statement. D7: Let D(P) be a Delaunay triangulation. For a given point p, e P, if the nearest neighbor to pt is point pj, then the line pi ,pj is an edge of D(P). We direct our attention next to a class of properties that deal with the magnitude of the internal angles of the triangles of Delaunay structures. This is an important area of study because of the need to avoid elongated triangles in certain applications such as constructing surface interpolations and finite element meshes. Fortune [38] gives a description of the angular characteristics of Delaunay structures in two and three dimensions. The reader is referred to a paper by Sibson [83] which discusses some of the angular characteristics of triangulations including algorithms for the construction of equiangular Delaunay structures. We present the following observation for T)(P) in the plane. D8: For every triangulation spanning a set of points P fR2, the Delaunay triangulation Tf(P) maximizes the minimum angle of any triangle. A generalization of this property does not exist for higher dimensions, although an elegant treatment of this Delaunay property in R3 is provided in [38]. Both Fortune and Aurenhammer [6] attribute this optimality criterion to Rajan [73]. For the next property, we first note that a containment sphere is the smallest sphere that entirely encloses a polyhedron. D9: Over all tetrahedral tessellations in R3, the Delaunay tessellations minimize the largest containment sphere for all tetrahedra in the tessellation. One implication of D8 and D9 is that we may think of Delaunay tessellations as being optimum structures in some sense. Aurenhammer [6] alludes to the fact that a Delaunay structure is the most compact of all tessellations. Articles by Fang and Piegl [32] and Watson [95] are two notable examples of the many articles on Delaunay tessellations in higher dimensions. 2.3 Graph-theoretic Properties of 'V(P) and D(P) Voronoi diagrams and Delaunay tessellations in R2 have a number of characteristics which are analogous to those of planar graphs. In feet, they can be represented as graphs. This alternate perspective has been particularly useful in 17 describing the characteristics of networks and for solving path planning problems. With many applications requiring this approach, considerable study has been given to the graph-theoretic properties of these geometric structures. The topological attributes of Voronoi and Delaunay graphs are the incidence relationships among the vertices, edges, and regions of the structure. Metric parameters such as the coordinates of vertices and the lengths of edges in a graph are also useful in many applications. These incidence and metric properties will be discussed in the next chapter where we will need to represent this information to store Voronoi diagrams in a data structure. The properties that follow provide bounds on the quantities of vertices, edges, and regions with respect to the overall structure. For additional insight into the graph- theoretic properties of Delaunay graphs, we examine how they are uniquely related to several other notable graphs including convex hulls, relative neighborhood graphs, and minimum spanning trees. 2.3.1 Abstract and Geometric Graphs Graphs are normally used to represent the topological relationships of a set of points and the edges joining them. The essential information in a graph is whether an edge (or in the case of directed graphs, two directed edges) exists between a pair of points. We begin by formalizing the definition of an abstract graph. Abstract Graph If we let P = {pj, p2,..., p) be a non-empty set, and {p,, py} be an unordered pair such that p} ,p2 e P, then the set of all distinct unordered pairs of points {iPi > Pj} I Pi 5 Pj Â£ P) is called the unordered product of P. This unordered product is denoted by P P. Also let L= {L{ ,L2,..., L} and/be a mapping of L into P % P, which defines the function/ (/) = {p,, py}. The paired sets P and L together with the function/form a graph, or more precisely an abstract graph. We use the common notation G(P, L,f) to represent an abstract graph. Sometimes, G(P, L) or simply G will also be used to denote an abstract graph. It is often useful to represent an abstract graph by a set of geometric points and the set of geometric line segments joining the points. To do this, we introduce the notion of a geometric graph. Geometric Graph Let Pg = (p/, p2,..., p} c lm be a set of n distinct points in m-dimensional Euclidean space. Also let Lg = {L,, L2,..., L) denote a set of non-self- 18 intersecting lines segments with end points in Pg. The line segments should satisfy the following three conditions: (1) Each line segment in Lg that forms a loop contains exactly one point of Pg. (By a loop, we mean a line segment L, that has both end points coinciding at one point p,, and Z, also does not intersect itself.) (2) Each line segment in Lg that does not form a loop contains exactly two points of Pg that are the end points of the line segment. (3) Each line segment in Lg does not have any common points with the exception of points in Pg. If these three conditions hold, we call the paired sets (Pg, Lg) a geometric graph and denote it as G(Pg, Lg). We refer to an isomorphism between an abstract graph G(P, L) and a geometric graph G(Pg, Zg) as a geometric realization of G(P, L). Fary [33] proved that any planar abstract graph can be embedded in the Euclidean plane where all edges of G(P, L) map to straight line segments of G(Pg, Lg) (i.e., a geometric realization). The concept of geometric graphs will be helpful in Chapter 3.1 where we define a data structure that is used to represent Voronoi diagrams. 2.3.2 Eulers Formula As a prerequisite to the properties that follow, we begin with Eulers celebrated formula for graphs and its analog for polyhedra. In 1758, Leonard Euler discovered a significant relationship between the numbers of vertices, edges, and regions for every planar graph: the number of edges is always two less than the number of regions and vertices combined. Letting nv,ne, and nr be the number of vertices, edges, and regions respectively, we have nv-ne + nr = 2, (2.7) which is Eulers formula. Eulers formula can be proved in a great variety of ways. A concise proof by Coxeter may be found in [19], while a much older treatment of the proof is included in the classic Mathematical Recreations and Essays by W.W. Rouse Ball [9]. 19 Euler also showed that every convex polyhedron can be represented as a planar graph. Thus, if we let V, E and F be the number of vertices, edges, and faces respectively, we now have the equation V-E + F= 2, (2.8) which is valid for any convex polyhedron. Around 1852, Eulers formula for polyhedra was generalized to any dimension by Schlaefli. The subsequent extension to the original relation is known as the Euler- Schlaefli formula and is formalized as follows. Let fT be a tessellation of a finite set P in Mm, and let n, be the number of /-dimensional faces in eT. We then have Â£(-iy,=i+(-ir (2.9) i = 0 valid for polytopes of any dimension. A number of other topological extensions to Eulers formula, including the history of their discovery, can be found in the delightful book Proofs and Refutations by Lakatos [51]. 2.3.3 Voronoi Graphs To begin exploring this subject, we first demonstrate how to derive a finite planar graph from a two-dimensional Voronoi diagram. Voronoi Graph Derivation Let Q= {q}, q2,..., q} be the vertices of V(P) and let E = {e,, e2,..., e} be the set of Voronoi edges in IfP) such that the first nc edges are infinite. As we know, infinite edges are not permitted in a finite graph. To allow for this, simply place a new point, say q0 at a sufficiently great distance from the convex hull of P. Each infinite edge of 'V(P) is then clipped at a certain point and its end extended to join with q0. This new set of modified edges along with the unmodified edges combine to form the set of edges Eb of the graph. Similarly, Q+I = QU {q0} forms the vertex set for the graph. In standard graph theory notation, we now have G(Q+I, Eb), which is the Voronoi graph induced from 'V(P). Figure 2.4 shows a Voronoi graph derived from a Voronoi diagram. 20 Figure 2.4 A Voronoi graph derived from a Voronoi diagram. VG1: Let nv,ne, and ng represent the number of vertices, edges, and generators (or regions) respectively in 'V(P) c M2, 2 < ng, then nv- ne + ng = 1 . (2.10) The deviation from the usual Euler formula (i.e., this sum equals 1 instead of 2) is a result of the additional vertex added to 'V(P) in order to make it a graph. As an extension of this relationand one that is analogous to the generalization of Eulers formulawe give the following property. VG2: Let nk represent the number of ^-dimensional Voronoi regions in an /w-dimensional "V(P). By using the Euler-Schlaefli formula (equation 2.9) we have I k = 0 (-l)*nA = (-ir. (2-11) If the Voronoi diagram that induces a Voronoi graph G(Q+,, Eb) is non- degenerate and has at least three generators, then each vertex of the graph will 21 have at least three edges joining it. By a simple counting argument, the number of edges in G(Q+I, Eb) is ne > 3 (nv + 1) / 2 . (2.12) By substituting this inequality into the Euler formula for Voronoi graphs (nv ne + ng =1) we derive our next property. VG3: Let nv,ne, and ng represent the number of vertices, edges, and generators respectively in 'V(P) cl2,3 nv< 2 ng 5 . (2.14) For an elegant derivation of this property, see [72]. The next couple of properties are derived by continuing counting arguments on preceding relations. VG4: Let nv,ne,ng, and nc represent the number of vertices, edges, generators, and unbounded polygons respectively in "V(P) c M2, 3 < ng. The following inequalities are valid: ne > 3 nv + nc 3 , (2.15) nv > 1/2 (ng nc) + 1 . (2.16) The next result can be derived directly from equation 2.13 above. VG5: For every planar Voronoi diagram, the average number of Voronoi edges per Voronoi region is not greater than six. This result does not generalize to higher dimensions. However, the upper and lower bounds on the maximum number of Voronoi vertices for V(P) in W are known. These rather complex inequalities and lengthy proofs may be found in [50]. 22 2.3.4 Delaunay Graphs Delaunay Graph Derivation A planar graph can be derived from a Delaunay triangulation quite simply as the structure itself is equivalent to a connected geometric graph. A Delaunay triangulation can be regarded as Delaunay graph whose nodes are given by the generating set P and set of edges, say Ed, given by the edges of the Delaunay triangles. In standard notation we have G(P, Ed) representing a Delaunay graph. Our first property of Delaunay graphs is an obvious one, but should not go without being stated. DG1: Every Delaunay graph G(P, Ed) is a dual of a Yoronoi graph G(Q+;, Eb). A property describing relationships between the number of Delaunay vertices, edges, external edges, and triangles can be gleaned from properties D2, DG1, and the fact that the external edges of T>{P) have a one-to-one correspondence with infinite Voronoi polygons of the dual 'V(F). DG2: Let T>(P) be a proper Delaunay triangulation over the set P. Also let ne be the number of Delaunay (or Voronoi) edges, nt be the number of triangles in D(P), nv be the number of vertices in D(P), and nc be the number of vertices on the boundary of CH(P). The following two equations will hold: nt = 2nv nc 2, (2.17) ne = 3 nv nc 3 . (2.18) Several interesting properties are known concerning the isomorphism of Delaunay triangulations to general triangulations constructed on P. The details may be found in [24], 2.3.5 Subgraphs of the Delaunay Triangulation We conclude this chapter with a few properties that describe the relationship of Delaunay triangulations to several other well-known graphs. The Euclidean minimum spanning tree (EMST), the relative neighborhood graph (RNG), and the Gabriel graph (GG) all have far-ranging applications. As we shall see, these graphs have a containment relationship to one another (i.e., there is a nesting of them) 23 and remarkably, they are all subgraphs of the Delaunay triangulation. For an excellent background reference on these notable graphs (and graph theory at large) see Applied Combinatorics by Roberts [74]. References for the specific graphs are given in Chapter 5.1 where we discuss their applications. The Gabriel Graph The Gabriel graph, denoted GG(P), for a set of points P is the graph in which the line pfPj is an edge if and only if the circle C0, defined by diameter p,pj, is an empty circle. Clearly, from properties D3 and D4 we notice the following relationship. DG3: Every edge in GG(P) is a Delaunay edge of T)(P). The necessary and sufficient condition for an edge p,pj of GG(P) to be an edge of T>(P), is that edge p,p} intersects the corresponding Voronoi edge. Figure 2.5 (b) shows the Gabriel graph for the same set of points as the Delaunay triangulation in Figure 2.5 (a). Relative Neighborhood Graph Denoted by RNG(P), the relative neighborhood graph is defined such that RNG(P) has an edge joining p, and pj if and only if dip,, pj) < min max {d(p,, pk), d(pj, pk)} , (2.19) k i ,j with d representing the Euclidean distance function. What equation 2.19 says is that, edge PiPj c RNG(P) if and only if there is no other point of P in the interior of the intersection of circles C, and C, both having radius d(pj, p). From this condition we also observe that an edge of RNG(P) is also an edge GG(P). So we have another property. DG4: All edges of RNG(P) are also edges of T>(P). The relative neighborhood graph on the set of points P of D(P) is presented in Figure 2.5 (c). 24 Euclidean Minimum Spanning Tree Represented by EMST(P), the Euclidean minimum spanning tree is defined as a tree graph with vertices in P, such that the sum of the Euclidean lengths of all the edges ofEMST(P) is the minimum of all possible trees which can be made from the vertices of P. Continuing the hierarchy, Okabe et al. [65] show that the edges ofEMST(P) are contained within the set of edges of RNG(P). Thus, we obtain the following property. DG5: The edges of EMST(P) are contained in the edge set of T>(P). Figure 2.5 (d) shows the EMST for the points P of T)(P). Figure 2.5 (a) Delaunay triangulation; (b) Gabriel graph; (c) relative neighborhood graph; (d) Euclidean minimum spanning tree; (e) convex hull. 25 In summary, for a set of points P in the Euclidean plane, the graphs GG(P), RNG(P), and EMST(P) are all subgraphs of T>(P). Furthermore, as implied in property D1 and shown in Figure 2.5 (e), the boundary of the convex hull dCH(P) is also a subgraph of T)(P). Figure 2.6 shows the containment hierarchy of these five graphs. Figure 2,6 The Delaunay triangulation and its subgraphs. In Chapter 4 we will show how Delaunay triangulations can be constructed by using any one of a variety of algorithms. The relationships stated above provide an effective method for constructing these graphs by first producing T)(P). The technique is straightforward; begin by constructing D(P) and remove the extraneous edges to produce the desired graph. 26 3. Representing Voronoi Diagrams in a Data Structure Before considering practical algorithms for the construction of Voronoi diagrams which follow in Chapter 4, we must first decide on a data structure which can adequately represent a Voronoi diagram. Moreover, we have to initially clarify exactly what is meant by the construction of the diagram. At first one may think that representing a Voronoi diagram is equivalent to drawing it on a sheet of paper or on a computer screen. These characterizations are visually intuitive and do reveal some of the general features about a diagrams structure. However, the numerical data required to draw the diagram does not necessarily contain enough information for this representation to be useful in a particular application. On the other hand, some representations can be overly detailed. Listing each and every edge and vertex on the boundary of every polygon in a Voronoi diagram requires a great deal of data storage space. Furthermore, it still may not convey all the data needed to fully describe the diagram What we need is a data structure that is concise, yet contains sufficiently detailed information that is easily retrievable. Examples of what may be required from such a data structure include: List the order of polygons and edges surrounding a particular vertex Retrieve the edge that is shared by two given polygons Determine which edges also join the two vertices of a particular edge 3.1 Winged-edge Data Structure We now describe a data structure that is quite suitable for representing Voronoi diagrams. It meets our requirements of sufficient detail, avoidance of unnecessarily large storage space, and ease of information retrieval. A thorough treatment of the concepts presented here may be found in [65] although the original work on the subject was first published in [11]. The primary function of the winged-edge data structure is to numerically represent planar geometric graphs. Also known as the polygon data structure, this form of geometric representation is widely-used in computer graphics and geometric modeling [10] and [46]. The winged-edge data structure holds explicit information concerning local incidence relations about the vertices, edges, and polygons of a Voronoi diagram (or any geometric graph). This data structure supports the rapid 27 retrieval of a geometric graphs particular characteristics. At the same time, it is not overly intensive on data storage requirements. As we mentioned above, the winged-edge data structure is used for storing information about planar geometric graphs. Recall that infinite edges are not permitted in geometric graphs. Voronoi diagrams however, do have edges that extend to infinity. As in the Voronoi graph derivation given in Section 2.3.3, we must bound these infinite edges, but without introducing an additional Voronoi vertex. Although it served our purposes for the graph-theoretic properties presented in Chapter 2, introducing an additional vertex here would not accurately characterize the Voronoi diagram we want to represent. To overcome this we introduce a closed curve which is sufficiently large to encompasses all the vertices qt and all the generating points p, of V(P). If we consider the infinite edges of y(P) terminating at this curve, as shown in Figure 3.1, then this bounding curve is divided into line segments We will consider these curved line segments among the set of Voronoi edges of the resulting geometric graph. We may also think of the vertices on this curve as points at infinity. Okabe et al. [65] calls this derived geometric graph the augmented geometric graph associated with V(P). Figure 3.1 Augmented geometric graph of a Voronoi diagram The winged-edge data structure representing an augmented geometric graph is now produced as follows. Convert the geometric graph to a directed geometric graph by assigning a direction to each edge. As you will see, this may be done 28 arbitrarily as the directions do not affect the nature of the information stored by the data structure. Next, the vertices of the augmented geometric graph are assigned a fixed linear order from 1 to nv. Likewise, the edges of this graph are given an ordering from 1 to ne. For the Voronoi diagram 'V(P) generated by set P = {pj, p2,..., p}, the corresponding augmented geometric graph will have n+ 1 regions. The additional polygon in this set is the outermost region (i.e., the region outside the bounding closed curve just introduced). This can polygon be considered the result of an added virtual generator , which gives the augmented geometric graph the generating set P* = {p,, p2,..., p, />}. The labeling for the polygons is straightforward. For a polygon p, of 'V(P), the associated augmented geometric graph polygon will be named polygon i, where i = 1,2,..., n, . With a labeled, directed geometric graph in place, we can now represent the incidence relations among the graphs vertices, edges, and polygons. To achieve this, we set up ten 1-dimensional arrays of integers. As we shall see, one array is required for the polygons, one array for the vertices, and eight arrays for the edges. These ten incidence arrays (numbered 1 through 10 below) contain the topological characteristics of ~V(F). They are defined in the following manner. For polygon /, (/ = 1, 2,..., n, <): (1) EdgeAroundPolygonjV] is the ordinal number of an edge on the boundary of polygon i. For vertex j, (j= 1,2,..., nv): (2) EdgeAroundVertex[/] is the ordinal number of an edge incident to vertex j. The edge arrays contain the majority of information about the graph. Figure 3.2 shows the eight geometric objects that are incident and adjacent to each edge in the graph. These eight objects are represented by arrays numbered 3 through 10. 29 CCWSuccessor[&] CWSuccessor[&] CWPredecessor[Â£] CCWPredecessor[Â£] Figure 3.2 The winged edge and its eight associated geometric objects. For edge k, (Â£= 1, 2,ne): (3) RightPolygon[Â£] is the ordinal number of the polygon to the right of edge k, (4) LeftPolygon[&] is the ordinal number of the polygon to the left of edge k, (5) StartVertex[Â£] is the ordinal number of the start vertex of edge k, (6) EndVertex[fc] is the ordinal number of the end vertex of edge k, (7) CWPredecessor [&] is the ordinal number of the edge adjacent to edge k clockwise around the start vertex, (8) CCWPredecessor[Â£] is the ordinal number of the edge adjacent to edge k counterclockwise around the start vertex, (9) CWSuccessor[Â£] is the ordinal number of the edge adjacent to edge k clockwise around the end vertex, 30 (10) CCWSuccessor[Â£] is the ordinal number of the edge adjacent to edge k counterclockwise around the end vertex. These ten arrays collectively represent the topological features of V(P). However, to completely characterize a Voronoi diagram, we also need to capture its metric qualities. For example, to be useful in a network application we may need to know the length of the Voronoi edges and perhaps the Cartesian coordinates of the set of generators. For this purpose, three additional arrays (numbered 11 through 13) are introduced to define the metric parameters of the Voronoi diagram. (11) TypeVertexD'] = 0 if vertex j is an ordinary point. TypeVertex[/] = 1 if vertex j is a point at infinity. This array holds binary information which determines the contents of the next two arrays. (12) X-VertexD']. (13) Y-Vert ex[/]. Both of these arrays may each contain two types of data as designated by array 11. The two types of data are: (a) If TypeVertex[/] = 0, then vertex j is an ordinary point and X-Vertex[/'], Y-Vertex[/'] will represent the x and y coordinates of this vertex. (a) If TypeVertexD-] = 1, then vertex j is a point at infinity and X-Vertex[/], Y-Vertex[/] will represent the x andy components of the unit vector that designates the direction of the associated (infinite) Voronoi edge. This set of thirteen arrays is collectively called the winged-edge data structure. In the following chapter on algorithms we will refer to the construction of a Voronoi diagram. This should be taken to mean the construction of the winged- edge data structure for the augmented geometric graph associated with V(P). 31 3.2 Geometric Information Retrieval With a winged-edge data structure in place, we can easily retrieve basic information about ^V{P). One of the most common requirements is to produce a list of edges and polygons that join at a particular vertex. The following algorithm accomplishes this. Algorithm 3.1: Retrieve geometric objects incident to a vertex Input: Winged-edge data structure of 'V{F) and vertex j. Output: Edge list Le and polygon list Lp that surround vertex j in a counterclockwise manner. Procedure: Step 1. Le empty list and Lp<- empty list. Step 2. k- EdgeAroundVertex[/'] and kStart k. Step 3. Add k to the tail of the list Le. Step 4. If/ = Start Vertex^] then add LeftPolygon[Â£] to the tail of Lp, and k- CCWPredecessor[Â£], else add RightPolygon[&] to the tail of Lp, and k CCWSuccessor[Â£]. Step 5. If k = kStart then return Le and Lp, else go to Step 3. In a single run of this algorithm, Steps 1 and 2 are performed only once. Steps 3, 4, and 5 are repeated as many times as there are edges (or polygons) incident to vertex j. The performance time of this algorithm is proportional to the degree of vertex j (which is always three if the Voronoi diagram is non-degenerate). Consequently, we can expect a 0(n) or linear performance from this algorithm. Two other useful sets of information about "V(P) are the list of edges and the list of vertices that surround a given polygon. Our next algorithm returns these lists. Algorithm 3.2: Retrieve the geometric objects surrounding a polygon Input: Winged-edge data structure of 'V(Pf) and polygon i. Output: Edge list Le and vertex list Lv that surround polygon i in a counterclockwise manner. Procedure: Step 1. Le empty list and Lv empty list. Step 2. k- EdgeAroundPolygon[z] and kStart k. 32 Step 3. Add k to the tail of the list Le. Step 4. If i LeftPolygon[&] then add EndVertex[fc] to the tail of Lv, and k *- CWSuccessor[Â£], else add StartVertex[&] to the tail of Lv, and k CWPredecessor[&]. Step 5. If k = kStart then return Le and Lv, else go to Step 3. In a single run of this algorithm, Steps 1 and 2 are done only once. Steps 3,4, and 5 are repeated as many times as there are edges on the boundary of polygon i. The performance time of this algorithm is proportional to the number of edges on the boundary of polygon i. Therefore, the performance of Algorithm 3.2 is also O(n). Take note that if polygon is used as input to Algorithm 3.2, the list of edges and vertices on the boundary of this infinite region is returned. The boundary of this infinite polygon is actually the convex hull CH(P) of the generating points of "V(P). This leads us to question whether we can retrieve a Voronoi diagrams dual from the data structure. It turns out that the winged-edge data structure is also convenient for retrieving a Voronoi diagrams corresponding Delaunay triangulation. To retrieve information about the Delaunay triangulation we must rename the geometric objects with their dual counterparts. To achieve this, simply exchange the roles of the vertices and polygons. This is, after all, the way in which a pair of duals are related. The reader is referred to [65] for complete details on extracting a Delaunay triangulation from a winged-edge data structure. 3.3 Complexity of Constructing V(P) It is not difficult to determine the lower bound on the time complexity for the representation, and as a result the construction, of Voronoi diagrams. This is done by a comparison to the well-known result that the lower bound for sorting n real numbers in increasing order is in the worst case 0{n log n). We show this as follows. Let X= {xy, x2,..., x} denote a set of n real numbers. Construct the set of plane coordinates P = {(x,, x,2), (x2, x/),..., (x, x2)} and build V(P). We can now obtain a cyclic sequence of points x,, x2,..., x on the boundary of CH(P). This may be done by using Algorithm 3.2 where we consider as output, the vertices of the polygon formed by P. Figure 3.3 shows the vertices of CH(P) incident to the parabolay = x2 and the corresponding points x, ,x2,..., x that are the projection 33 of the vertices of CH(P) downward in the negative y direction to x-axis. The points ofX now lie along the x-axis in increasing order. This implies that if a Voronoi diagram generated from set P could be constructed faster than 0(n log n) time, then we would be able to sort X faster than 0(n log n) time. Clearly, this contradiction shows that any algorithm for building "V(P) where ||P || = n requires at least 0(n log ri) time in the worst case. Figure 3.3 Using the convex hull to sort real numbers. If we consider the average case performance, we see that any construction algorithm can create a winged-edge data structure in only 0(n) time. Therefore, 0{ri) is the performance for constructing a Voronoi diagram in the average case. As we shall see in the next chapter, the performances given here for both the average and worst case time complexity are strict bounds on the algorithms for Voronoi diagram construction. 34 3.4 Other Data Structures Used To Represent V(jP) As a final note on representing geometric graphs, we mention two commonly-used data structures that also are used for representation in computational geometry. It will suffice here to provide only the high-level descriptions of these data structures; the interested reader may seek specific details in the papers referenced. The first is called the doubly-connected-edge list, or DCEL. Initially outlined in [60], DCEL uses eight arrays to represent the topological information about a geometric graph. This data structure also supports the efficient storage and retrieval of Voronoi diagrams with an overall performance of 0(n). A concise description of DCEL and its retrieval algorithms can be found in [72]. Another efficient method that may be used for storing Voronoi diagrams (and geometric objects in general) is the quad-edge data structure [45]. While more complex than either the winged-edge data structure or the DCEL, this representation has the advantage of simplicity when used in conjunction with some Voronoi diagram construction algorithms. Its use is somewhat analogous to the winged-edge data structure since the representation of the polygons and vertices require minimal storage space. For example, the vertices of a polygon are assigned (in the quad-edge data structure) to an arbitrary edge among the set of edges that bound the polygon defined by these vertices. Moreover, an edge pointer conveniently provides access to this list of edges. The polygons of a Voronoi diagram are also represented and retrieved in the same concise manner. If this method is used to represent a Voronoi diagram, the corresponding Delaunay triangulation is automatically encoded without incurring additional complexity; simply interpret the vertices as polygons and the polygons as vertices. For a short overview of the quad-edge data structure, see [67]. 35 4. Algorithms For Constructing Voronoi Diagrams This chapter discusses computer algorithms for the construction of Voronoi diagrams. The variety of methods that are used for creating Voronoi diagrams are varied and perhaps as fascinating as the diagrams themselves. As mentioned in the previous chapter, we will consider the construction of a Voronoi diagram as a representation in a data structure (e.g., the winged-edge data structure). For starters we present a simple algorithm that follows directly from an alternative definition of the Voronoi diagram. This naive algorithm constructs the Voronoi polygons one at a time. Although inefficient for practical purposes, it is insightful to examine a brute-force method as it can illustrate the shortcomings of an obvious approach to building a Voronoi diagram. Next, we provide the details of three well-known construction algorithms. These algorithms are consistently efficient in most applications; one has an average time complexity of 0(ri) and two have a worst-case performance of 0(n log n). They do however, vary in degree of difficulty of implementation. The three algorithms we examine are as follows. Incremental Method This popular algorithm works by successively adding new generating points to an existing Voronoi diagram. The introduction of each generator grows an additional Voronoi polygon in the diagram. This method is powerful, easy to understand, and relatively simple to implement. Divide-and-Conquer This fundamental technique is used often for designing efficient algorithms. It also adapts well to the construction of Voronoi diagrams. The divide and conquer algorithm recursively constructs subsets of the diagram and then merges the subsets together. The details of the merging process are rather involved. Consequently, this algorithm is somewhat difficult to implement. Plane Sweep The plane sweep is a clever technique used to solve many two-dimensional geometric problems. A sweep line is passed across the plane (e.g., from left to right) whereby the Voronoi diagram is constructed along this line. This algorithm 36 is simple in principle, although a mapping of the original Voronoi diagram to a transformation must precede the sweep. The difficulty of implementation falls somewhere between that of the incremental method and the divide-and-conquer technique. Following each construction method, we analyze the algorithmic performance. The average and worst-case performances are derived from the subsequent steps of each algorithm. A table of comparative performance summarizes these results. We continue our discussion of construction algorithms with a look at implementation issues and the ways in which we can produce a numerically valid computer program. Real-world applications often involve data sets which can only be accurately represented by floating-point numbers. With computers capable of a limited finite-precision arithmetic, geometric algorithms are very susceptible to inconsistencies caused by numerical errors. The algorithms presented have several common pitfalls of numerical imprecision which must be avoided to produce a valid construction. Initially, we make a couple of assumptions to help avert constructional instability. Later in this chapter we examine implementation techniques that allow us to abandon these simplifying assumptions and thereby produce valid constructions in a real-world situation. The two assumptions now imposed are: Empty Circle Assumption To avoid degeneracy, we assume that no four Voronoi generating points he on a common circle. Recall that our previous non-cocircularity assumption (in Chapter 2) stated that there are no empty circles which pass through four generator points. The assumption here is more general; it avoids circles passing through any four generator points. This assumption is needed because in some algorithms, subsets of the generator set may be used to construct the entire diagram. Numerical Precision Assumption In our initial description of algorithms we assume precision arithmetic is available for all calculations. Given the tendency towards numerical instability, it is conventional [65] to consider geometric algorithms in a theoretically ideal world in which precision arithmetic is the norm. The algorithms presented in this chapter are for constructing Voronoi diagrams in the Euclidean plane. Formidable difficulties exist in higher-dimensional construction methods. The algorithms for producing them do not always follow directly from the two-dimensional case. We consider the details of this subject 37 outside the scope of this thesis, but do include a brief overview of higher- dimensional construction algorithms and their performance. Some authors prefer to first construct the Delaunay tessellation for a given set of generating points (e.g., see [38] and [95]). Algorithm performance for constructing a Delaunay tessellation is, on average and in the worst case, equivalent to constructing the corresponding Voronoi diagram. As indicated in Chapter 3, either structure can be obtained in 0(n) time from the corresponding dual (i.e., they are linear-time transformable). There is no apparent savings in performance (and in fact, storage space), but there are however, some implementation advantages to first producing a Delaunay tessellation. The primary benefit is that the regions of a non-degenerate Delaunay tessellation are simplices, whereas Voronoi regions are generally multi-faceted polytopes. Geometric structures that have a constant number of edges (or facets) surrounding each region are generally more convenient to work with (e.g., three edges per Delaunay triangle in R2, four faces per Delaunay tetrahedron in R3, etc.). Therefore, we conclude this chapter with a brief discussion of algorithms for constructing the Delaunay tessellation. 4.1 The Polygon Construction Algorithm Before introducing the first algorithm, we state an alternative definition of the Voronoi diagram that forms the basis of this construction. The justification for this definition is the well-known fact that any convex polygon in the plane is the intersection of the half planes generated by the edges of the polygon. Applying this notion, we construct each Voronoi polygon V(pi) individually, thereby producing the entire 'V(P) structure. Voronoi Diagram (alternative definition) Let P= {pj ,p2, ...,/?}c R2, 2 < n j*i, ij e I be a set of Voronoi generators. Then region denoted by (4-!> is the Voronoi polygon associated with p,. Furthermore, the set V(i>) = {V(p,), V(p2),.... V(p)} (4.2) is the Voronoi diagram generated by set P. 38 We can now outline a construction method that builds the Voronoi diagram from its Voronoi polygons one by one. Algorithm 4.1: The polygon construction method Input: Set of n Voronoi generatorsp, ,p2, p - Output: Voronoi diagram V(P) = {V(p,), V(p2),..., V(p)}. Procedure: Step 1. For each / = 1,2,..., n, repeat 1.1. Generate n 1 half planes H(pt ,p}), \ < j < n, j i. 1.2. From the set of H(pt, Pj), construct the common intersection for each Voronoi polygon V(pt). Step 2. Return the set {V(p,), V(p2),..., V(p)} as output and terminate. This method is works because it is essentially a restatement of the above definition, but we cannot call this a practical algorithm. As we now show, the performance of the polygon construction method is inefficient. Constructing the half planes for two points p, and p, (all of Step 1) requires only constant time. Therefore, the time required to construct n 1 half planes is n 1 times a scalar, say x. To produce the intersection of the half planes, first construct the intersection of two half planes, giving a polygon with two edges (although an unbounded polygon). Next, construct the intersection of this polygon with the third half plane, and so on. At the kth iteration of successively adding half planes we must determine the intersection of a &-sided polygon with the next half plane that is introduced. In the worst case, this step (Step 1.2) requires time proportional to k to determine whether any of the k sides of the polygon intersects the line of the half plane. Therefore, producing the intersection of n 1 half planes requires a time proportional to 1 + 2 +...+ (n 2). The sum of this series is obviously (n 2)(n l)/2. Introducing another constant y, we have y(n 2)(n 1) for the time required to construct n 1 half planes. After repeating this procedure for all V(p), the total time required by this algorithm is T(n) = n{x(n 1) +y(n 2){n 1)} = 0(n3). (4.3) Several optimization techniques have been applied to the polygon construction algorithm which slightly improves its performance. The construction of the intersection of n 1 half planes can be brought down to 0(n log n) time [72]. This decreases the overall performance of the algorithm to 0(n2 log n) which is still impractical for most applications. 39 The output of this method is a just a set polygons which does not contain explicit topological information about the structure as a whole. If we need to extract various kinds of information from a Voronoi diagram, a more detailed construction method should be chosen. Although the polygon construction method is not recommended for general purposes, it can be useful for special types of simulation where the goal is to generate many samples of ^V(P) to obtain statistical data (e.g., the average number of edges per polygon or the average number of edges per vertex). 4.2 The Incremental Method The most important procedure for the construction of Voronoi diagrams is the incremental method. Conceptually simple, it has become one of the most often- used algorithms in this area of computational geometry. It begins by first constructing a Voronoi diagram from a minimal set of two or three generator points. The diagram evolves by successive introducing additional generator points into the existing diagram and constructing the polygons associated with these newly added points. Formally, we have the following procedure. LetVm m = 1,2,..., n represent a Voronoi diagram built from the first m generator pointsp7 ,p2, ...,pm The incremental method proceeds by transforming V,_j to Vi for every i. When i = mwe have constructed the complete Voronoi diagram on m points. Figure 4.1 shows an evolving Voronoi diagram ^V6 constructed from six generator points pj through p6. Also shown is a new polygon being formed by the introduction of an additional generator p7. The new polygon is constructed by first locating the generator point of the region that p7 is in, which in this example is p,. The perpendicular bisector of these two generators intersects two Voronoi edges of the Voronoi polygon V(p,). Label these two intersection points x7 and x2 such that generator p7 is to the left of line xt x2 when one considers the line directional fromx; to x2. The line segment x} x2 will partition V(p,) into two polygons associated with two generators p; and p7 respectively. This newly constructed line segment is a Voronoi edge common to polygons p, and p7. 40 Figure 4.1 A Voronoi diagram on six points with a new generator p7 being introduced. We continue by adding edges to the boundary of polygon p7. The perpendicular bisector of p2 and p7 will intersect an edge of Voronoi polygon V(p2) at a point, say x5. Similarly, we proceed to construct line segments in a counterclockwise manner around generator p7. As above, these line segments are constructed as perpendicular bisectors between p7 and its neighboring generators. When the starting point x, is reached, the polygon surrounding p7 is the Voronoi polygon V(p7) generated by p7. The edges constructed along the way are the Voronoi edges of V(p7). As a final cleanup step, the substructure contained within polygon V(p7) is deleted, thereby producing the Voronoi diagram V7. This description of the incremental method is a basic outline of the steps required to construct a Voronoi diagram. To produce a correct algorithm an additional procedure must be included to avoid geometric inaccuracies in the resulting diagram. A common problem with the incremental method is that the procedure used to create the Voronoi polygon boundaries only performs well when the polygon of the added generator is bounded by existing Voronoi polygons. The edges of a unbounded polygon under construction may not necessarily intersect an 41 existing polygons edge. This normally requires an additional subroutine, although the extra processing can be avoided by using the following simple technique. Let V(P) be a Voronoi diagram with generator set P = {p3 ,p2 , Begin by renumbering the generators as p4 ,p5,..., p. This introduces three more generators to the original set of n points (i.e., n is now three greater than it was). The coordinates of the three additional generators p, ,p2, and p3 are chosen to form a triangle large enough to contain all n of the original generators of "V(P). Construct the simple Voronoi diagram for p, ,p2, and p3. Add to this three- generator Voronoi diagram the original set of n generators, one by one, according to the incremental procedure described above. We see that the convex hull of Pi ,p2, , Pm is always the enclosing triangle p, p2 p3 for any m = 4, 5,..., n. Now recall from Chapter 2 that a Voronoi polygon V(p,) is infinite if and only if p, is on the boundary of CH(P). Thus, the points on the convex hull (i.e., the triangle Pi P2 Pi) are precisely the generators that form the infinite polygons of 'V(P). By using this procedure it is assured that a Voronoi diagram 'Vm, under incremental construction, is always finite. Following the construction of this augmented Voronoi diagram, it may be necessary to eliminate the additional three points and retrieve the diagram generated by the original set P. This can be done if the triangle formed by the additional three generators is sufficiently large. Figure 4.2(a) shows the central portion of the Voronoi diagram in 4.2(b) magnified to where the additional three generators have no effect. To achieve this, first assume that the generators p4,p5.....pr. are all contained within the unit square S {(x, ,y) | 0 ^ x, y < 1}. We can now select the additional three generators to be: Pl = (0.5, 3^ + 0.5) p = (-3^ + 0.5, -3^? + 0.5) 2 4 4 p = (3-^ + 0.5, -3 + 0.5) (4.4) 4 4 These points form a triangle that is large enough to entirely enclose the unit square defined above and along with it, the generators p4,p5,..., p. Other choices for Pi ,p2, and p3 can be made providing that triangle p} p2p3 entirely contains S. It can be seen from Figure 4.2(b) that none of the edges of the Voronoi polygons 42 associated with the additional three generators intersect the boundary of the unit square S. Therefore, once the augmented Voronoi diagram is constructed in this manner we can delete the outer three polygons without affecting the diagram we originally intended to create. Pi A Figure 4.2 Avoiding infinite polygons: (a) Voronoi diagram for original generators; (b) same Voronoi diagram with three additional generators. 43 For the incremental method to be efficient and achieve a lower bound performance of 0(n log n), another subroutine is required in the construction algorithm. We provide a general description of this optimizing procedure, the details of which, may be found in the literature (see references below). Recall that to construct the edges of a polygon associated with a new generator we must first locate the closest existing generator. In general, the time complexity of this task is 0(rr) as follows. Let Pi and pm be the old and new generators respectively. For p, to be the closest generator to pm it must satisfy d{Pi, pj Z d(pj, pj for all; = 1, 2,..., m 1. (4.5) If the distance from pm to every other generator is calculated, then p, can be found. However, this is not efficient as it requires 0(1 + 2 +...+ (n 1)) = 0(n2) time. This would increase the average time complexity of the incremental method to 0(n2). To maintain linear performance on average, a nearest neighbor search is used to find the location of the closest generator. Algorithm 4.2: Nearest neighbor search Input: m generator points p, ,p2,..., pm Voronoi diagram "Vm_,, and an initial guess p,, 1 < / < m 1. Output: Generator nearest to pm . Procedure: Step 1. Among the generators whose Voronoi polygons are adjacent to V(p,), locate the generator pj with minimum distance to pt, taken over all generators pk. The distance function is defined as d(pj, pm) = mint d(pk, pm). (4.6) Step 2. If d(pi, pm)< dipj, pm) then return p,, else Pj pj and go to Step 1. It helps to make a good initial guess at pt as this will greatly influence the algorithms performance. To facilitate this, incremental algorithms sometimes utilize a bucketing technique. The basis of the bucketing technique is to first partition the space containing the generators into a regular grid. The bucketing is implemented with a data structure known as a quaternary tree which uses the grid to maintain an approximately uniform distribution of these sites throughout the polygon insertion process. 44 The bucketing technique can conveniently support a dual purpose here. Recall that during incremental construction of a Voronoi diagram, the substructure contained within each newly added Voronoi polygon is deleted from the diagram. If the substructure is large, its components cannot be located and removed in linear time. Property VG5 from Chapter 2 states that the average number of Voronoi edges per Voronoi polygon is not greater than six. However, this does not necessarily indicate that the average number of edges for a newly added Voronoi polygon will be six or less. To address this performance concern in an implementation, the bucketing technique described above can be simultaneously used for this task as well. The details of quaternary tree implementation are outside the scope of this thesis. The reader is referred to [63] and [64] for a description of quaternary trees and [65] for implementation details on applying the bucketing technique to the incremental construction method. In summary, we provide the incremental construction algorithm It is assumed that the interested reader will research the specific implementation details. Algorithm 4.3: Incremental method for constructing 'V(F) Input: A set {p4,p5,..., p) of n 3 generator points contained within the unit square S defined above. Output: Voronoi diagram VfP), P = {p,, p2,..., p} withp}, p2 andp3 included as the additional generators defined in equation 4.4. Procedure: Step 1. Find positive integer k such that 4* is closest to n, divide S into 4* square buckets (i.e., regular grid), and construct a quaternary tree to facilitate the bucketing technique (see references below). Step 2. Within the quaternary tree data structure, reorder the generators to bep4, p5,..., p. Step 3. Construct Voronoi diagram using the additional three generators p,, p2 and p3 defined in equation 4.4. At this step the winged-edge data structure used to represent the Voronoi diagram is created. Step 4. For each m = 4, 5,..., n repeat: 4.1. Using Algorithm 4.2 with the initial guess provided by the quaternary tree data structure, locate the generator pt nearest to pm. 45 4.2. 4.3. 4.4. Step 5. Locate points x, and x2 that are the intersections of the perpendicular bisector ofpt and pm with the boundary of m- Using the procedure described above, extend the sequence of line segments (x,x2, x2x3,..., to form the edges of the Voronoi polygon generated by pm . Delete from the substructure enclosed within the sequence of Voronoi edges of pm Name the resulting diagram 'Vm and update the representation of the Voronoi diagram in the winged-edge data structure. Return 'V='Vn. To determine the overall performance of the incremental method, the individual steps are examined. Steps 1 and 2 can be executed in 0(n) time. Steps 3 is performed simply in 0(1) time. Step 4 (4.1,..., 4.4) can also be achieved in 0(1) time on average for each m. Likewise, Step 5 is only 0(1). Hence, Algorithm 4.3 achieves an overall performance of 0(ri) on average. It should be noted that this average linear performance is met if the generator point set is distributed randomly in the unit square, which is often the case in real-world applications. Also note that the efficiency of performance is largely due to the use of the bucketing technique which brings the average time complexity of the incremental method from 0(2) down to 0(ri). For an early paper on applying bucketing and quaternary trees to the construction of Voronoi diagrams, see [64]. In this paper, Ohya, Iri, and Murota report linear performance for Voronoi diagram construction over a variety of generator point distributions. 46 4.3 The Divide-and-Conquer Method The divide-and-conquer method is one of the most fundamental ways of solving complex problems. It works by recursively dividing a problem into two or more subproblems which are easier (or faster) to solve individually. The solution to the overall problem is accomplished by merging the solutions of the subproblems as the algorithm backs out of the recursive steps. The Voronoi diagram, despite its apparent complexity, is very well-suited to construction by divide-and-conquer. The method described here exploits various structural properties of the Voronoi diagram which facilitates the merging of the subproblems in linear time. Prior to setting up the divide-and-conquer algorithm we first sort the set of generator points in increasing order by their x coordinates. The sorting can be achieved in 0{n log n) time by employing any optimal sorting routine (i.e., merge sort, heap sort, quick sort, etc.). As usual, the set of generators is represented by P-{pi,p2, , Pn}- As we will presently see, this ordering is required to find the median x-coordinate for each recursive step. This median may be interpreted as a vertical line which divides the Voronoi diagram into two distinct sets of generators. Recall that two assumptions were initially imposed to avoid numerical and algorithmic instability. For the divide-and-conquer paradigm we introduce a third assumption which uniquely defines the ordering of P. Generator Alignment Assumption Among the set of generators P = {pt, p2 no two generators align vertically. That is, given two generators p, = (xa, yb) and pj = (xc, yd), a* c for all Pi, pj e P. We now formally state the divide-and-conquer algorithm. Algorithm 4.4: Divide-and-conquer method for constructing 'V(P) Input: The ordered set P = {p, ,p2, ...,pn} of Voronoi generators and the number n of points in the set P. Output: The Voronoi diagram 'V(P) generated from set P. Procedure: Step 1. If n < 3, then construct 'V(P) directly and go to Step 3. Step 2. Else do 2.1. Let t be the integral part of nil. Divide P into Ph={pi,...,p,} andPR = {pt+I,...,p} wherePL andPR 47 2.2. 2.3. 2.4. Step 3. denote the left and right generators of the set P. Construct the left Voronoi diagram "VL for generators PL using Algorithm 4.4 (i.e., use this algorithm recursively). Construct the right Voronoi diagram "VR for generators PR using Algorithm 4.4. Merge VL and 'VK into the Voronoi diagram V generated by point set P. For this step, use Algorithm 4.6 which is stated below. Return 'V. The crucial step in this algorithm is the merging of the two Voronoi diagrams. This is done by constructing a merge chain, that is, the set of new Voronoi edges and Voronoi vertices that result by the merging of "VL and "VR. The success of this procedure depends on the topological properties of the merge chain and subsequently, how rapidly the merge chain between the two diagrams can be constructed. The other steps pose no difficulties in implementation, nor do they inflict excessive performance requirements on the algorithm as a whole. On the other hand, the merge chain construction can become a performance bottleneck on the entire algorithm if it is not implemented properly. We proceed by describing how this can be accomplished efficiently. Figure 4.3 shows two Voronoi diagrams "VL and "VR constructed from the left and right generators, PL and PR respectively. The left generators PL are shown as filled circles, whereas the right generators PR are shown as unfilled circles. Likewise, the solid lines represent 'VL and the dashed lines represent "VR. The merge chain (denoted by M) is shown as a heavy solid line. To merge the two diagrams we must generate new Voronoi edges and Voronoi vertices that are common to both diagrams. These new edges and new vertices (i.e., merge chain M) originate from the interaction between the Voronoi polygons of the generators of "VL and VR. Notice that the individual segments of Mare always perpendicular bisectors between a left and right generator as they should be in the composite Voronoi diagram. A remaining last step is to delete the superfluous line segments that intersect the merge chain. 48 Note that every horizontal line in the plane intersects the merge chain in exactly one point. This is a result of the merge chain being vertically monotonic. To illustrate this, Figure 4.4 shows an impossible merge chain. If a horizontal line existed that intersected M in two or more points, some generator, say st e PL , would have a larger x-coordinate than some generator, say t, e PR. This would contradict the assumption that PL and PR are entirely separated by a median vertical line. 49 The merge routine proceeds by moving a point upward from the bottom of the entire set of generators that form the two Voronoi diagrams being merged. This point, say z, is always incident with the merge chain and traces a polygonal line that runs between the left generators and right generators. At any time during the merging process point z is always equidistant to two generators, one in PL and one in PK. The path of point z can be thought of as the locus of points equidistant to the left and right sets of generators. To produce the combined diagram, this polygonal line must be constructed. The final step is to delete the extraneous line segments of "VL and VR that lie to the right and left, respectively, of the polygonal line. Before the polygonal line can be constructed, an additional subroutine is required. Since the polygonal line is built from the bottom upward, the bottommost edge needs to be determined. Likewise, the merging process terminates when the topmost edge is encountered, so it too must be found. To find them we rely on an algorithmic technique for locating a line that joins two distinct convex polygons with the condition that all other points of both polygons are on the same side of the line. Let Â£/be a convex polygon. One way to define U is by a cyclic set of vertices that surround the polygon in a counterclockwise manner. We denote this cyclic list of s 50 elements as {u,, u2,..., us, Uj}. For a vertex w, e U, let cnext[w;] be the immediate predecessor (i.e., clockwise next) of w, and ccnext[w,] be the immediate successor (i.e., counterclockwise next) of w,. If UL is a convex polygon to the left of UR, then the x coordinates of w, e UL are smaller than the x coordinates of w, e UR, for all i = 1, 2,s. Given vertex v e UL and vertex we UR, the directed line segment CS(v, w) is called a common support of UL and UR if all the vertices of UL and UR are on the same side of CS(v, w) or are on this line. When the vertices are above CS(v, w) we have the lower common support, denoted LCS(v, w); when below CS(v, w), we have the upper common support, denoted UCS(v, w). Figure 4.5 shows two convex polygons, UL and UR, along with the lower and upper common supports. Upper Common Support Figure 4.5 Lower and upper common supports of two convex polygons The following algorithm locates the lower common support for two convex polygons. The upper common support can be found with a slight modification to the algorithm. 51 Algorithm 4.5: Find the lower common support Input: Convex polygons UL and UR such that the maximum x coordinate for all vertices in UL is smaller than the minimum* coordinate for all vertices in UR . Output: A pair of coordinates consisting of vertex v in UL and vertex w in UR such that CS(v, w) forms the lower common support LCS(v, w) of UL and UR . Procedure: Step 1. Step 2. 2.1. 2.2. Step 3. Find the vertex v in Uh with the largest x coordinate, and the vertex w in UR with the smallest x coordinate. Perform 2.1 and 2.2 alternately until v and w do not change. Repeat v cnext[v] until vertex cnext[v] is lower than LCS(v, w). Repeat w ccnext[w] until vertex ccnext[w] is lower than LCS(v, w). Return LCS(v, w). To obtain the performance of this algorithm, let n represent the total number of vertices in Ut and UR combined. Step 1 requires linear time as the initial pair of vertices, v and w, can be located by scanning the two lists of generator points. Since the two cyclic lists have a combined total of n elements, the complexity for Step 2 is also 0(n). Therefore, the performance of Algorithm 4.5 is 0(ri). A 0(log n) algorithm for finding common supports is presented in [69], but the savings does not contribute to the performance of the overall divide-and-conquer algorithm which runs at 0(n log ri). With this subprocedure in hand, we can now describe the merging process. The objective is to construct the polygonal line of new Voronoi edges and vertices. Let LCS(ph ,pp) denote the lower common support for CH(PL) and CH(PR). The perpendicular bisector of line pLpR, represented by B(pL ,pR), will be the lowermost segment of the polygonal line. From Figure 4.3 we can see that pL pR is not an edge of either PL or but it is an edge of PL uPR Therefore, the merge chain begins with the pair of points (pL pR). In the description that follows, it is helpful to refer to Figure 4.3. Let z be a point that travels upward along the perpendicular bisector between PL and PR, and let z_ be a point at infinity that is incident to this bisector. Also let z; denote the first point at which z initially intersects a Voronoi edge e, of either Voronoi diagram or "VR (see Figure 4.3). This segment of the bisector z_jz} that 52 z has thus far traveled becomes the first Voronoi edge of merge chain. Likewise, point zt is the first Voronoi vertex on the merge chain. As it continues to move, point z then alters a Voronoi region at point z,. If e is an edge of X it modifies that Voronoi diagram. Likewise, if e is an edge of X, it modifies that Voronoi diagram If e is an edge of X, thenpL is replaced by the left generator that is on the opposite side of e. Similarly, if e is an edge of X, thenpRis replaced by the other right generator that produced edge e. Notice that z, is incident to the perpendicular bisector between the updated pair of generators pL and pR Let z resume its upward movement from zh this time along the new bisector. It will cross another Voronoi edge of either X or X Denote this intersection point z2. This provides the next new Voronoi edge and the next new Voronoi vertex. In this manner, z will continue through all vertices z, until reaching the pair pL, pR that defines the upper common support. From here point z heads off as the perpendicular bisector to line pLpR As in the case of the first (and bottommost) edge, this final (and uppermost) Voronoi edge is also a semi-infinite ray that can be thought of as incident with point z. With the foregoing explanation, we can now state the merge procedure. Algorithm 4.6: Merging two Voronoi diagrams Input: Two Voronoi diagrams X and X constructed from generator sets PL and PR, respectively, with the condition that the generators of VL have smaller x coordinates than those of X . Output: Voronoi diagram V constructed from the union of the two sets of generators, that is, PL u PR . Procedure: Step 1. Step 2. Step 3. Step 4. 4.1. 4.2. 4.3. 4.4. 4.5. Construct the convex hulls of PL and PR . Using Algorithm 4.5, locate the lower common support LCSiy, w) and the upper common support UCS(v, w). w0 ~ the point on B(pL pp) at infinity in the negative y direction, and / 0. Repeat Use Step 4.5 below (i.e., 4.1 and 4.5 are equivalent). i i + 1. Locate the point aL (other than w,w) at the intersection of B(Pi, Pr) and the boundary of V(pL). Locate the point % (other than w,_,) at the intersection of B(Pl 5 Pr) and the boundary of V(pR). If aL has a smaller y coordinate than aR, 53 Then w,<- aL and PL - the generator on the other side of the Voronoi edge containing at, Else w,- aR, and PR *- the generator on the other side of the Voronoi edge containing aR . Until UCS(pL, pR) is reached. Step 5. m - /. wm+, ~ the point on B(pL pp) at infinity in the positive y direction. Step 6. Add the polygonal line (w0w,, Wj\v2,..., wmwm+I). Step 7. Delete from "VL the edge segments to the right of the polygonal line and delete from VR the edge segments to the left of the polygonal line. Step 8. Return the (merged) Voronoi diagram "V. As usual, to assess performance we determine the complexity of each step. Step 1 can be done in linear 0(ri) time since we can construct the Voronoi diagrams "VL and "VR as a direct retrieval from the winged-edge data structure (recall Algorithm 3.2). Step 2 also requires 0(n) by using Algorithm 4.5. Step 3 is done in constant time. Step 4 is performed in at most 0(n) time because there are at most n Voronoi edges added as a result of the polygonal line. Steps 4.1 and 4.5 only need constant time. Steps 4.2 and 4.3 have a time complexity proportional to the number of edges on the individual boundaries of "VL and 'VR. However, these Voronoi diagrams are not bounded by a constant. To decrease the total number of edges considered in Steps 4.2 and 4.3 for all repetitions of Step 4, we can use an optimizing technique. We only examine the generators on the side of the edge (being constructed) that is opposite the direction in which the edge has turned at the previous vertex. If the edge turns to the left at the vertex, then examine the edges to the left of the given edge. Similarly, if the edge turns to right at the vertex, then examine the generators to the right of the given edge. We thereby avoid repeated checks of the irrelevant edges throughout the merge process. Using this simple technique, the overall performance of Step 4 is 0{ri). Step 5 is done in constant time. Steps 6, 7, and 8 each require 0(n) time. Hence, the overall performance of Algorithm 4.6 is 0(n). This is not bad for such a seemingly complex procedure. 54 We can now determine the performance of the main divide-and-conquer algorithm (i.e., Algorithm 4.4). Steps 1,2.1, and 3 are done in constant time. As we have just shown, Step 2.4 has a performance of O(ri). In Steps 2.2 and 2.3 the main algorithm invokes itself recursively, and as such, are the limiting factors on performance. For a set of n generators, let T{n) represent the time required by the Algorithm 4.4. Since the algorithm calls itself twice with half of the original input as the new input, we have the recurrence T(n) = 2T{n/2) + 0(n) (4.7) where 0(n) in this equation represents the performance of all steps other than the recursive steps. Considering that 71(1) and T(2) are constant, equation 4.7 implies that 71(n) = 0(n log n). As noted in Chapter 3.3, the best performance any Voronoi diagram construction algorithm can do is 0(n log n). Therefore, the divide-and-conquer method is optimal insofar as worst-case performance is concerned. The divide-and-conquer approach to Voronoi diagram construction has a rich history of study. Shamos and Hoey [82] first attacked the problem using this technique and a number of variations on the original theme have since been presented. It was determined in [63] that the divide-and-conquer method is not optimal on average as it requires 0(n log n) time. However, the average performance can be improved if the generators are uniformly distributed in the plane. Given this condition, Dwyer in [26] reported an average performance of 0(n log log n) while Katajainen and Koppinen [48] claim an average 0(n) variation on the divide-and-conquer approach. 55 4.4 The Plane Sweep Method In a 1986 paper by Fortune [34], a fundamental and elegant technique of computational geometry known as the plane sweep method was applied to the construction of Voronoi diagrams. This method uses a vertical sweepline, which passes (i.e., sweeps) across the plane from left to right. This technique is used to scan the plane in order to locate geometric objects within it. As the sweepline moves, it encounters objects such as points and edges. At every encounter, information is gained which can be used to solve a portion of the problem at hand. In this way, the plane sweep method can be used to construct geometric objects such as Voronoi diagrams and Delaunay triangulations. This technique is conceptually simple, but must be implemented carefully to produce correct results. As in the case of other Voronoi diagram construction methods, the plane sweep technique requires a bit of up-front processing. This preprocessing involves a planar transformation (i.e., mapping) of the set of Voronoi generator points to a corresponding isomorphic Voronoi diagram. Without this initial transformation, the plane sweep method cannot be directly applied to Voronoi diagram construction as follows. As the sweepline moves across the plane, along it, the Voronoi vertices and Voronoi edges are constructed. However, some vertices and edges surrounding a given polygon (and the polygon itself) are encountered by the sweepline before it has reached the generator of that polygon. What is needed is a way to predict the locations of the generators in the half plane that has not yet been swept. To overcome this, we develop a mapping which places each generator Pi at the leftmost point of the Voronoi polygon F(p,). This enables us to use the plane sweep method because the generator of a Voronoi polygon will now be the first object of that polygon to be encountered during the plane sweep. As usual, let P= {p, ,p2, ,p) be a set of n generators for Voronoi diagram 'V(P). In addition to the three assumptions stated previously in this chapter, we now add a fourth assumption which simplifies the description of the plane sweep algorithm. Horizontal Alignment Assumption Any generator point and any vertex on the boundary of the associated polygon do not align horizontally. The transformation function (i.e., mapping) is now defined as follows. Given a point p, define rip) min {dip, p{) | 1 the largest empty circle centered at p. Now let x(p) and yip) represent the x 56 coordinate and y coordinate, respectively, of point p. For a point p, a mapping function (p is defined from M2 to M2 as If/? Â£ V(pj), then
A point p is unchanged in this mapping if/? is a generator. That is, (pip) =/? if and respectively.
segment of a hyperbola that opens to the right and has p{ as its leftmost point. Each
to a hyperbolic segment.
diagram. Determining the topological structure of V(P) is equivalent to
point generator /?,. This implies that the images of all Voronoi vertices and The following two figures shows the mapping function (p applied to a Voronoi diagram generated from five points. Figure 4.7 shows the original diagram "V(P), whereas Figure 4.8 shows the transformed diagram cp(V(P)). Figure 4.7 Voronoi diagram 'V(F) constructed from five generators. 58 Figure 4.8 The image under map
Voronoi diagram in Figure 4.7.
point of any Voronoi polygon in (p(V(i>)). Therefore, cp(^) has exactly one
one of these events. Two types of events may occur, as can be seen in Figure 4.8.
(p(l/(P)) can divided into two segments at the leftmost point. Denote the upper
Procedure:
(q>(V(/0), Kipj.pi,
the intersection of Kipj, pj) and K(pj, pj).
ffiPi, pj), h(Pj, Pk), and h.
the plane sweep algorithm at work. For insight into the correctness of the
proved by Fortune in [34] and [35].
(2) At each Voronoi vertex in
right and the other two edges extend to the left.
The performance of Algorithm 4.7 is now determined as follows. Step 1 requires |