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:
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 moredimensional spaces and for different metrics.
Voronoi diagram as a graph
We can look at the Voronoi diagram like it is a nonoriented 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_{n1}, 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_{i1}, 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 codeIn 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 queueAs the event queue, we can use a priority queue, where events are sorted by "y" coordinate. It should allow us to delete items.
Borderline sequenceWe can represent a borderline sequence as a doublylinked 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 halfarcs and a new arc is added between them. The "right leaf" becomes an inner node and gets 2 children. Left child represents a left halfarc of a previous arc. Right child is a subtree, which contains a newly added arc on the left and right halfarc 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:
 For sites lying in the line, it is evidently true.
 Let's suppose, that sites are not in the line.
 $V$: number of vertices in Voronoi diagram
 $E$: number of edges
 $N$: number of inner faces = number of sites
$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 20times faster, than my implementation in ActionScript 3.
Great writeup 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 reimplement it in Java? Do you object to me using it?
May 8th, 2012Hello,
May 8th, 2012thank you. This whole thing was an idea of Steven Fortune and I am just describing it. Of course you can use it.
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)
https://docs.google.com/open?id=0B6oy2F9QV_YGcy1ibUN2R0poWmM
Greetings
May 11th, 2012Ahoj, 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, 2012Hello,
Very interesting code. Is somebody convert it to c# ?
Regards
June 4th, 2012Seems 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
July 7th, 2012points = [];
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);
seb. 49, there is a good conversion of this algorythn to C# at stackoverflow.
July 10th, 2012Some guy has implemented it in C# in 2005 already.
Oh excuse me that is not stackoverflow but codeproject.
July 10th, 2012Here is a link:
http://www.codeproject.com/Articles/11275/FortunesVoronoialgorithmimplementedinC
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, 2012Trying to port this to Unity (C#) lets see how it goes.
August 16th, 2012Thanks 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, 2012Hi, 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,
January 12th, 2013Rus
Hello,
beachline is an array of [... arc, edge, arc, edge, arc, ...]
January 12th, 2013By “replace par” I mean replacint the arc by subarray inside a beachline. You have to replace current arc by 3 arcs + 2 edges between them.
Oh! Suddenly everything makes sense! Thank you so much! (!)
January 12th, 2013Amazing 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, 2013Beautiful 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.
June 29th, 2013Thanks again
PS: I love the comments in you native language would be nice to have the meaning for us
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, 2013Hi 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, 2013Seems I am wrong. After running the sample a few times with different numbers of sites these glitches still appear. :/
July 28th, 2013Dear 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, 2013Ivan,
Great writeup! 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 nearvertical.)
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 specialcase vertical and nearvertical edges.
August 6th, 2013Thank 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[...] Author: Mikko JohanssonDetailed Fortune’s Algorithm implementation: http://blog.ivank.net/fortunesalgorithmandimplementation.htmlIncremental Triangle InsertionAnother method is to incrementally insert one point at a time, [...]
October 8th, 2013[...] Detailed Fortune’s Algorithm implementation: http://blog.ivank.net/fortunesalgorithmandimplementation.html [...]
October 9th, 2013I 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, 2013Hi Ivan,
Sorry about the stupid question in the last email. I just figure out that.
Thanks.
November 14th, 2013Hi, 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.
In your ‘AddParabola’ function:
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, 2013hi 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 wish now i’d read the comments before troubleshooting.
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!
January 8th, 2014dunk.
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, 2014curious…
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!
January 11th, 2014dunk.
the Voronoi library i ended up using for my project is this one:
March 11th, 2014http://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.