Fortune’s algorithm and implementation

This text was created as a credit work at subject Algorithms nad data structures II, MFF UK. It contains an implementation in C++ (at the end) and ActionScript 3.0 (in previous posts).

You can check out my implementations:

• Voronoi diagram in ActionScript 3
• Voronoi diagram in JavaScript
• Voronoi diagram in C++

About Voronoi diagram

Let us have a set $P$, consisting of $n$ points in a plane. Let us call these points sites. Now let’s make a function $V : P \rightarrow \mathcal{P}(plane)$, so it attaches all points, that are the closer to this site, than to any other site, i. e. $dist(x, s) \leq min_{z\in P}\left(dist(x,z)\right) \rightarrow x \in V(s)$. Let’s call these $V(s)$ cells. Each site has a cell of the closest points around it.

Note, that cells are not disjoint. There are some points in the plane, which have the same distance to 2 or more sites. We call these points a Voronoi diagram. Intersection of two cells is an empty set, line segment, ray or line.

Similarly, we can define a Voroni diagram for more-dimensional spaces and for different metrics.

Voronoi diagram as a graph

We can look at the Voronoi diagram like it is a non-oriented graph. Its vertices are intersections of 3 or more cells and edges are intersections of 2 cells. We can add 1 more “virtual” vertex and connect all the infinite edges with it, to get a standard graph.

Delaunay triangulation

Let us have a set of points $P$ in a plane. We can create many triangulations on them. Delaunay triangulation for a set $P$ of points in the plane is a triangulation $DT(P)$ such that no point of $P$ is inside the circumcircle of any triangle in $DT(P)$. It can be shown, that Delaunay triangulation is a dual graph to Voronoi diagram (as a graph) for any set $P$ of sites.

Algorithms for computing the Voronoi diagram

• Divide and conquer – algorithm divides the points into right and left part, it recursively computes a Voronoi diagram for these two parts and finally it “merges” these two parts by computing voronoi lines between them.
• Incremental algorithm – it counts a Voronoi diagram for two sites. Then it takes other sites, one by one, and edits current diagram.

These algorithms have many disadvantages – difficult implementation, nontrivial merging of diagrams, numerical inaccuracy.

Fortune’s algorithm

This algorithm is based on “sweeping” the plane. There is a line $L$ moving from bottom to the top, wich proceeds each point it passes. Let’s think, which points under $L$ already belong to Voronoi diagram.

If the point $x$ has the same distance to $a$ and $b$ (all 3 lies under a sweep line), this point $x$ doesn’t have to be in diagram, because there can be next site right above the sweep line, which is closer. We can decide a Voronoi diagram only for points, which are closer to any site under the line than to the line $L$. These points lie under a parabloa, site $a$ is its focus point and the sweep line is its directrix, due the definition of a parabola.

The border, which says, for which points we can draw a diagram, is made of an arcs of parabolas. We call it a beachline or borderline. Intersections of these parabola’s arcs have the same distance to some 2 sites (2 focii of parabolas) and to the sweep line. Thus, when moving a sweepline upwards, these intersections draw the diagram. The beachline is a sequence of arcs $p$ and intersections $x$, which we can call “edges”, because they draw the edge of the Voronoi graph.

Beachline $ p_0, x_0, p_1, x_1, \ldots , x_{n-1}, p_n $ .





The beachline changes while the sweepline moves upwards. Beachline can be changed only at 2 occasions:

• Site event – the sweepline passes a new site, which adds a new parabola into the beachline. The beachline is divided into 2 arcs and 2 edges (2 parts of final edge) start to come in the opposite directions under this new site. Thus, the parabola $p_i$ is replaced by a sequence $p_i, x_i, p_{i+1}, x_{i+1}, p_{i+2}$, where $p_i$ a $p_{i+2}$ have a focus of a previous parabola and $p_{i+1}$ has its focus in a newly added site.
• Circle event – an arc disappears, two neighbour arcs “squeeze” it. A new edge between neighbours is started. Thus, a sequence $x_{i-1}, p_i, x_i$ is replaced by a new edge $x_i$. A place, at which arc disappears and new edge begins, is in the same distance to all 3 focii (its focus and neighbours’ focii). Thus, it is the circumcenter of the circumcircle, defined by these 3 focii – sites.

We will store these events in a sorted queue or heap, where they are sorted by “y” coordinate. We see, that the beachline changes only at a few moments (events), so we can move it discretely from one evnt to another.

Code of an algorithm

Input is set of points (sites). Output is set of edges or polygons.

Show the code – Hide the code

for each site
   {
      create a site event e, 
      e.point = current site, insert e into queue
   }
   while queue is not empty
   {
      e = get the first event from the queue
      if its site event : AddParabola( e.point )
      else : RemoveParabola( e.parabola );
   }
   // done!!!  

   function AddParabola ( point u )
   {
      par = arc under point u;
      if (par has its circle event, when it is removed form the beachline)
         remove this event form the queue
      new arcs a, b, c;
      b.site = u;
      a.site = c.site = par.site; // site of arc is a focus of arc
      xl, xr  = left and right edges, which comes from point on par under u
      xl is a normal to  (a.site, b.site);
      xr is a normal to (b.site, c.site);
      replace par by the sequence a, xl, b, xr, c
      CheckCircleEvent(a);
      CheckCircleEvent(c);
   }
   function RemoveParabola ( Parabola p )
   {
      l = an arc lef to p;
      r = an arc on the right from p;
      if (l or r have their Circle events) remove these events from the queue
      s = the circumcenter between l.site, p.site and r.site
      x = new edge, starts at s, normal to (l.site, r.site)
      finish two neighbour edges xl, xr at point s
      replace a sequence xl, p, xr by new edge x
      CheckCircleEvent(l);
      CheckCircleEvent(r);
   }
   function CheckCircleEvent(Parabola p)
   {
      l = arc on the left to p;
      r = arc on the right to p;
      xl, xr = edges by the p
      when there is no l  OR no r  OR  l.site=r.site  RETURN
      s = middle point, where xl and xr cross each other
      when there is no s (edges go like\ /) RETURN
      r = distance between s an p.site (radius of curcumcircle)
      if s.y + r is still under the sweepline  RETURN
      e = new circle event
      e.parabola = p;
      e.y = s.y + r;
      add e into queue 
   }

In many places, smart people say that it is enough just to check, if a top point of the circle is above the sweep line. But these two edges can go like \/ (or rotated like < and never meet each other).

Implementation

Event queue

As the event queue, we can use a priority queue, where events are sorted by “y” coordinate. It should allow us to delete items.

Borderline sequence

We can represent a borderline sequence as a doubly-linked list. Inserting and deleting items is easy, but searching an item takes a linear time. Steven Fortune came up with a great idea, to represent this sequence as a binary tree.

• Each leaf in the tree represents an arc. Each inner node is an edge between two arcs.
• If there is only one arc, the tree is only one leaf, root

• When we add a new arc into the tree (when a sweep line rolls over the new site), the “right leaf” (the arc under the site being added) splits into two half-arcs and a new arc is added between them. The “right leaf” becomes an inner node and gets 2 children. Left child represents a left half-arc of a previous arc. Right child is a subtree, which contains a newly added arc on the left and right half-arc of previous arc on the right.

  
• Removing an arc is also very easy. Its parent represents the left or the right edge, some “higher predecesor” is the next edge. If currently removed arc is the left child, we replace the predecesor by the right child and vice versa. We attach the new edge to the “higher predecesor”.

Complexity

Theorem: For $n \geq 3$ sites, Voronoi diagram contains at most $2n – 5$ vertices and $3n – 6$ edges.
Proof:

1. For sites lying in the line, it is evidently true.
2. Let’s suppose, that sites are not in the line.

Let’s denote these variables:

• $V$: number of vertices in Voronoi diagram
• $E$: number of edges
• $N$: number of inner faces = number of sites

From Euler’s formula we know, that:
$V – E + N = 2$
Because Voronoi diagram contains infinite edges, let’s create a new vertex “infinity” and connect all these edges to it. Now it is a planar graph.
$(V + 1) – E + N = 2$
In Voronoi graph we know, that every vertex has a degree at least 3 (including “infinity” vertex). An edge si between two vertices, so
$3 * (V + 1) \leq 2 * E$
By a simple algebraic modification, we get a previous theorem.

Primitive operations, such as searching for an item in the tree or in the queue, removing from the queue, all it can be done in $O(log(n))$. At each event, we do a constant number of these primitive operations.

The number of site events is $N$, circle events are at most $2N – 5$. At each of them, we do $c$ primitive operations. Thus, final complexity is $O(n*log(n))$.

Implementation in C++

HERE is my implementation, written in C++ and commented in english. It includes EXE file, which runs an animation of vertices and draws the diagram (using OpenGL). The speed is breathtaking – program can compute a diagram for 1000 vertices 50 times in a second. It is almost 20-times faster, than my implementation in ActionScript 3.