In this chapter we will deal with the problem known as the "post office problem".
It is about dividing the territory of the city into areas around individual post offices so that in each area the local post office is the closest. Another example is the determination of the catchment areas of hospitals by distance. The mathematical formulation of the problem is as follows:
Let $P=\{p_1,p_2,\dots,p_n\}$ be a set of points in the plane. The Voronoi diagram
(briefly V-diagram) is a plane subdivision which has $n$ faces
$$V(p_i)=\{q\in \R^2;\ \d(q,p_i)\leq \d(q,p_j) \text{ for all }j\in \{1,\dots,n\}, j\ne i \},$$
where $\d(q,p)$ denotes the classical Euclidean distance of points $q$ and $p$.
Figure 9.1 An example of the Voronoi diagram.
The set of points in the plane having the distance to a point $ p_i $ less than or equal to the distance to the point $p_j \ne p_i$ is a halfplane
The face $V(p_i)$, called the Voronoi cell, is the intersection of all such halfplanes for $j\ne i$
$$V(p_i)=\bigcap_{j\ne i}h(p_i,p_j).$$
According to Chapter 5 we can find the intersection of $n-1$ halfplanes in time
$O(n\log n)$. To find all Voronoi cells by this way would need time at least
$O(n^2\log n)$. We will show an algorithm for finding the Voronoia diagram with running time $O(n\log n)$. The fact that such running time is real is supported by the following statement:
Theorem 9.1
Any Voronoi diagram for the set of $n\geq 3$ points (which do not lie on one line) has at most $2n-5$ vertices and at most $3n-6$ edges.
Proof
Suppose that among the given points there are no three different that lie on one line. Denote $m$ and $h$ the number of vertices and edges of the
V-diagram, respectively. If we add a vertex $v_{\infty}$ "in infinity" which is an intersection of all the edges of the V-diagram that are halflines, we create a planar graph with $n$ faces, $m+1$ vertices and $h$ edges. See the figure.
Figure 9.2 The planar graph made of the Voronoi diagram.
The graph satisfies the Euler Formula
$$m+1-h+n=2.$$
The degree of every vertex is at least 3, the sum of all degrees is $2h$. Hence
$$2h\geq 3(m+1).$$
Substituting into this inequality $h=m+n-1$ from the Euler Formula we get
$$2m+2n-2\geq 3(m+1),$$
which gives $m\leq 2n-5$. Similarly after substituting $m=h-n+1$ we obtain
$$h\leq 3n-6.$$
\(\Box\)
For a finite set of points $P$ and any point $q$ in the plane not lying in $P$ we define a disc $C_P (q)$ with the center $q$ and the radius equal to the distance of the point $q$ from the set $P$
where $\d(q, P)$ is the distance from $q$ to its closest point from $P.$
Note that the circle bounding this disc contains at least one point of the set $P$.
Obviously we have:
Lemma 9.2
A point $q$ lies on an edge of the V-diagram if and only if there are at least two points of the set $P$ lying on the boundary of $C_P(q)$.
A point $r$ is a vertex of the V-diagram if and only if there are at least three points of the set $P$ lying on the boundary of $C_P(r)$.
Figure 9.3 Illustration to Lemma 9.2.
Sweep line method
One of the possible algorithms for the construction of the V-diagram uses the sweep line method, but in a way that is somewhat more sophisticated than those with which we have met so far. The sweep line $l$ moves from the top to the bottom and creates the V-diagram above it. This diagram does not occur in the entire halfplane $l^+$ above the line $l$, but only in a part of it, which is the union of areas above certain parabolas. How does this area look like?
A point $q\in l^+$ is lying in this area if its distance from a point $p_i\in P\cap l^+$ is less or equal to its distance from the line $l$
Then for all points $p_j\in P$ from the halfplane $l^-$ under the line $l$
$$\d(q,p_j)\geq\d(q,l)\geq\d(q,p_i).$$
So the point $q$ from the area determined by the inequality (1) does not lie
in the Voronoi cell $V(p_j)$ of the point $p_j\in P\cap l^-$.
Denote $\alpha(p,l)$ the parabola formed by the points which have the same distance from the point $p$ as from the line $l$. The point $p$ is the focus of the parabola and
the line $l$ is called the directrix of the parabola. Let $\alpha^+(p,l)$ is the set of points lying above or on the parabola $\alpha(p,l)$. The area above the sweep line $l$ where the V-diagram will be formed, is the union
$$\bigcup_{p_i\in P\cap l^+}\alpha^+(p_i,l).$$
The boundary of this area, formed by arcs of parabolas, is called the beach line.
In this case the tree $\mathcal T$ connected with this sweep line method will
be a balanced binary tree that maintains the order of arcs of the beach line in its leaves. Breakpoints of the beach line, i.e., the points where two consecutive arcs of the beach line with the focuses $p_1$ and $p_2$ meet, have the same distance from $p_1$ and
$p_2$, as
$$\d(q,p_1)=\d(q,l)=\d(q,p_2),$$
and so they lie on edges of the Voronoi diagram. These breakpoints are internal nodes of the tree $\mathcal T$. Their description of the form $\langle p_1\,:\,p_2\rangle$
means that the node corresponds to the intersection of the parabolas with
the focuses $p_1$ and $p_2$ (both have the sweep line $l$ as the directrix) where the arc of the parabola with the focus $p_1$ is on the left and the arc of the parabola with focus $p_2$ is on the right. Using such a tree, for a point $p$ on the sweep line $l$,
you can find a beach line arc that is above the point $p$. The following figure shows an example of a beach line and a corresponding tree.
Figure 9.4 A beach line.
Figure 9.5 A balanced binary tree corresponding the the previous beach line.
Events
The answers to the following two questions are of fundamental importance:
In which position of the sweep line does a new arc appear on the beach line?
In which position of the sweep line does an arc disappear from the beach line?
The following two lemmas answer these questions. Their proofs are technical, so we will not do them, and rather limit ourselves to illustrating them with animations.
Lemma 9.3
New arcs in the beach line come into being just when the sweep line passes through some point of the set $P$.
Such a point is called site event. The following animation illustrates this situation. The sweep line $l$ goes through the point $p\in P$. Let $\alpha$ be the arc of the beach line above the point $p$ with the focus at the point $p_1\in P$. After the sweep line passes through the point $p$, the arc $\alpha$ splits into two arcs $\alpha_1$ and $\alpha_2$ and a new arc $\beta$, a part of the parabola with focus $p$, appears between them in the beach line.
The intersections of the arcs $\alpha_1$ and $\beta$ and of the arcs $\beta$ and
$\alpha_2$ are points of the edge of the Voronoi diagram which separates the
Voronoi cell $V(p_1)$ from $V(p)$. The change in the beach line corresponds to the following change of the tree $\mathcal T$:
Figure 9.7 The change of $\mathcal T$ corresponding to the previous site event.
The internal nodes $\langle p_1\, :\, p_2\rangle$, $\langle p_1\,:\,p\rangle$ a $\langle p\,:\,p_1\rangle$ of the tree $\mathcal T$ describe breakpoints of the beach line
using the focuses of the corresponding parabolas in the order from left to right, as we have explained above. After changing the tree $\mathcal T$ we have to rebalance it.
Now consider three consecutive arcs $\alpha$, $\beta$ and $\gamma$ of the beach line. Let $p_1$, $p_2$, and $p_3$ be gradually the focuses of their respective parabolas. Consider the circle curcumscribed to the triangle $p_1p_2p_3$. Let us call $s$ its center and denote $q$ the lowest point of this circle, that is, the point with the minimal coordinate $y$. If the point $q$ lies under the sweep line $l$, we call it a circle event for the arc $\beta$.
When the sweep line $l$ passes this point, the center $s$ of the curcumcircle of the triangle $p_1p_2p_3$ has the same distance from these points as from the sweep line
$l$, so it lies on all three arcs $\alpha$, $\beta$ and $\gamma$. When moving the sweep line $l$ under $q$, the arc $\beta$ disappears from the beach line. The center $s$ of the circle is a vertex of the Voronoi diagram in which the edges separating the cell $V(p_1)$ from the cell $V(p_2)$, $V(p_2)$ from $V(p_3)$ and $V(p_1)$ from
$V(p_3)$ meet.
Figure 9.9 The change of $\mathcal T$ corresponding to the previous circle event.
Lemma 9.4
An existing arc of the parabola disappears from the beach line if and only if the sweep line passes through its circle event.
Algorithm
The start of the algorithm corresponds to the position of the sweep line $l$ above all points. The tree $\mathcal T$ is empty, and the event queue $\mathcal Q$
contains all the points of the set $P$ (site events) arranged lexicographically from top to bottom and from left to right. When the sweep line passes through an event, depending on whether it is a site or circle event, we do the following actions:
We add an arc to the tree $\mathcal T$ and edges into the double-connected edge list for the V-diagram,
We remove an arc from the tree $\mathcal T$ and we add an edge and a vertex to the V-diagram,
From the queue $\mathcal Q$ we remove the event that we have just passed, eventually we also remove the circle events that correspond to the triples of consecutive arcs from which one has disappeared after passing through the event. The order of arcs changes, and the current circle events are counted from three consecutive arcs in the new order. New circle events are added to the queue.
We do this procedure until the queue $\mathcal Q$ is empty. However, the tree
$\mathcal T$ will not be empty. Its current internal nodes determine the edges of the V-diagram that are halflines. This allows us to find a rectangle $R$, in which all the points of the set $P$, all the vertices of the V-diagram and the break points corresponding to these remaining internal nodes of the tree $\mathcal T$ are located.
The course of the algorithm in the simplest possible situation of a set of three points is demonstrated by the following animation.
AlgorithmVoronoiDiagram($P$)
Input. A set $P:=\{p_1,\ldots, p_n\}$ of point sites in the plane.
Output. The Voronoi diagram Vor($P$) given inside a bounding box in a doubly-connected edge list $\mathcal D$.
Initialize the event queue $\mathcal Q$ with all site events, initialize an empty status structure $\mathcal T$ and an empty doubly-connected edge list $\mathcal D$.
while $\mathcal Q$ is not empty
do Remove the event with largest $y$-coordinate from $\mathcal Q$.
if the event is a site event, occurring at site $p_i$
thenHandleSiteEvent($p_i$)
elseHandleCircleEvent($\gamma$), where $\gamma$ is the leaf of $\mathcal T$ representing the arc that will disappear
The internal nodes still present in $\mathcal T$ correspond to the half-infinite edges of the Voronoi diagram. Compute a bounding box that contains all vertices of the Voronoi diagram in its interior, and attach the half-infinite edges to the bounding box by updating the doubly-connected edge list appropriately.
Traverse the half-edges of the doubly-connected edge list to add the cell records and the pointers to and from them.
HandleSiteEvent($p_i$)
If $\mathcal T$ is empty, insert $p_i$ into it (so that $\mathcal T$ consists of a single leaf storing $p_i$) and return. Otherwise, continue with steps 2–5.
Search in $\mathcal T$ for the arc $\alpha$ vertically above $p_i$. If the leaf representing $\alpha$ has a pointer to a circle event in $\mathcal Q$, then this circle event is a false alarm and it must be deleted from $\mathcal Q$.
Replace the leaf of $\mathcal T$ that represents $\alpha$ with a subtree having three leaves. The middle leaf stores the new site $p_i$ and the other two leaves store the site $p_j$ that was originally stored with $\alpha$. Store the tuples $\langle p_j, p_i \rangle$ and $\langle p_i, p_j \rangle$ representing the new breakpoints at the two new internal nodes. Perform rebalancing operations on $\mathcal T$ if necessary.
Create new half-edge records in the Voronoi diagram structure for the edge separating $\mathcal V (p_i)$ and $\mathcal V (p_j)$, which will be traced out by the two new breakpoints.
Check the triple of consecutive arcs where the new arc for $p_i$ is the left arc to see if the breakpoints converge. If so, insert the circle event into $\mathcal Q$ and add pointers between the node in $\mathcal T$ and the node in $\mathcal Q$. Do the same for the triple where the new arc is the right arc.
HandleCircleEvent($\gamma$)
Delete the leaf $\gamma$ that represents the disappearing arc $\alpha$ from $\mathcal T$. Update the tuples representing the breakpoints at the internal nodes. Perform rebalancing operations on $\mathcal T$ if necessary. Delete all circle events involving $\alpha$ from $\mathcal Q$; these can be found using the pointers from the predecessor and the successor of $\gamma$ in $\mathcal T$. (The circle event where $\alpha$ is the middle arc is currently being handled, andhas already been deleted from $\mathcal Q$).
Add the center of the circle causing the event as a vertex record to the doubly-connected edge list $\mathcal D$ storing the Voronoi diagram under construction. Create two half-edge records corresponding to the new breakpoint of the beach line. Set the pointers between them appropriately. Attach the three new records to the half-edge records that end at the vertex.
Check the new triple of consecutive arcs that has the former left neighbor of $\alpha$ as its middle arc to see if the two breakpoints of the triple converge. If so, insert the corresponding circle event into $\mathcal Q$ and set pointers between the new circle event in $\mathcal Q$ and the corresponding leaf of $\mathcal T$. Do the same for the triple where the former right neighbor is the middle arc.
Time complexity
Let the set $P$ for which we implement the algorithm have $n$ points. After the first site event in the queue, the beach line will be a single parabola. When passing through other site events, the number of arcs in the beach line increases at most by two. Therefore, the total number of arcs that sometimes occur in the beach line will not exceed $2n-1$. In a binary balanced tree with the maximum of $2n-1$ leaves, the operations of addition of two leaves, removing one leaf and rebalancing the tree
take the time
$$O(\log(2n-1))=O(\log n).$$
Since each arc in the beach line arises once and disappears at most once, the number of executed events will be $O(2n-1) = O(n)$. In each event, we execute a finite number of steps which last a constant time or time $O(\log n)$. This gives the proof of the following statement:
Theorem 4.1
The Voronoi diagram for $n$ points can be computed using the above sweep line algoritm in $O(n\log n)$ time, using $O(n)$ storage.
An implementation of the algorithm using Matlab has been done by Jana Tomšů in the bachelor thesis (in Czech) available at
