Lecture 16: Fortune's Algorithm and Voronoi diagrams

Properties of the Voronoi diagram: 

We first revisit some observations about the structure of Voronoi
diagrams in the plane.

Voronoi edges: Each point on an edge of the Voronoi diagram is
equidistant from its two nearest neighbors p_i and
p_j. Thus, there is a circle centered at such a point such
that p_i and p_j lie on this circle, and no other
site is interior to the circle.

Voronoi vertices: It follows that vertex at which three Voronoi cells
V(p_i), V(p_j), and V(p_k) intersect is
equidistant from all sites. Thus it is the center of the circle
passing through these sites, and this circle contains no other sites
in its interior.

Degree: If we assume that no four points are cocircular, then the
vertices of the Voronoi diagram all have degree three.

Convex hull: A cell of the Voronoi diagram is unbounded if and only if
the corresponding site lies on the convex hull. (Observe that a point
is on the convex hull if and only if it is the closest point from some
point at infinity.)

Size: If n denotes the number of sites, then the Voronoi diagram is a
planar graph (if we imagine all the unbounded edges as going to a
common vertex infinity) with exactly n faces. It follows from Euler's
formula that the number of Voronoi vertices is at most 2n - 5 and the
number of edges is at most 3n - 6. 

Computing Voronoi Diagrams: There are a number of algorithms for
computing Voronoi diagrams. Of course, there is a naive
O(n² log n) time algorithm, which operates by computing
V(p_i) by intersecting the n - 1 bisector halfplanes
h(p_i, p_j), for j != i. However, there are much
more efficient ways, which run in O(n log n) time. Since the convex
hull can be extracted from the Voronoi diagram in O(n) time, it
follows that this is asymptotically optimal in the worst case.

Historically, O(n²) algorithms for computing Voronoi
diagrams were known for many years (based on incremental
constructions). When computational geometry came along, a more
complex, but asymptotically superior O(n log n) algorithm was
discovered. This algorithm was based on divide and conquer. But it was
rather complex, and somewhat difficult to understand.  Later, Steven
Fortune invented a plane sweep algorithm for the problem, which
provided a simpler O(n log n) solution to the problem. It is his
algorithm that we will discuss. Somewhat later still, it was
discovered that the incremental algorithm is actually quite efficient,
if it is run as a randomized incremental algorithm. We will discuss
this algorithm later when we talk about the dual structure, called a
Delaunay triangulation.

Fortune's Algorithm: Before discussing Fortune's algorithm, it is
interesting to consider why this algorithm wasn't invented much
earlier. In fact, it is quite a bit trickier than any plane sweep
algorithm we have seen so far. The key to any plane sweep algorithm is
the ability to discover all ``upcoming'' events in an efficient
manner. For example, in the line segment intersection algorithm we
considered all pairs of line segments that were adjacent in the sweep
line status, and inserted their intersection point in the queue of
upcoming events. The problem with the Voronoi diagram is that of
predicting where the upcoming events will occur. The reason is that a
site that lies ahead of the sweep line may generate a Voronoi vertex
that lies behind the sweep line. It is these ``unanticipated events''
that make the design of a plane sweep algorithm challenging.

Fortune made the clever observation of rather than computing the
Voronoi diagram through plane sweep in its final form, instead to
compute a ``distorted'' but topologically equivalent version of the
diagram. This distorted version of the diagram was based on a
transformation that alters the way that distances are measured in the
plane. The resulting diagram had the same structure as the Voronoi
diagram, but its edges were parabolic arcs, rather than straight line
segments. Once this distorted diagram was generated, it was an easy
matter to ``undistort'' it to produce the correct Voronoi diagram.

Our presentation will be different from Fortune's. Rather than distort
the diagram, we can think of this algorithm as distorting the sweep
line. Actually, we will think of two objects that control the sweeping
process. First, there will be a horizontal sweep line, moving from top
to bottom. We will also maintain a x monotonic curve called a beach
line (I guess because it looks like waves rolling up on a beach.) The
beach line is formed from parabolic arcs. As the sweep line moves
downward, the beach line moves downward as well, but the beach line's
shape depends on the locations of the sites. We will see that the
Voronoi diagram ``grows out'' of the beach line.


Figure 57: Plane sweep Voronoi diagrams.

The Beach Line: In order to make these ideas more concrete, recall
that the problem with ordinary plane sweep is that sites that lie
below the sweep line may affect the diagram that lies above the sweep
line. To avoid this problem, we will maintain only the portion of the
diagram that cannot be affected by anything that lies below the sweep
line. To do this, we will ``cut'' the halfplane lying above the sweep
line into two regions: those points that are closer to some site p
above the sweep line than they are to the sweep line itself, and those
points that are closer to the sweep line than any site above the sweep
line. The set of points q that are equidistant from the sweep line to
their nearest site above the sweep line is called the beach line.

(Hey, you, reader!  THIS NEXT PARAGRAPH IS THE WHOLE TRICK.
everything else is just implementation)

Observe that for any point q above the beach line, we know that its
closest point cannot be affected by any site that lies below the sweep
line. Hence, the portion of the Voronoi diagram that lies above the
beach line is ``safe'' in the sense that we have all the information
that we need in order to compute it (without knowing abouts what
points are still to appear below the sweep line).

(GOT IT?, ok, now go on...)

What does the beach line look like? We know from high school geometry
that the set of points that are equidistant from a site lying above a
horizontal line and the line itself forms a parabola that is open on
top. With a little analytic geometry, it is easy to show that the
parabola becomes ``skinnier'' as the site becomes closer to the line,
the parabola degenerates into a vertical ray emanating from the
site. (You should work out the equations to see why this is so.)

Thus, the beach line consists of the lower envelope of these
parabolas. (By the way, if we were interested in computing the Voronoi
diagram of line segments, the beach line is the boundary of the
Voronoi cell of the sweep line itself.) Because the parabolas are x
monotone, so is the beach line. Also observe that the vertex where two
arcs of the beach line intersect, which we call a breakpoint, is a
point that is equidistant from two sites and the sweep line, and hence
must lie on some Voronoi edge. In particular, if the beach line arcs
corresponding to points p_i and p_j share a common
breakpoint on the beach line, then this breakpoint lies on the Voronoi
edge between p_i and p_j. From this we have the
following important characterization.

Lemma: The beach line is an x monotone curve made up of parabolic
arcs. The breakpoints of the beach line lie on Voronoi edges of the
final diagram.  Fortune's algorithm consists of simulating the growth
of the beach line as the sweep line moves downward, and in particular
tracing the paths of the breakpoints as they travel along the edges of
the Voronoi diagram. Of course, as the sweep line moves the parabolas
forming the beach line change their shapes continuously. As with all
plane sweep algorithms, we will be interested in simulating the
discrete event points where there is a ``significant event'', that is,
any event that changes the topological structure of the beach line. It
turns out these significant events will be of two varieties:

Site events: When the sweep line passes over a new site a new arc will
be inserted into the beach line.

Vertex events: (Sometimes called circle events.) When the length of
a parabolic arc shrinks to zero, the arc disappears and a new Voronoi
vertex will be created at this point.

The algorithm consists of processing these two types of events. As the
Voronoi vertices are being discovered by vertex events, it will be an
easy matter to update a DCEL for the diagram as we go, and so to link
the entire diagram together. Next time we will show how these two
events are discovered and processed.


Fortune's Algorithm finally

Fortune's Algorithm: Last time we introduced Voronoi diagrams, which
for given a set of sites partitions the plane into regions according
to which site is the closest. We also introduced Fortune's algorithm,
a plane sweep algorithm for computing Voronoi diagrams in the plane.
Recall that (because of the unanticipated events) rather than just
sweeping a straight line, Fortune's algorithm sweeps two objects, the
standard horizontal sweep line and curved line called the beach
line. The beach line is the set of points that are equidistant from
some virtual sweep line and the nearest site. It consists of pieces of
parabolic arcs, joined end to end at points called breakpoints. The
breakpoints lie on the Voronoi diagram.


Fortune's algorithm simulates the growth of the beach line as the
sweep line moves downward, and in doing so, traces the paths of the
breakpoints as they travel along the edges of the Voronoi
diagram. Last time we noted that there are two types of events, site
events when the horizontal sweep line sweeps over a new site and
vertex events when two breakpoints on the beach line collide with each
other, resulting in the generation of a vertex of the Voronoi
diagram. Today we discuss how to handle these two types of events.


Site events: A site event is generated whenever the horizontal sweep
line passes over a site. As we mentioned before, at the instant that
the sweep line touches the point, its associated parabolic arc will
degenerate to a vertical ray shooting up from the point to the current
beach line. As the sweep line proceeds downwards, this ray will widen
into an arc along the beach line. To process a site event we will
determine the arc of the sweep line that lies directly above the new
site. (Let us make the general position assumption that it does not
fall immediately below a vertex of the beach line.) We then split this
arc of the beach line in two by inserting a new infinitesimally small
arc at this point. As the sweep proceeds, this arc will start to
widen, and eventually will join up with other edges in the
diagram. (See the figure below.)


Figure 58: Site events.

It is important to consider whether this is the only way that new arcs
can be introduced into the sweep line. In fact it is. (We will not
prove it). As a consequence of this proof, it follows that the maximum
number of arcs on the beach line can be at most 2n - 1, since each new
point can result in creating one new arc, and splitting an existing
arc, for a net increase of two arcs per point (except the first).

The nice thing about site events is that they are all known in
advance. Thus, after sorting the points by y coordinate, all these
events are known.

Vertex events: In contrast to site events, vertex events are
generated dynamically as the algorithm runs. As with the line segment
plane sweep algorithm, the important idea is that each such event is
generated by objects that are neighbors on the beach line. However,
unlike the segment intersection where pairs of consecutive segments
generated events, here triples of points generate the events.

In particular, consider any three consecutive points p_i,
p_j, and p_k whose arcs appear consecutively on
the beach line from left to right. (See the figure below.) Further,
suppose that the circumcircle for these three sites lies at least
partially below the current sweep line (meaning that the Voronoi
vertex has not yet been generated), and that this circumcircle
contains no points lying below the sweep line (meaning that no future
point will block the creation of the vertex).


Consider the moment at which the sweep line falls to a point where it
is tangent to the lowest point of this circle. At this instant the
circumcenter of the circle is equidistant from all three sites and
from the sweep line. Thus all three parabolic arcs pass through this
center point, implying that the contribution of the arc from
p_j has disappeared from the beach line. In terms of the
Voronoi diagram, the bisectors (p_i, p_j) and
(p_j, p_k) have met each other at the Voronoi
vertex, and a single bisector (p_i, p_k) remains.


Figure 59: Vertex events.

Sweep line algorithm: We can now present the algorithm is greater
detail. The main structures that we will maintain are the following:


(Partial) Voronoi diagram: The partial Voronoi diagram that has been
constructed so far will be stored in a DCEL. There is one technical
difficulty caused by the fact that the diagram contains unbounded
edges. To handle this we will assume that the entire diagram is to be
stored within a large bounding box. (This box should be chosen large
enough that all of the Voronoi vertices fit within the box.)


Beach line: The beach line is represented using a dictionary
(e.g. a balanced binary tree or skip list). An important fact of the
construction is that we do not explicitly store the parabolic
arcs. They are just their for the purposes of deriving the
algorithm. Instead for each parabolic arc on the current beach line,
we store the site that gives rise to this arc. Notice that a site may
appear multiple times on the beach line (in fact linearly many times
in n). But the total length of the beach line will never exceed 2n-1.

Between each consecutive pair of sites p_i and
p_j, there is a breakpoint. Although the breakpoint moves as
a function of the sweep line, observe that it is possible to compute
the exact location of the breakpoint as a function of p_i,
p_j, and the current y coordinate of the sweep line. Thus,
as with beach lines, we do not explicitly store breakpoints. Rather,
we compute them only when we need them.

The important operations that we will have to support on the beach line are

(1) Given a fixed location of the sweep line, determine the arc of the
beach line that intersects a given vertical line. This can be done by
a binary search on the breakpoints, which are computed ``on the
fly''. (Think about this.)

(2) Compute predecessors and successors on the beach line.

(3) Insert an new arc p_i within a given arc p_j,
thus splitting the arc for p_j into two. This creates three
arcs, p_j, p_i, and p_j.

(4) Delete an arc from the beach line.

It is not difficult to modify a standard dictionary data structure to
perform these opera tions in O(log n) time each.

Event queue: The event queue is a priority queue with the ability both
to insert and delete new events. Also the event with the largest y
coordinate can be extracted. For each site we store its y coordinate
in the queue.


For each consecutive triple p_i, p_j,
p_k on the beach line, we compute the circumcircle of these
points. (We'll leave the messy algebraic details as an exercise, but
this can be done in O(1) time.) If the lower endpoint of the circle
(the minimum y coordinate on the circle) lies below the sweep line,
then we create a vertex event whose y coordinate is the y coordinate
of the bottom endpoint of the circumcircle. We store this in the
priority queue. Each such event in the priority queue has a cross link
back to the triple of sites that generated it, and each consecutive
triple of sites has a cross link to the event that it generated in the
priority queue.

The algorithm proceeds like any plane sweep algorithm. We extract an
event, process it, and go on to the next event. Each event may result
in a modification of the Voronoi diagram and the beach line, and may
result in the creation or deletion of existing events.  Here is how
the two types of events are handled:

Site event: Let p_i be the current site. We shoot a vertical
ray up to determine the arc that lies immediately above this point in
the beach line. Let p j be the corresponding site.  We split this arc,
replacing it with the triple of arcs p_j, p_i,
p_j which we insert into the beach line. Also we create new
(dangling) edge for the Voronoi diagram which lies on the bisector
between p_i and p_j. Any old triple that involved
p_j as its center arc is deleted from the priority queue,
and we generate new events for each of the possible three new triples
that result from this insertion.

Vertex event: Let p_i, p_j, and p_k be
the three sites that generate this event (from left to right). We
delete the arc for for p j from the beach line. We create a new vertex
in the Voronoi diagram, and tie the edges for the bisectors
(p_i, p_j), (p_j,p_k ) to it,
and start a new edge for the bisector (p_i, p_k)
that starts growing down below. Finally, we delete any events that
arose from triples involving this arc of p_j, and generate
new events corresponding to consecutive triples involving
p_i and p_k (there are two of them).


The analysis follows a typical analysis for plane sweep. Each event
involves O(1) processing time plus a constant number accesses to the
various data structures. Each of these accesses takes O(logn) time,
and the data structures are all of size O(n). Thus the total time is
O(n log n), and the total space is O(n).

Applets:


nice
animation but fixed points.