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 code - Hide 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.

Old comments (closed because of spam)

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 :P)

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. Andrés said:

Amazing explanation! I am using your C++ library and I don’t understand what width and height exactly is. I put a value of width and height 500 but the edges starts and ends go through this limit and I get some bad lines. If I go up to 10000 it gets better but still some wrong lines when drawing it (not in OpenGL).

May 29th, 2013
16. Maa said:

Beautiful code
it look like Fortune / Shane O’Sullivan
but well structured, symetrical and clean, Bravo,
Your page on fortune is the best also..

Somebody mentionned bad line I got the same troubles with Fortune / Shane O’Sullivan that I started with before finding this one (too bad).
first you have to remove identical site (they cause trouble),
then you have to move to double instead of float
infact you even have to remove the too close point because they end up being like identical point and generate roudn off which produce these lines.

I am temped to restart from scratch with this.
Thanks again
PS: I love the comments in you native language would be nice to have the meaning for us

June 29th, 2013
17. isharko said:

Hey Ivan, thanks for the code! I’m trying to implement it in Cinder, a c++ framework and noticed little glitches appear occasionally in the diagram, where edges are drawn across the entire surface at strange angles. This vine illustrates that. Any ideas? Thanks,

https://vine.co/v/hmeIxZrE23h

July 27th, 2013
18. Sebastian Ziebell said:

Hi Ivan,

first of all thanks for the excellent explanation of the Voronoi diagrams, and also for sharing the C++ code. The last couple of days I tried to get my head around the algorithm and got the sample application running based on the changes fern17 made.

As isharko mentioned above I had similar strange artifacts / lines that cross the whole window (see http://twitpic.com/d5742m) for a single frame. After debugging the code and having an initial situation, where the artifacts would appear (sites did not move) I figured out the problem. Inside the function “InsertParabola” in the condition “if(root->isLeaf && root->site->y – p->y < 1)" it needs to be a 0, not a 1. My guess is, there are site constellations that do not invoke the necessary circle events when this condition is in place. I hope this is correct, at least for me it fixes everything and I could not observe any strange behavior anymore.

Best regards and thx again

July 28th, 2013
19. Sebastian Ziebell said:

Seems I am wrong. After running the sample a few times with different numbers of sites these glitches still appear. :-/

July 28th, 2013
20. Ivan Kuckir said:

Dear friends,

I am very sorry for the poor quality of my implementations. I wrote them when I was 19 and I knew nothing about the programming. Please, read just the pseudocode and make your own implementation. I will try to review it in the future, but I can’t promise anything.

July 30th, 2013
21. Tom Johnson said:

Ivan,

Great write-up! And nice code!

I think the issue with the C++ code that people are encountering is probably due to how line intersections are handled in Voronoi::GetEdgeIntersection() and Voronoi::FinishEdge(). The code in these functions uses the linear function coefficients “f” and “g”, which correspond to the equation:

y = f*x + g

The problem is that this representation of a line is not defined if the line is vertical (and leads to numerical computation problems if the line is near-vertical.)

A simple example is if you were working with two sites: (1, 1) and (2, 1). This should give you a single vertical “edge” at x=1.5. But for this example, the values for “f” and “g” as calculated by the C++ code end up being “inf” and “-inf”, respectively. (At least on my computer.) And thus, subsequent computation depending on “f” and “g” is all garbage.

A better approach might be to use the vector representation to calculate the intersections. Or at least special-case vertical and near-vertical edges.

August 6th, 2013
22. Steve said:

Thank you Sebastian Ziebell. I have been looking for a fortune solver that worked well and didn’t consume too much memory, processing power, etc. I found Ivan Kuckir’s and it works great, but I was having the same problem you had. I changed the 1 to a 0 in the insertParabola function just like you suggested and everything started working perfectly. Thank you and Ivan.

August 16th, 2013
23. How to Use Voronoi Diagrams to Control AI | Gamedevtuts+ said:

[...] Author: Mikko JohanssonDetailed Fortune’s Algorithm implementation: http://blog.ivank.net/fortunes-algorithm-and-implementation.htmlIncremental Triangle InsertionAnother method is to incrementally insert one point at a time, [...]

October 8th, 2013
24. How to Use Voronoi Diagrams to Control AI - Game Development KB said:

[...] Detailed Fortune’s Algorithm implementation: http://blog.ivank.net/fortunes-algorithm-and-implementation.html [...]

October 9th, 2013
25. Jay said:

I am studying your code of c++ version. Thank you so much for your contribution.

I have a question, it may be stupid. But I am confused that, in your voronoi diagram, your points/sites are moving. But, where does those movement comes from? I did not see you generate the points again and again… Would you please tell me about it?

Thanks you.

November 6th, 2013
26. Jay said:

Hi Ivan,

Sorry about the stupid question in the last email. I just figure out that.

Thanks.

November 14th, 2013
27. Thales said:

Hi, thank you for the pseudocode. It is much easier to understand than anything else I could find on the web.

I do not understand one line in particular though.

par = arc under point u;
if (par has its circle event, when it is removed form the beachline)
remove this event form the queue

What do you mean by that if statement?

December 25th, 2013
28. dunk said:

hi Ivan, thanks for the code! it’s the only sample i’ve found online that returns a graph of both the initial sites and the circle events.

i too saw the same artifacts that isharko mentioned.

Maa actually told us all the solution one post before: “first you have to remove identical site (they cause trouble),”

i can report that duplicate coordinates in the input vor::Vertices structure cause problems.
the following modified snippet from around line 36 in Voronoi.cpp will check for duplicates and skip them:
———————————————————————
VEvent * e;
VPoint checkDup(-1,-1);
while(!queue.empty())
{
e = queue.top();
queue.pop();
// Duplicate entries in queue are not handled correctly so check for them here.
if(e->point->x == checkDup.x and e->point->y == checkDup.y) continue;
checkDup.x = e->point->x; checkDup.y = e->point->y;
ly = e->point->y;
if(deleted.find(e) != deleted.end()) { delete(e); deleted.erase(e); continue;}
if(e->pe) InsertParabola(e->point);
else RemoveParabola(e);
delete(e);
}
———————————————————————

thanks again!
dunk.

January 8th, 2014
29. dunk said:

ok, that fixed some issues but not all.
Tom Johnson appears to be correct in that the way the VEdge construction calculates f leads to a loss of precision as the denominator approaches 0 in the following:
f = (b->x – a->x) / (a->y – b->y) ;

will post a solution when i have one.

January 8th, 2014
30. dunk said:

curious…
so i tried using an alternative way of calculating the intersection point of line segments that would work for lines parallel to the x axis.

things improved but there is still occasional corruption.

i wonder if it is happening when a 4th point happens to be in the circle during a circle event…?

meh.
i had some fun with this but must move on.
i’ve found a more robust implementation of this algorithm in the “voronoi” folder of this project:
http://www.sourceforge.net/projects/mapmanager

have fun!
dunk.

January 11th, 2014
31. dunk said:

the Voronoi library i ended up using for my project is this one:
http://www.boost.org/doc/libs/1_53_0_beta1/libs/polygon/doc/voronoi_main.htm .
after trying several it was the only one i found that did not lead to occasional graph corruption.

March 11th, 2014