Citation
Voronoi diagrams

Material Information

Title:
Voronoi diagrams properties, algorithms, and applications
Creator:
Romano, Don
Publication Date:
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 )

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 represents a Delaunay tessellation with n triangular regions Tt.
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 ne < 3 ng 6 , (2.13)
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
only if p £ P. In the discussion and figures that follow, we represent the Voronoi
vertex q, the Voronoi edge e, the Voronoi polygon V(pt), and the Voronoi diagram
Vas (piq), (pie), respectively.
Let Pi and pj denote two generators such that xipi) > xipj) as shown in Figure 4.6.
Also let e be the Voronoi edge separating p, and pj. Function

segment of a hyperbola that opens to the right and has p{ as its leftmost point. Each
hyperbolic segment corresponds to a Voronoi edge. Consequently, the image
under cp, of each Voronoi polygon V(pi), is the intersection of hyperbolically
bounded half planes of M2.
Figure 4.6 Under function to a hyperbolic segment.
57


We now show how the transformation

diagram. Determining the topological structure of V(P) is equivalent to
determining that of cp(V(P)). It should be observed that the for all generators p,
other than the leftmost one, the transformed polygon point generator /?,. This implies that the images of all Voronoi vertices and
Voronoi edges associated with p, fall to the right of pt when mapped under
function (V(P)).
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.
To further describe Figure 4.8, take note that the image of a Voronoi vertex is
referred to as a Voronoi vertex of (p("V(P)) and the image of a Voronoi edge is
referred to as a Voronoi edge of (p(V(P). Also note that by the horizontal
alignment assumption, the image point of any Voronoi polygon in (p(V(i>)). Therefore, cp(^) has exactly one
Voronoi edge extending to the right and two Voronoi edges extending to the left.
As the vertical sweepline moves from left to right, generators and vertices are
encountered. The structure of the Voronoi diagram one of these events. Two types of events may occur, as can be seen in Figure 4.8.
One type happens when the sweepline hits a generator as at lines Sj ,s2,s3, and s4.
The other type of event occurs when the sweepline hits a Voronoi vertex, as at line
s5 (i.e., the sweepline hits cp(g4)).
59


As we have seen, the perpendicular bisector (i.e., Voronoi edge) of two generator
points Pi and pj is mapped by function

(p(l/(P)) can divided into two segments at the leftmost point. Denote the upper
segment by Kip,, pj) and the lower segment by Kip,, pj). Call these the upper and
lower segments half hyperbolas.
Prior to stating the plane sweep algorithm, a couple of data structure preliminaries
are required. First, we need to represent the structure of (p( V(P)) as it appears
along the sweepline. To facilitate this, a list L is used to record the regions and
boundary edges which alternate in occurrence along the sweepline. Next, a set of
points Q in the plane is required to represent the potential points at which future
events may happen along the sweepline. We initialize Q with the set of generators.
Voronoi vertices which are candidates for future events are added to Q as they are
located by the sweepline. The following algorithm was first presented by Fortune
in [34]. The treatment here is due to Okabe et al. [65].
Algorithm 4.7: Plane sweep method for constructing "V(P)
Input: Set P = {p,, p2,..., p) of n generators.
Output: Voronoi diagram Procedure:
Step 1.
Step 2.
Step 3.
Step 4.
4.1.
4.2.
4.3.
^(Pj, pk)) on L by h = Kip,, pj) or h = h+ipt, pj.
Q-P.
Select and delete the leftmost point, say pt, from Q.
L - the list consisting of a single region of (p(VfP)).
Repeat (4.1,4.2, and 4.3)
Select and delete the leftmost point w from Q.
If w is a generator, say w =p,,
then (do 4.2.1,..., 4.2.3)
4.2.1. Locate the region cp( V(pjj) in L containing p,.
4.2.2. Replace (q>(V(/0), Kipj.pi, 4.2.3. Add to Q the intersection of h (pj, pi) with the
immediate lower half hyperbola in L, and the
intersection Kipj, pj) with the immediate upper half
hyperbola in L.
If w is an intersection, say w = cp(qj,
then (do 4.3.1,..., 4.3.4) (it is assumed here that the intersection of Kipj, pj) and K(pj, pj).
4.3.1. Replace the subsequence iKipi, pj), (p( Vipjj),
60


4.3.2. Delete from Q any intersections of h{pi, pj) or
ff(pj, pj) with others.
4.3.3. Add to Q any intersections of h with its immediate
upper half hyperbola and its immediate lower half
hyperbola in L.
4.3.4. Mark ffiPi, pj), h(Pj, Pk), and h.
Until Q is empty.
Step 5. Report all half hyperbolas ever listed in L, all the Voronoi
vertices marked in Step 4.3.4, and the incidence relations
among them.
To understand how Algorithm 4.7 works, it is helpful to step through it while
referring to a diagram of the resulting the plane sweep algorithm at work. For insight into the correctness of the
algorithm we mention several characteristics of proved by Fortune in [34] and [35].
(1) The generator pt, other than the leftmost generator, lies at the
leftmost point of the Voronoi region (2) At each Voronoi vertex in right and the other two edges extend to the left.
(3) For every Voronoi vertex in cp("V(P)), the two edges extending to
the left will be adjacent in L at some time during the plane sweep.
This insures that the points inserted in Steps 4.2.3 and 4.3.3 include
all the Voronoi vertices.
Note that this algorithm constructs the image of Voronoi diagram 'V(P) under the
mapping function information as the normal (i.e., not transformed) Voronoi diagram 'V(P). All the
relationships and properties that we expect from a normal Voronoi diagram can be
found in the transformed structure. If the normal Voronoi diagram is desired,
perhaps for a visual representation, it can be obtained by applying the inverse
function of The performance of Algorithm 4.7 is now determined as follows. Step 1 requires
assigning the generator set P to Q in sorted order. This clearly requires 0(n log n)
61


time. The basic operations on O can be facilitated through the use of a priority
queue which can be implemented with a heap data structure. Using this data
structure, the operations on Q (with n points) can be performed in 0(log n) time.
So Step 2 takes <9(log n). Step 3, a simple assignment, is constant time.
Since the number of Voronoi edges and Voronoi vertices are both of 0(n), the
total number of elements that can ever be listed in L is of 0(n). Similarly, since
there are only 0(n) pairs of Voronoi edges that can be adjacent to each other in L,
the total number of elements ever placed in Q is of 0{ri). This implies that there
can be at most 0(n) candidates for Voronoi vertices. As a result, Steps 4.1,4.2,
and 4.3 each have a performance time of 0(log n), but notice how they are
repeated 0{ri) times in Step 4. Therefore, the total time complexity of Algorithm
4.7 is 0(n log n).
4.5 Comparative Performance of Construction Methods
The following table summarizes the performance results of the four Voronoi
diagram construction algorithms given above. As previously noted, certain
optimization techniques can sometimes be applied to improve algorithmic
performance. Table 4.1 indicates the worst-case and average time complexity of
the basic methods without regard to additional optimization procedures.
Table 4.1 Algorithm Performance
Polygon Construction Incremental Method Divide-and- Conquer Plane Sweep Method
Worst-Case 0(n3) Average 0(n)* 0{n) 0(n log n) 0(n log n)
* The average linear performance for the polygon construction method is attained
through significant optimization of the algorithm.
62


4.6 Numerical Accuracy and Implementation Issues
The numerical limitations of computers, combined with degenerate or nearly
degenerate input data, can cause considerable topological inconsistencies to result
in Voronoi diagrams constructed by the foregoing algorithms. These pitfalls can
undermine the seemingly best of algorithms. This is not only a problem for
Voronoi diagrams, but also for geometric algorithms in general.
At the beginning of this chapter two assumptions were imposed to avoid structural
degeneracy and numerical errors during algorithm execution. Later, two additional
assumptions were introduced to add convenience to the construction process and
further avert inaccuracies. However, algorithms composed in such a manner may
not necessarily produce a valid or accurate Voronoi diagram.
Figure 4.9 illustrates a common topological inconsistency that is often encountered
in the incremental method. Recall that the incremental algorithm works by growing
the boundaries of newly added polygons to an existing structure. These boundaries
are, by definition, a sequence of perpendicular bisectors between adjacent
generator points. If the bisectors pass too close to a Voronoi vertex, the resulting
geometric relationships can be misread. This could result in the polygon boundary
not ending up at the starting point as it should. If this type of inaccuracy occurred,
the incremental algorithm may fail to produce a diagram of sufficient fidelity that
may be required by the specific application.
Figure 4.9 Topological inconsistency caused by
numerical errors.
63


Two types of corrective techniques can be used to avoid construction
inaccuracies: the topology-oriented approach and numerical precision measures.
By employing one or both of these techniques, the assumptions that were
previously imposed can be removed. Voronoi diagrams that are representative of
real-world applications can then be constructed with the assurance that they are
valid and exact. We touch briefly upon these two methods and provide the reader
with references from which to gather more detailed information.
When inferior precision arithmetic is all that is available, numerical errors are
inevitable. The reliance then, is upon combinatorial computation; a topology-
oriented approach that is commonly used when high numerical precision is
unavailable or impractical for the application at hand. This technique is based upon
avoiding degeneracy by adding redundant points to the initial generator set. This
augmented generator set does not alter the topological structure of the Voronoi
diagram, but if implemented correctly, can avert degeneracy. This is accomplished
by strategically placing the additional generators on existing Voronoi edges. If
these generators are properly placed, there will be exactly three generators incident
to a circle whose center is a Voronoi vertex (i.e., the condition necessary for all
non-degenerate Voronoi diagrams). Following construction of a valid diagram, the
additional points can be ignored (or removed) without impacting the integrity of
the overall structure.
The topology-oriented approach was first proposed in [86]. A publicly-available
program called VORONOI2, based upon this technique and written in FORTRAN,
was implemented by Sugihara and Iri [87]. An illuminating example of the
topology-oriented approach is also provided in [65].
The second, and perhaps more commonly used method, is to employ specific
numerical representations. One obvious approach is to use integers instead of reals
(i.e., floating point numbers). The input set of generators may then be represented
by a vector of integers. By doing so, extended-precision integer arithmetic can be
used, but at a considerable cost in overhead. The use of modular arithmetic over a
collection of finite fields has been proposed as an alternative [38], although little
information is available concerning actual results.
The natural alternative to extended-precision arithmetic is floating point arithmetic.
When used, floating point arithmetic avoids certain bit-complexity issues,
although non-degeneracy becomes much more stringent. This may be overcome by
the topology-oriented approach or as Fortune [38] suggests, by slightly perturbing
64


the exact positions of the Voronoi generators to allow for a floating-point
implementation.
Exact arithmetic seems to be the only known way to obtain a completely reliable
implementation, although there are disadvantages in inconvenience and the
potential cost in processor time. Multiple-word integer arithmetic has also been
considered for Voronoi diagram construction. With this type of implementation, a
considerable degradation in performance was reported by Edelsbrunner and Miicke
[28]. Karasik et al. [47] proposed the use of adaptive-precision rational arithmetic
to implement the divide-and-conquer algorithm in two dimensions. This
implementation proved to be considerably faster than rational arithmetic, but
somewhat slower than floating point arithmetic.
There has been significant research toward applying numerical analysis techniques
to geometric algorithms in general. See [36] and [46] for further investigations on
this subject. Numerical stability in Delaunay triangulation algorithms have also
been looked at in [37].
In summary, floating point arithmetic is the most common implementation choice.
Although not guaranteed to be completely reliable, it seems adequate enough for
most applications. Regardless of the arithmetic used in the implementation,
Fortune [38] recommends constructing the Delaunay triangulation rather than the
Voronoi diagram. The reason is to avoid manipulating computed values; Voronoi
vertices are sometimes misrepresented (i.e., poorly determined) when constructed
by certain implementations such as exact arithmetic.
4.7 Constructing the Voronoi Diagram in Higher Dimensions
Higher-dimensional Voronoi diagrams are important in a number of applications,
some of which are discussed in the next chapter. However, the construction of a
Voronoi diagram (and Delaunay tessellations) in three or higher dimensions is
significantly more complex than in the two-dimensional case. A contributing factor
is the increase in the size and complexity of the structure. It was determined by
Klee [50] that the number of &-faces for n generators in w-dimensional space is on
the order of for 0 < k < m. Clearly, this will not provide the 0(ri)
average performance that is enjoyed in two dimensions.
The commonly-used techniques of two-dimensional construction do not
necessarily extend to higher dimensions without substantial modification. For
example, the plane sweep method described in this Chapter 4.4 does not generalize
65


to three dimensions. Moreover, an implementation of the plane sweep in three
dimensions (i.e., space sweep) is impractical in that it would require at least
quadratic complexity [38]. For the algorithms that are practical, the
implementation details can be much more involved, post-checks are sometimes
required to insure correctness, and as one may expect, the run-time performance of
construction worsens.
We briefly describe a practical and popular method for higher-dimensional Voronoi
diagram construction known as the lift up transformation. This approach
generalizes nicely to n-dimensional space, is somewhat manageable to implement,
and has a tolerable time complexity if the diagrams are not overly large. The basis
for this method is best seen by first stating a theorem on Delaunay triangulations.
A derivation leading up to this theorem can be found in [67].
Theorem: The Delaunay triangulation of a set of points in two dimensions
is precisely the projection to the xy-plane of the lower convex hull of the
transformed points in three dimensions, transformed by mapping upwards
to the paraboloid z = x2 + y2.
This result says that a Delaunay triangulation can be obtained by mapping the set
of generators to a structure in three dimensions (i.e., a paraboloid), followed by a
projection back into two dimensions. After the Delaunay triangulation is produced,
we can easily derive the corresponding Voronoi diagram in linear time. Curiously,
this theorem generalizes to any dimension. For example, to construct a Delaunay
tessellation for a set of m-dimensional generators, first lift the point up to a
hyper-paraboloid in m+1 dimensions and then project the image back to a
hyperplane in m dimensions. The performance of this method in three dimensions
has a worst-case lower bound of 0{r^).
This technique was first presented by Brown [15] and later exploited by
Aurenhammer and Edelsbrunner [7]. A detailed explanation of this innovative
technique is given in [65]. Other treatments may be found in [6] and [38].
A number of other algorithms have been devised for constructing higher-
dimensional Voronoi diagrams including the popular incremental method. This
method is not to be confused with the two-dimensional incremental method
described in Chapter 4.2, although similarities exist between the two algorithms.
Sometimes called the beneath-beyond method, it begins by initially constructing
the convex hull of a small number of points. Additional points are added
66


incrementally, to give the convex hull and in turn the Delaunay tessellation for the
entire set of points. The performance for n generators in m dimensions is
<9(l_(m+I)/2j +n jQg ny ^ description can be found in [72].
Other higher-dimensional construction algorithms include a naive method, a
divide-and-conquer method, and several other incremental construction methods.
For a sample listing with references, see [65].
4.8 Algorithms for Constructing the Delaunay Tessellation
As mentioned earlier, the Delaunay tessellation is readily obtainable once the
Voronoi diagram is in hand. As we also noted, they are linear time transformable.
Although the algorithms of this chapter have focused on the Voronoi diagram, a
number of construction methods exist for the Delaunay tessellation. Some
Delaunay construction methods are a direct application of the very same
procedures used to construct the Voronoi diagram. The incremental, divide-and-
conquer, and plane sweep algorithms previously examined in this chapter are easily
adaptable to Delaunay triangulation construction. With a preference towards
initially producing the Delaunay tessellation, Fortune [38] applies these standard
construction methods to that end.
There is however, one interesting method specifically designed for the Delaunay
triangulation. This technique is alternately known as the flipping method [38] or
the swapping procedure [65]. The idea is that, for a given set of points in the
plane, initially produce any triangulation. Then examine in turn, each pair of
adjacent triangles. If a pair of triangles, say abc and acd, are not locally Delaunay,
flip (i.e.. swap) the diagonal of the quadrilateral formed by the triangle pair. For
example, if quadrilateral abed had diagonal ac it would be replaced by the opposite
diagonal bd. The condition of being locally Delaunay means to satisfy the min-max
angle criteria (see Property D8, in Chapter 2.2). With respect to performance, not
much is known about the number of flips required to transform a random
triangulation into the Delaunay triangulation. However, it has been determined that
at least rr flips are required for the transformation. This places a bound of 0(n2) on
the best we can hope for in the performance of the flipping method.
As a point of interest, we mention a Delaunay triangulation algorithm that,
although intolerably inefficient, can have a very concise code implementation.
ORourke shows in [67] how the lift up transformation (i.e., the theorem stated in
the previous section) can be used to construct the Delaunay triangulation in less
than 30 lines of code! Although, this particular implementation has an inferior
67


performance of 0(n4), it is an example of how deeper geometric insights can lead
to compact code.
68


5.
Applications of Voronoi Diagrams
Voronoi diagrams are used to model and solve many types of problems concerning
the spatial relationships between a set of points in ^-dimensional space. With an
extensive history of study and an astounding assortment applications, it is not
surprising that Voronoi diagrams have been used in a variety of wide-ranging
disciplines. The applications presented in this chapter are intended as an overview;
a survey that will help the reader appreciate the usefulness of these important
geometric structures.
The Voronoi diagram has a number of generalizations and related constructs. An
application may require a particular type of diagram to solve the problem at hand.
Moreover, implementation details can be complex and are often specific to the
type of structure used to model the application. It is beyond the scope of this thesis
to elaborate upon the lower-level details of the many types of applications; indeed,
this would comprise a lengthy volume in itself. Adequate references are provided
for the reader to obtain additional information about a particular application of
interest.
The need to model a diverse assortment of proximity problems has given rise to
many variations on the Voronoi diagram. Several generalizations of the Voronoi
diagram have proven quite useful in solving problems where the basic diagram fails
to provide sufficient information about the system being modeled. Related
structures such as the Delaunay tessellation, the Delaunay subgraphs, and the
convex hull, are often used when appropriate.
The applications presented are divided into three general categories: mathematics,
the natural sciences, and computer science. Each of these major areas contribute
different insights into the nature and use of Voronoi diagrams.
5.1 Applications in Mathematics
The concept of partitioning space into regions associated with sites is probably
much older than the first mathematical formalization of the subject. The fact that
some naturally occurring structures closely resemble Voronoi diagrams would lead
us to think that early scientists may have thought of the concept. In the
seventeenth century, Descartes published Principia Philosophiae which included a
69


drawing of his conception of the disposition of matter in the solar system.
Although this depiction bears a great resemblance to a Voronoi diagram, no
mention of this subject survives from his time.
The earliest known investigation of Voronoi diagrams was due to Dirichlet [25],
who in 1850 exploited the concept to prove the unique reducibility of quadratic
forms. In 1908, Voronoi [94] generalized Dirichlets result to higher dimensions.
As a tribute to their original work on the subject, the structure is most often
referred to as the Dirichlet tessellation or the Voronoi diagram. The subsequent
mathematical interest in Voronoi diagrams gave rise to a rich area of study which
laid the groundwork for the computer implementations and numerous applications
which were to follow. This section touches upon some of the major applications of
Voronoi diagrams that are a result of the mathematicians perspective on the
subject.
Voronoi Diagrams of Regularly Placed Generators
The first applications of Voronoi diagrams were in the study of geometrical
crystallography. This in turn motivated considerable mathematical interest into the
nature of how to uniformly tile space. A tiling of Rn is a covering of the 77-
dimensional space by closed sets, the interiors of which are pairwise disjoint. A
crystallographic group is a group of motions which map the tiling onto itself (i.e., a
homomorphism). The elements of this type of tiling are necessarily congruent
convex polytopes, known formally as stereohedra. The enumeration of all
stereohedra is one of the fundamental problems of crystallography and the
enumeration is still incomplete [6].
A subclass of stereohedra, characterized by regularly placed sites, is called
plesiohedra. Voronoi regions are well-suited to represent plesiohedra, and as a
powerful mathematical utility, they contributed to much progress in this field.
Delaunay et al. [23] provided a complete listing of all plesiohedra in R2 and defined
24 types of plesiohedra in M3. The geometric characteristics of plesiohedra can be
extremely complex, especially in higher dimensions. Engel [30] discovered three-
dimensional plesiohedra with 38 facets. The proven upper bound for the number of
facets of a plesiohedron in R3 is 390. To produce plesiohedra with a large number
of facets, regularly placed generators are chosen, followed by a computer
construction of parts of the resulting Voronoi diagram [30]. For more information
on the subject of tilings see [43] and [44].


Other applications of Voronoi diagrams to regularly placed sites include numerical
integration [8], sphere packing [75], and lattice systems [18].
Voronoi Diagrams of Irregularly Placed Generators
As noted above, regularly placed generators in ^-dimensional space correspond to
Voronoi diagrams composed of congruent convex polytopes. To the contrary,
generators with irregular placements form Voronoi diagrams composed of non-
congruent convex polytopes. Many problems in mathematics require modeling a
set of points with varying distance between adjacent points in the set. The notion
of applying Voronoi diagrams to solve these problems has motivated research in
several areas including the study of sphere packings. Rogers [75] investigated the
criteria for dense sphere packings. While doing so, the following question arose:
How large is the smallest Voronoi polytope with k facets and constructed on
generators having minimum distance two? Muder in [58] and [59] gives an
approximation to this problem in two dimensions and also provided results on
related quantities.
Considerable effort in the study of combinatorics has focused on the quantities of
the number of regions, edges, and vertices of Voronoi diagrams in various
dimensions as a function of the number of generators points. (We gave some of the
bounds on these relationships in Chapter 2.) The overall size of a Voronoi
diagram, in some sense, determines the amount of storage space needed to
represent the structure. The ratio of this quantity to the cardinality of the set of
input generators is an important number when implementations are required for
large diagrams. Preperata [71] showed that the size of a Voronoi diagram in!3 is
0(rr). The size of Voronoi diagrams in higher dimensions was completely analyzed
by Klee [50]. They determined that Voronoi diagrams in !" are equivalent to
certain convex hulls in Rn+1. The lift up method given in Chapter 4.7 of this paper
uses this very result to construct Voronoi diagrams in any dimension. This concept
led to the observation that various combinatorial properties including the size of an
(+1 )-dimensional structure apply directly to the ^-dimensional structure. Using
this notion, Seidel [80] gave exact bounds on the numbers of individual faces of
^-dimensional Voronoi diagrams.
Generalized Voronoi Diagrams
The Voronoi diagrams that we have examined thus far are called ordinary Voronoi
diagrams. Although planar and higher-dimensional structures were presented, the
diagrams were defined with a couple of conditions which result only in ordinary
71


Voronoi diagrams. Recall from the definition of the Voronoi diagram in Chapter 1,
that the Euclidean plane is partitioned so that every point in the plane is assigned
to its nearest generator point with respect to the Euclidean metric (i.e., our
conventional distance metric). If the generators are each assigned a particular
weight, and the Euclidean distance is replaced by a function which uses the weight
as a parameter, the resulting structure is known as a weighted Voronoi diagram.
The ordinary Voronoi diagram can be thought of as a weighted Voronoi diagram
with the usual Euclidean metric as the distance function. The distance function can
be the Euclidean metric weighted by a multiplicative constant to give what is called
a multiplicatively weighted Voronoi diagram. It may be defined by an additive
weighting to give an additively weighted Voronoi diagram. Combining the two
gives a compoundly weighted Voronoi diagram. When the distance function has as
parameters, the square of the Euclidean distance and the generator weight, the
resulting structure is called a power diagram.
Weighted Voronoi diagrams have been used in numerous applications, perhaps
more so than the ordinary kind. See [65] for an excellent mathematical description
of these structures and many examples of their use. Aurenhammer [5] provides
properties, construction algorithms, and applications for the power diagrams.
Another important family of generalized Voronoi diagrams are higher-order
Voronoi diagrams. An order-k Voronoi diagram for n sites is a partition of Mm into
regions such that any point within a fixed region has the same k closest generators.
In this definition, order indicates the number of points constituting the generator
and higher implies more than one point, not a higher dimensional space. We can
think of an order-A: Voronoi diagram in much the same way as an Voronoi ordinary
diagram except that the resulting regions can have any number (including 0) of
generators within the regions boundary. The regions are convex polytopes which
collectively cover the entire space of Mm. Ordinary Voronoi diagrams are those
with k = 1, that is, they are order-1 Voronoi diagrams.
Higher-order Voronoi diagrams have many properties not found in the ordinary
type of diagram. As with the class of weighted Voronoi diagrams, many
applications rely on order-& diagrams for a solution; the ordinary Voronoi diagram
may not adequately represent problems of considerable complexity.
A common application of order-fc Voronoi diagrams is the farthest-point problem.
Given a distinct set of n points in Mm, assign all locations in the space to the
farthest member(s) of the point set. The resulting structure tessellates Mm into
72


regions associated with elements of the point set. This partitioning of lRm is called
the farthest-point Voronoi diagram. Constructions were achieved in 0{n log n)
performance by Shamos [81], who used a divide-and-conquer scheme analogous
to the algorithm for the ordinary closest-point Voronoi diagram.
Higher-order Voronoi diagrams have been investigated by many researchers. Early
papers on the subject include Miles [56], Shamos and Hoey [82], and Lee [52].
Much is known about the characteristics of order-fc Voronoi diagrams in the plane.
On the other hand, the attributes of these complex structures in higher dimensions
remains one of the important open problems in combinatorial geometry.
Generalized Voronoi diagrams are not just limited to the weighed and higher-order
varieties. Voronoi diagrams can be constructed for spaces in which obstacles such
as polygons or other shapes replace the usual generator points. Other Voronoi
structures have been defined by altering the underlying space in which they live.
Voronoi diagrams on the surface of a sphere, cone, torus, and polyhedra have been
investigated (see [16] and [97]). On the more esoteric side of this research, the
behavior of Voronoi regions in Riemann manifolds was studied by Ehrlich and Im
Hof [31].
Stochastic Properties of Voronoi Diagrams
In our definition of the Voronoi diagram "V introduced in Chapter 1, the set of
generators P = {p,, p2,..., pn } had the conditions that n was finite and that all
generators were distinct. A generalization of "V can be made if the generators are
located in Mm according to the homogeneous Poisson point process. This structure
is known as a Poisson Voronoi diagram and its corresponding dual, the Poisson
Delaunay tessellation, denoted "VP and T)P respectively. Useful information
including distributions and correlations of various attributes such as the number of
sides, number of faces, area, volume, and total edge length may be derived from
these structures. These and other unique properties make the Poisson tessellations
an attractive modeling tool for certain planar and higher-dimensional applications
were random process plays an important role.
One of the earliest motivations for studying the stochastic properties of Poisson
Voronoi diagrams resulted from the relevance they have to chemical and physical
processes. As previously mentioned, crystallography and metallurgy extensively
exploited Voronoi diagrams; Poisson analysis of these structures has added
significant results to these fields. Meijering [55] used a Poisson field of points to
determine the volume, total boundary area, and the total edge length of the
73


polyhedral regions of a Voronoi diagram in M3. Variances for the polyhedral
regions and plane sections through the regions was given by Gilbert [42].
A number of researchers used empirical methods such as simulations to derive
results. Kiang [49] produced a formula for the random size distribution of Voronoi
regions in M1. Miles [56] provided a thorough study of Voronoi diagrams induced
by a homogeneous Poisson point process in the plane. Other combinatorial and
mensuration results were given by Crain [20] who observed frequencies of edges
per polygon and Dwyer [27] who gave the expected number of vertices of the
ordinary Voronoi diagram in n dimensions if the generators are uniformly
distributed in a hypersphere. Also, several results have also been obtained for
Voronoi diagrams with generators introduced in random order instead of drawn
from a probability distribution.
Subgraphs of the Delaunay Triangulation
The subgraphs of the Delaunay triangulation are of great mathematical importance
having a multitude of applications in path planning, combinatorial optimization,
network analysis, operations research, and, of course, graph theory. The most
fundamental subgraph of the Delaunay triangulation is the convex hull of its points
(see Property Dl). In Chapter 2 we defined three other notable subgraphs; the
Gabriel graph (see Gabriel and Sokal [40]), the relative neighborhood graph (see
Toussaint [90] and Urquhart [93]), and the Euclidean minimum spanning tree (see
Rohlb [76], Osteen and Lin [68], and Rosencrantz et al. [77]).
Another subgraph of the Delaunay triangulation that deserves mention is the
nearest neighbor graph, a structure that is used to solve the all nearest neighbors
problem (see Chapter 5.3 below). For a set of points, the nearest neighbor graph
(denoted NNG) has an edge between two points if and only if the endpoints (i.e.,
vertices) of the edge are nearest neighbors to each other. This graph is, in fact, a
subgraph of the Euclidean minimum spanning tree, which extends the hierarchy of
Delaunay subgraphs shown in Figure 2.6. Properties of the nearest neighbor graph
were investigated by Satty [78] and Pielou [70].
Several other graphs are closely related to the Delaunay triangulation and rely
upon it for suboptimal constructions. A Steiner tree is actually shorter than the
Euclidean minimum spanning tree if additional vertices are allowed. Garey,
Graham, and Johnson [41] showed it to be NP-complete and even with todays
computing technology, we are unable to solve problems with more than about 25
points [72]. The traveling salesperson problem (TSP) is perhaps the most well-
74


known graph theory problem because of its tremendous practical significance.
Although JVP-complete, it is easy to prove that a Euclidean minimum spanning tree
can be used to obtain an approximate TSP tour whose length is less than twice the
length of a shortest tour. A concise proof is given in [67].
5.2 Applications in the Natural Sciences
The natural sciences were the first areas of study to use Voronoi diagrams outside
of the pure mathematical domain in which they were discovered. During the early
part of this century, Voronoi diagrams were applied to crystallography,
meteorology, mining, and metallurgy. By mid-century, they modeled ecological
systems and geographic data, and were used in urban planning. Voronoi diagrams
are currently used in fields of research which range from the structure of molecules
to the distribution of galaxies in the universe.
Researchers in several fields of study have independently rediscovered the Voronoi
diagram. Consequently, the structure was given different names specific to the
respective discipline. Some of these alternate names are still in common use today.
The following applications are a few of the many to be found in the natural science
literature.
Crystallography
As mentioned above, the earliest application of the Voronoi diagrams was in the
study of crystallography, a subject dominated by German researchers around the
turn of this century. This came about from the initial developments of the Voronoi
diagram concept; defining the relationships between a set of regularly place points.
Niggli [61] referred to a Voronoi region as Wirkungsbereich (i.e., domain of
action, field of activity, or area of influence), a name still in common use today.
The application was an attempt to answer the most basic question of
crystallography: which domains of action completely fill R2 or R3 with congruent
regions? Articles by Delaunay [21] and Nowacki [62] presented significant work
on this subject.
Physicochemical Analysis
Equilibrium and other properties in physicochemical systems are a result of the
spatial distribution of distinct sites such as atoms or molecules. These properties
can be modeled by partitioning the space between the sites using nearest neighbor
analysis. The resulting structures are Voronoi diagrams. Wigner and Seitz [96]
75


used this concept in metallurgy, after which, the diagrams became known as
Wigner-Seitz zones. Spatial analysis was used by Frank and Kasper [39] to
investigate the composition of alloys, Brostow and Sicotte [13] to study the
behavior of liquid argon, Augenbaum and Peskin [4] in their research on hydro-
dynamic codes, and a number of other physicochemical researchers in related
fields.
Spatial Interpolation
Consider sites distributed over a finite surface (e.g., a country). To each site,
assign an observed numerical value (e.g., annual precipitation). Spatial
interpolation involves finding a function which best represents the entire surface
and which predicts values for locations other than the original set of data sites. If
we think of the sites as Voronoi generators, then the spatial interpolation can be
modeled as a finite generalized Voronoi diagram.
Meteorology and geography are two fields which often encounter problems
requiring spatial interpolation. In meteorology, Thiessen [89] used Voronoi
diagrams to determine estimates of the average precipitation over large areas. This
was perhaps the first use of Voronoi diagrams outside of the mathematical domain.
His technique was later studied by others who referred to it as Thiessens
method. Nowadays, Thiessen polygon is the preferred term over Voronoi polygon
in both meteorology and geography.
In the geographical sciences, Thiessens method has been applied to urban
planning by Snyder [85], to cartography by Arnold and Milne [2], and to facility
location by Okabe et al. [66]. Among the numerous applications in this broad field,
Suzuki and Iri [88] studied how to recover sites from a given subdivision of a
geographical area.
Applications in Ecology
In the field of ecology, plant and animal species are studied with respect to their
surroundings. The spatial patterns of the competition for space, food, and light are
a natural fit for Voronoi diagram representation. Brown [14] examined the density
of trees in a forest. In his research, an individual tree was assigned a region, called
the area potentially available to the tree. Brown had essentially rediscovered the
Voronoi diagram. Mead [54] used the same concept for plants in general and
denoted the Voronoi regions as plant polygons. Another paper by Mollison [57]
76


used Thiessen polygons to investigate ecological contact models and the spread of
epidemics.
5.3 Applications in Computer Science
By the middle part of this century, many problems had been identified to which
Voronoi diagrams offered a promising means of solution. However, the lack of
simple and efficient methods for constructing Voronoi diagrams limited their use in
many applications. The earliest diagrams were created by Euclidean construction,
that is, drawn with a straight-edge and compass. By the 1970s, progress in
computer construction and representation of Voronoi diagrams enabled
researchers to tackle applications that were formerly impractical.
For some problems such as nearest neighbor queries, the construction of the
Voronoi diagram is sufficient; the information required to solve the problem is
readily available in the representation of the structure itself. However, in many
applications the construction of the Voronoi diagram is followed by subsequent
processing and algorithms to produce the desired solutioa Efficient methods of
representation and construction were examined in Chapters 3 and 4. This section
outlines some of the classical problems of computer science that rely upon
constructed Voronoi diagrams for a solution. We present applications for the
planar case only; extensions to higher dimensions (where applicable) can be found
in the referenced articles.
Nearest Neighbor Query
Problems involving nearest neighbor relationships are the most fundamental and
useful applications of the Voronoi diagram. The simplest variation is the nearest
neighbor query also known as the post office problem. Given a set P of n sites in
the plane (post offices), locate a site closest to a given query point p (location of a
person). By computing and comparing all n distances, we get a trivial solution that
has a 0(n) performance. However, if this operation is repeated for a large number
of points Pi, the linear performance gives way to quadratic complexity. For
repeated queries, a more efficient approach is to first produce a Voronoi diagram
with generator set P. After the diagram is constructed in 0(n log n) time, point
locations within it are supported in 0(log n) time. Thus, the post office problem
can be solved with Voronoi diagrams in logarithmic query time. See Edelsbrunner
et al. [29] for an efficient solution to the point location problem.
77


Largest Empty Figures
In Chapter 2, we described the largest empty circle problem (see Property VI0). It
is stated as follows: Given a set P of n distinct points in the plane, find the largest
empty (i.e., contains no generators) circle whose center is in the interior of CH(P).
Shamos and Hoey [82] found a 0(n log n) solution to this problem by showing
that the center of such a circle is either a vertex of the Voronoi diagram of P or the
intersection of a Voronoi edge and the boundary of the CH(P). For an algorithm,
see ORourke [67].
This problem has important applications in facility location, operations research,
and industrial engineering. A typical application is to find a good location for a
new store, that is, to place it at the farthest distance from existing ones. Toussaint
[91] provides a unique application: to locate a nuclear reactor as far away from
centers of population as possible.
The largest empty circle problem can be generalized to include geometric figures
other than the circle. Given n points in the plane, place a figure of prescribed
shape (e.g., a rectangle) such that its area is maximized, but no point is covered by
the figure. An example of an application would be to find the largest rectangle
which can be salvaged from a flawed sheet of material such as fabric or wood.
Chazzelle et al. [17] provided a 0(n log3 n) solution for rectangles with sides
parallel to the planar axes.
A related problem worth mentioning is to find the smallest enclosing circle for a
set of points in the plane. Toussaint and Bhattacharya [92] proved that this circle is
centered at a vertex of the farthest point Voronoi diagram (see Chapter 5.1). As
with most Voronoi diagram applications in the plane, a 0(n log ri) solution exists.
Path Planning and Collision Avoidance
An important family of applications involves the motion of an object (e.g., a robot)
through a environment of obstacles. To minimize the risk of collision, the moving
object should remain as far away as possible from the obstacles during its travel
through the space. If the obstacles are points, then the ordinary Voronoi diagram
will suffice. If the obstacles are polygons or other configurations, then a
generalized Voronoi diagram is needed to determine the path with the least
likelihood of a collision. A further refinement on these type of problems deals with
finding the shortest path through the set of obstacles. For an excellent overview of
collision avoidance and path planning see [67]. We continue our examination of
78


this topic by describing a couple of fundamental problems in computational
geometry.
The all nearest neighbor problem was briefly mentioned in Chapter 1 and related
to the nearest neighbor graph (NNG) in Chapter 5.1. It is stated as follows: given a
set P of n points in the plane, find a nearest neighbor of each. The construction of
the NNG for the set P provides a 0(n log n) solution as the NNG is obtained in
linear time from the Delaunay triangulation on P. The all nearest neighbor problem
is used in collision detection, point pattern analysis (see Boots and Getis [12]), and
other areas of computer science.
Another problem that addresses one of the most fundamental questions in
computational geometry is the closest pair problem. Closest pair problems have
considerable relevance in collision detection applications such as robot motion and
air traffic control. It is easily stated as follows: given a set P of n points in the
plane, find two points, say pt and pj, whose mutual distance is the smallest among
all pairs of points in P. Clearly, a quadratic time solution can be accomplished by
considering each pair of points, but this is suboptimal. It is obvious that any pair of
closest points must be in consecutive sorted order to give an optimal 0(n log n)
solution. An alternative to this approach is to construct the Voronoi diagram for
set P, produce a nearest neighbor graph (see above), and perform a comparison of
the edge lengths. This also achieves the optimal time complexity of 0(n log n) as
shown by Shamos and Hoey [82].
It should be evident that the motion of a disk (e.g., a circular robot) between two
given points in a region of obstacles is simply the generalized Voronoi diagram for
the obstacles. The translational motion of a convex planar robot (i.e., a convex
polygonal cross section) was examined by Leven and Sharir [53] who determined a
0(n log n) performance for constructing the path. When rotational motion is
included, a third degree of freedom must be added to the model as in the example
of a line segment (e.g., a ladder) that is allowed to rotate during the travel. For this
type of complex problem, a Voronoi diagram can be defined such that at any time,
the line segment is equidistant from its two closest obstacles. In this case closeness
is with respect to the minimum distance between the segment and the obstacle. An
investigation of this problem was undertaken by Sifrony and Sharir [84] wherein it
was shown that this type of diagram could be constructed (and a path planned) in
0(t log n), where t = CHji1) is a parameter which represents the length of the line
segment. Schwartz and Yap [79], and Alt and Yap [1] provide an introduction to
this subject and cover the use of Voronoi diagrams in motion planning.
79


When the object is taken to be a single moving point, we have a special case of
collision-free motions called path planning. Considerable research has been
devoted to finding the shortest path through a region of polygonal obstacles.
One useful application is to determine the locations, for each of m fixed sites
within a polygonal factory floor, that is reachable along a shortest path.
Aronov [3] found that by using a generalized Voronoi diagram for n obstacles,
a 0((n + rri) log2 (n + mj) performance could be achieved.
80


References
[1] Alt, H. and Yap, C.K. (1990) Algorithmic aspects of motion planning: A tutorial (Part
2). Algorithms Review, 1, 61-77.
[2] Arnold, D.B. and Milne, W.J. (1984) The use of Voronoi tessellations in processing soil
survey results. IEEE Computer Graphics Applications, 4, 22-30.
[3] Aronov, B. (1989) On the geodesic Voronoi diagram of point sites in a simple polygon.
Algorithmica, 4, 109-140.
[4] Augenbaum, J.M. and Peskin, C.S. (1985) On the construction of the Voronoi mesh on
a sphere. Journal of Computational Physics, 59, 177-192.
[5] Aurenhammer, F. (1985) Power diagrams: properties, algorithms, and applications.
SIAM Journal on Computing, 16(1), 78-96.
[6] Aurenhammer, F. (1991) Voronoi diargams: a survey of a fundamental data structure.
ACM Computing Surveys, 23(3), 345-405.
[7] Aurenhammer, F. and Edelsbrunner, H. (1984) An optimal algorithm for constructing
the weighted Voronoi diagram in the plane. Pattern Recognition, 17,251-257.
[8] Babenko. V.F. (1977) On the optimal cubature formulae on certain classes of continuous
functions. Analysis Mathematics., 3, 3-9.
[9] Ball. W.W.R. (1962) Mathematical Recreations and Essays, 11th edition, Macmillan,
New York. NY.
[10] Ballard. D.H. and Brown, C.M. (1982) Computer Vision, Prentice-Hall, Englewood
Cliffs. NJ.
[11] Baumgart, B. (1975) A polyhedron representation for computer vision. IFIPS
Conference Proceedings, 44, 589-596.
[12] Boots, B.N. and Getis, A. (1988) Point Pattern Analysis, Sage Scientific Geography
Series, Volume 8, Sage Publications, Newbury Park, CA.
[13] Brostow, W. and Sicotte, Y. (1975) Coordination number in liquid argon. Physica, 80A,
513-522.
[14] Brown, G.S. (1965) Point density in stems per acre. New Zealand Forestry Service
Research Notes, 38, 1-11.
81


[15] Brown, K.Q. (1979) Voronoi diagrams from convex hulls. Information Processing
Letters, 9(5), 223-228.
[16] Brown, K.Q. (1980) Geometric transforms for fast geometric algorithms. Ph.D. Thesis,
Report CMU-CS-80-101, Camegie-Mellon University, Pittsburgh, PA.
[17] Chazelle, B., Drysdale, R.L. and Lee, D.T. (1986) Computing the largest empty
rectangle. SIAM Journal on Computing, 15, 300-315.
[18] Conway, J.M. and Sloane, NJ.A. (1982) Voronoi regions of lattices, second moments of
polytopes, and quantization. IEEE Transactions on Information Theory, 28, 211-226.
[19] Coxeter, H.S.M. (1973) Regular Polytopes, 3rd edition, Dover, New York.
[20] Crain, I.K. (1978) The Monte Carlo generation of random polygons. Computers and
Geosciences, 4, 131-141.
[21] Delaunay, B.N. (1932) Neue Darstellung der geometrischen. Kristallographie Zeitschrift
fur Kristallographie, 84, 109-149.
[22] Delaunay, B.N. (1934) Sur la sphere vide. Bulletin of the Academy of Sciences of the
U.S.S.R., Classe des Sciences Mathematiques et Naturelles, 7(6), 793-800.
[23] Delaunay, B.N., Dolbilin, N.P. and Stogrin, M.I. (1980) Combinatorial and metric
theory of planigons. Proceedings of the Steklov Institute of Mathematics, Issue 4, 109-
MO.
[24] Dillencourt, M.B. (1990) Toughness and Delaunay triangulations. Discrete and
Computational Geometry, 5, 573-601.
[25] Dirichlet, G.L. (1850) Uber die reduction der positeven quadratischen formen mit drei
unbestimmten ganzen zahlen. Journal fur die Reine und Angewandte Mathematik, 40,
209-227.
[26] Dwyer, R.A. (1987) A faster divide-and-conquer algorithm for constructing Delaunay
triangulations. Algorithmica, 2, 137-151.
[27] Dwyer, R.A. (1989) Higher-dimensional Voronoi diagrams in linear expected time.
Proceedings of the 5th Annual Symposium on Computational Geometry, 326-333.
[28] Edelsbrunner, H. and Miicke, E.P. (1987) Simulation of simplicity: a technique to cope
with degenerate cases in geometric algorithms. Proceedings of the 4th Annual ACM
Symposium on Computational Geometry, 118-133.
[29] Edelsbrunner, H., ORourke, J. and Seidel, R. (1986) Constructing arrangements of
lines and hyperplanes with applications. SIAM Journal on Computing, 15, 341-363.
82


[30] Engel, P. (1981) Uber Wirkungsbereichsteilungen mit kubischer Symmetric. Zeitschrift
fur Kristallographie, 154, 199-215.
[31] Ehrlich, P.E. and Im Hof, H.C. (1979) Dirichlet regions in manifolds without conjugate
points. Commentarii Mathematici Helvetici, 54, 642-658.
[32] Fang, T-P. and Piegl, L.A. (1995) Delaunay trianglulations in three dimensions. IEEE
Computer Graphics and Applications, 15, 62-69.
[33] Fary, I. (1948) On straight-line representation of planer graphs. Acta Sci. Math.,
Szeged, 11,229-233.
[34] Fortune, S. (1986) A sweepline algorithm for Voronoi diagrams. Proceedings of the 2nd
Annual ACM Symposium on Computational Geometry, 313-322.
[35] Fortune, S. (1987) A sweepline algorithm for Voronoi diagrams. Algorithmica, 2, 153-
174.
[36] Fortune, S. (1989) Stable maintenance of point-set triangulation in two dimensions.
Proceedings of the 3ffh Annual Symposium on the Foundations of Computer Science,
494-499.
[37] Fortune, S. (1991) Numerical stability of algorithms for Delaunay triangulations and
Voronoi diagrams in two dimensions. Manuscript, AT&T Bell Labs.
[38] Fortune, S. (1992) Voronoi diagrams and Delaunay triangulations. Computing in
Euclidean Geometry, Lecture Notes Series on Computing, Vol. 1, World Scientific, NJ,
193-233.
[39] Frank, F.C. and Kasper, J.S. (1958) Complex alloy structures regarded as sphere
packings. Acta Crystallographica, 11, 184-190.
[40] Gabriel, K.R and Sokal, R.R. (1969) A new statistical approach to geographic variation
analysis. Systematic Zoology, 18, 259-278.
[41] Garey, M.R., Graham, R.L. and Johnson, D.S. (1977) The complexity of computing
Steiner minimal spanning trees. SIAM Journal of Applied Mathematics, 32, 835-859.
[42] Gilbert, E.N. (1962) Random subdivisions of space into crystals. Annual Mathematical
Statistics, 33, 958-972.
[43] Griinbaum, B. and Shephard, G.C. (1980) Tilings with congruent tiles. Bulletin of the
American Mathematical Society, 3, 951-973.
[44] Griinbaum, B. and Shephard, G.C. (1987) Tilings and Patterns, Freeman and Co., New
York, NY.
83


[45]
Guibas, L. and Stolfi, J. (1985) Primitives for the manipulation of general subdivisions
and the computation of Voronoi diagrams. ACM Transactions on Graphics, 4(2), 74-
123.
[46] Hoffmann, C.M. (1989) Geometric and Solid Modeling: An Introduction, Morgan
Kaufinann, San Mateo, CA.
[47] Karasick, M., Lieber, D. and Nackman, L. (1990) Efficient Delaunay triangulations
using rational arithmetic. ACM Transactions on Graphics, 10(1), 71-91.
[48] Katajainen, J. and Koppinen, M. (1988) Constructing Delaunay triangulations by
merging buckets in quadtree order. Fundamenta Informaticae, 11, 275-288.
[49] Kiang, T. (1966) Random fragmentation in two and three dimensions. Zeitschriftfur
Astrophysik, 64, 433-439.
[50] Klee, V. (1980) On the complexity of ^-dimensional Voronoi diagrams. Archiv de
Mathematik, 34, 75-80.
[51 ] Latakos. 1.(1988) Proofs and Refutations, Cambridge University Press, Cambridge,
England.
[52] Lee. D.T. (1982) On ^-nearest neighbor Voronoi diagrams in the plane. IEEE
Transactions on Computing, 31, 478-487.
[53] Leven. D. and Sharir, M. (1987) Planning a purely translational motion for a convex
object in two-dimensional space using generalized Voronoi diagrams. Discrete and
Computational Geometry, 2, 9-31.
[54] Mead. R. (1966) A relationship between individual plant-spacing and yield. Annals of
Botany. N.S., 30, 301-309.
[55] Meijering, J.L. (1953) Interface area, edge length, and number of vertices in crystal
aggregates with random nucleation. Philips Research Report 8,270-290.
[56] Miles. R.E. (1970) On the homogenous planer Poisson process. Math. Biosci, 6, 85-127.
[57] Mollison, D. (1977) Spatial contact models for ecological and epidemic spread. Journal
of the Royal Statistical Society B, 39, 283-326.
[58] Mudder, D.J. (1988a) How big is an w-sided Voronoi polygon?, Manuscript, MITRE
Corp., Bedford, MA.
[59] Mudder, D.J. (1988b) Putting the best face on a Voronoi polyhedron. Proceedings of the
London Mathematical Society, 56, 329-348.
[60] Muller, D.E. and Preperata, F.P. (1978) Finding the intersection of two convex
polyhedra. Theoretical Computer Science, 7(2), 217-236.
84


[61 ] Niggli, R. (1927) Die topologische Strukturanalyse. Zeitschrift fur Kristallographie, 69,
391-415.
[62] Nowacki, V.W. (1976) Uber allgemeine Eigenschaften von Wirkungsbereichen.
Zeitschrift fur Kristallographie, 143, 360-385.
[63] Ohya, T., Iri, M. and Murota, K. (1984a) Improvements of the incremental method for
the Voronoi diagram with computational comparisons of various algorithms. Journal of
the Operations Research Society of Japan, 27, 306-336.
[64] Ohya, T., Iri, M. and Murota, K. (1984b) A fast Voronoi-diagram algorithm with
quaternary tree bucketing. Information Processing Letters, 18(4), 227-231.
[65] Okabe, A., Boots, B. and Sugihara, K. (1992) Spatial Tessellations: Concepts and
Applications of Voronoi Diagrams, John Wiley, Chichester, England.
[66] Okabe, A., Yoshikawa, T., Fujii, A. and Oikawa, K. (1988) The statistical analysis of a
distribution of activity points in relation to surface-like elements. Environment and
Planning A, 20, 609-620.
[67] ORourke, J. (1993) Computational Geometry in C, Cambridge University Press,
Cambridge, England.
[68] Osteen, R.E. and Lin, P.P. (1974) Picture skeletons based on eccentricities of points of
minimum spanning trees. SIAM Journal on Computing, 3, 23-40.
[69] Overmars, M.H. and van Leeuwen, J. (1981) Maintenance of point configurations in the
plane. Journal of Computers and System Sciences, 23, 116-204.
[70] Pielou, E. (1977) Mathematical Ecology, Wiley-Interscience, New York, NY.
[71] Preperata, F.P. (1977) Steps into computational geometry. Report R-760, Coordinated
Science Lab., University of Illinois, Urbana, IL, 23-24.
[72] Preperata, F.P. and Shamos, M.I. (1985) Computational Geometry: An Introduction,
Springer-Verlag, New York, NY.
[73] Rajan, V.T. (1991) Optimality of the Delaunay triangulation in Md. Proceedings of the
7th Annual ACM Symposium on Computational Geometry, 357-372.
[74] Roberts, F.S. (1984) Applied Combinatorics, Prentice-Hall, Englewood Cliffs, New
Jersey, NJ.
[75] Rogers, C.A. (1964) Packing and Covering, Cambridge University Press, Cambridge,
England.
[76] Rohlb, F.J. (1973) Hierarchical clustering using the minimum spanning tree. The
Computer Journal, 16, 93-95.
85


[77] Rosencrantz, D. J., Stearns, R.E. and Lewis II, P.M. (1977) An analysis of several
heuristics for the traveling salesman problem. SIAM Journal on Computing, 6, 563-581.
[78] Satty, T.L. (1970) Optimization in Integers and Related Extremal Problems, McGraw-
Hill, New York, NY.
[79] Schwartz, J. and Yap, C.K. (1986) Advances in Robotics, Lawrence Erlbaum
Associates, Hillside, NJ.
[80] Seidel, R. (1982) The complexity of Voronoi diagrams in higher dimensions.
Proceedings of the 2ffh Annual Allerton Conference on CCC, 94-95.
[81] Shamos, M.I. (1978) Computational Geometry, Ph.D. dissertation. Yale University, New
Haven, CT.
[82] Shamos, M.I. and Hoey, D. (1975) Closest-point problems. Proceedings of the 16th
Annual IEEE Symposium on Foundations of Computer Science, 151-162.
[83] Sibson, R. (1978) Locally equiangular triangulations. The Computer Journal, 21, 243-
245.
[84] Sifrony, S. and Sharir, M. (1986) A new efficient motion-planning algorithm for a rod
in polygonal space. Proceedings of the 2nd Annual Symposium on Computational
Geometry, 178-186.
[85] Snyder, D.E. (1962) Urban places in Uraguay and the concept of a hierarchy.
Northwestern University Studies in Geography, 6, 29-46.
[86] Sugihara, K. and Iri, M. (1988) Geometric algorithms in finite-precision arithmetic.
Abstracts of the 13lh International Symposium on Mathematical Programming, 196.
[87] Sugihara, K. and Iri, M. (1989) VORONOI2 reference manual: Topology-oriented
version for the incremental method for constructing Voronoi diagrams. Research
Memorandum RMI 89-04, Department of Mathematical Engineering and Information
Physics, University of Tokyo, Tokyo, Japan.
[88] Suzuki, A. and Iri, M. (1986) Approximation of a tessellation of the plane by a Voronoi
diagram. Journal of the Operations Research Society ofJapan, 29, 69-96.
[89] Thiessen, A. H. (1911) Precipitation averages for large areas. Monthly Weather Review,
39, 1082-1084.
[90] Toussaint, G.T. (1980) The relative neighborhood graph of a finite plane set. Pattern
Recognition, 12,261-268.
[91] Toussaint, G.T. (1983) Computing largest empty circles with location constraints.
International Journal on Computer Information Science, 12, 347-358.
86


[92] Toussaint, G.T. and Bhattacharya, B.K. (1981) On geometric algorithms that use the
farthest-point Voronoi diagram. Report 81-3, McGill University, Montreal, Quebec,
Canada.
[93] Urquhart, R. (1980) A note on the computation of subgraphs of the Delaunay
triangulation. Report, University of Glasgow, U.K.
[94] Voronoi, G. (1908) Nouvelles applications des parametres continus a la theorie des
formes quadratiques deuxieme memoire, recherches sur les parallelloedres primitifs.
Journal fur die Reine und Angewandte Mathematik, 134, 209-227.
[95] Watson, D.F. (1981) Computing the n-dimensional Delaunay tessellation with
application to Voronoi polytopes. The Computer Journal, 24(2), 198-287.
[96] Wigner, E. and Seitz, F. (1933) On the constitution of metallic sodium. Physical
Review, 43, 804-810.
[97] Yap, C.K. (1987) An 0(n log n) algorithm for the Voronoi diagram of a set of simple
curve segments. Discrete and Computational Geometry, 2, 365-393.
87