# 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:

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 codeHide the code

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.

1. Joe said:

Great write-up on this. Thanks for making it available. I’ve read a lot of Voronoi sites and papers recently and yours is easier to follow than most of them.

Is there a license for the source code? Can I re-implement it in Java? Do you object to me using it?

May 8th, 2012
2. Ivan Kuckir said:

Hello,
thank you. This whole thing was an idea of Steven Fortune and I am just describing it. Of course you can use it.

May 8th, 2012
3. fern17 said:

Hi! Thank you for this, it was really helpful!

One note: I tried to compile the C++ implementation on GNU/Linux with g++ and failed. I edited a few files of your project in order to compile it..
Here are the changes if you want to share it (and make it multiplatform )

Greetings

May 11th, 2012
4. Peter said:

Ahoj, som doktorand na stavebnej fakulte SvF STU BA, ak som to spravne pochopil, Ty si z MFF UK ? Dostal sa mi do ruk podklad ohladom Voronoiovych diagramov a z hladiska mojej dizertacnej prace vyzeraju byt pouzitelne ako vzor nosnej konstrukcie striech priestorovych skrupin. Nie som programator, takze pochopenie algoritmu je pre mna narocne. Bol by si ochotny zodpovedat mi zopar otazok ohladom tejto problematiky ?

Dakujem

May 15th, 2012
5. seb.49 said:

Hello,

Very interesting code. Is somebody convert it to c# ?

Regards

June 4th, 2012
6. Joakim said:

Seems there’s a collinear bug with the javascript version (at least).

Here’s 2 point sets that will trigger a bad delaunay triangulation.

// Set 1
points = [];
points.push(new Point(2, 2));
points.push(new Point(2, 4));
points.push(new Point(8, 4));
points.push(new Point(2, 6));
v = new Voronoi();
v.Compute(points, 10, 10);

// Set 2
points = [];
points.push(new Point(2, 2));
points.push(new Point(4, 2));
points.push(new Point(6, 2));
points.push(new Point(4, 8));
v = new Voronoi();
v.Compute(points, 10, 10);

July 7th, 2012
7. Freeliner said:

seb. 49, there is a good conversion of this algorythn to C# at stackoverflow.
Some guy has implemented it in C# in 2005 already.

July 10th, 2012
8. Freeliner said:

Oh excuse me that is not stackoverflow but codeproject.
http://www.codeproject.com/Articles/11275/Fortune-s-Voronoi-algorithm-implemented-in-C

July 10th, 2012
9. Sean Esopenko said:

Thanks for the easy to read source code. I’m working on understanding Vornoi diagrams and your code’s easy to understand and well laid out.

Thanks for sharing!

July 15th, 2012
10. Bjorn said:

Trying to port this to Unity (C#) lets see how it goes.

August 16th, 2012
11. Daniel Hunt said:

Thanks for this great writeup. I’m writing an implementation of Fortune’s algorithm in Java for a project and this has been invaluable in aiding my understanding of the algorithm. The analysis at the end was also very helpful in boosting my understanding of the bounds on performance – I haven’t found another site yet with that kind of analysis. Thanks!

October 15th, 2012
12. Rus said:

Hi, I’ve found your descriptions to be the single most descriptive descriptions on the web, however, there one section, of the psudocode, I just cannot understand.

replace par by the sequence a, xl, b, xr, c

Towards the end of the add parabola function. Is this simply saying, create the edges based off of the parabola?

Thanks,
Rus

January 12th, 2013
13. Ivan Kuckir said:

Hello,

beachline is an array of [... arc, edge, arc, edge, arc, ...]
By “replace par” I mean replacint the arc by sub-array inside a beachline. You have to replace current arc by 3 arcs + 2 edges between them.

January 12th, 2013
14. Rus said:

Oh! Suddenly everything makes sense! Thank you so much! (!)

January 12th, 2013
15. What To Consider When Designing A Certificate of Completion said:

How do I make the page look like a certificate on Microsoft Word?

May 24th, 2013