## Introduction

Voronoi diagrams are quite useful tool in computational geometry. It has have been used, for instance, to
calculate the area per tree in the forest, to figure out where the poisoned wells where in a city based on the address of victims and so on. Not to go too
much into detail, it generally appears were you want to find out who is closest to whom. A collection of problems were Voronoi diagrams are used, is given below:

- Collision detection
- Pattern recognition
- Geographical optimization
- Geometric clustering
- Closest pairs algorithms
- k-nearest-neighbor queries

If
you are interested in even more examples you can go to the website here:

The history behind the Voronoi
diagram is quite complicated as it has been rediscovered several times. The
first written description of something that looks like a Voronoi diagram was
done by Descartes in the middle if the 17th century. The first to consider then
general idea was Dirichlet in 1850, but was not given a rigid mathematical
treatment before Voronoi's article in 1908 (for a more complete history you should read the article by Aurenhammer). As an important side note, Delaunay
(1934) proved that the Voronoi diagram had a sister (or brother if you will),
called Delaunay triangulation, which usually can be constructed from a Voronoi
diagram and vice versa.

In the following decades
after the initial discovery there were several approaches for constructing a
Voronoi diagram, but most of them were usually quite slow O(n^2). There were an
algorithm that could do it in O(n*log n) time, but it were quite difficult to
implement, so it was barely used. It was not until Fortune discovered the sweep
line algorithm in 1985 (published in 1987), which also runs on O(n*log n) time,
that a generally accepted outline for constructing a Voronoi diagram was
implemented. It is this algorithm I'm going to try to explain in detail in this
article. However a brief overview of the different algorithms for constructing
the Voronoi diagram is given below:

- Sweep line
- Divide-and-conquer
- Spiral-search
- Three dimensional convex hull
- Incremental

The
reason for writing this article about Fortune's sweep line algorithm is to
clarify the implementation of a Voronoi diagram. It is mostly done to clarify
things for me, as I try to implement the algorithm, with the hope that it is
interesting for others as well. This article would hopefully increase your
understanding of the implementation. And yes, I do know the implementation done
by Ben in
C#, also available on the Codeproject site. Some of the code used in this article even came from his
implementation.

## Background

What is a Voronoi diagram anyway? Well, there are several different types of Voroni diagram, but they all share some common attributes that are quite general in style and can be explained by using simple metaphors. One can easily describe a Voronoi diagram as the area around a point, where the center point is the closest one in the surrounding polygon, as illustrated below:

For a simple case of just two point, we can find the line that is equally far from the two points easily. It is done by constructing two circles that has a radius that is greater than half the distance between the two points. An illustration of the halfway line between the point A and B is given below:

Another way of looking at it is by assuming that all the points, are wooden sticks were all have the same height. If you throw a perfect blanket over all the sticks, assuming that the sticks a high enough so the blanket doesn’t touch the ground, the valleys would give you all the lines in the Voronoi diagram. The reason can be "proven" with the illustration of a rope tied on two equally high poles. The halfway point would be found were the rope is at its lowest point, and thus is a complementary way of finding the halway distance as the illustration of the two circles above.

Hopefully the image is simple for everyone to grasp, but it is not convenient to construct the diagram from this analogy, or in Steven Fortune words: "It is notoriously difficult to obtain a practical implementation of an abstractly described geometric algorithm”.

A mathematically way of looking at the diagram is to place cones (or a pointy hat if you’d like, see picture above) with slopes at the side in a 45 degree angle with the center point in each point. This would give overlapping sides that give you the straight lines between the points in the Voronoi diagram. This would in fact mean that you will always find the overlapping line exactly at the half way point between any of the closest points. The result from the intersection lines from the cones above is the Voronoi diagram below:

Fortunes sweep line algorithm is based on these principles, but instead of placing the cones with the point at the top, as the cones are turned upside down. The sides of the cone tilted in 45 degrees and an intersection sweep line (the plane is named pi in the picture below) is tilted in 45 degrees opposite to the tilted cones. All angles are given relative to the plain at z = 0. This approach is illustrated below were the sweep line constantly moves from 0 to x = infinite (or rather the boundary rectangle that is certain to contain all the points):

It is quite clear for anyone that has had the explanation of the different geometrical curves (see the article about Conic slides from Wikipedia) that the intersection line between the sweep plain and the cone would form an parabolic equation. This intersection line is dependent on the relative point at the center of the cone P(x,y) and the position of the sweep line L. We will call the line that intersection line created by the sweep-plain pi and the cones, the beach line from now on.

The algorithm must maintain the current location (y-coordinate) of the sweep line. It stores, in left-to-right (or top to bottom, as you can perform sweeps in both X and Y direction) order the set of sites that define the beach line. It is however important to know that the algorithm never needs to store all the parabolic arcs of the beach line. We will maintain a sweep-line status
and we are interested in simulating the discrete event points where
there is a “significant event”, that is, any event that changes the
topological structure of the Voronoi diagram and the beach line.

Here the major problem in the construction of data structures is met. Namely: "When and how the topological structure of parabolic front is changing?" The answer is hidden in events. Two kinds of events appear. The first is a point event and the second a circle event.

### Site events

When the sweep line passes over a new site a new arc will be inserted into the beach line.

This looks like a complicated thing to calculate, but given that a parabolic equation for a given point Pj and the y coordinate for the sweep line, you can calculate the distance from the new point on the sweep line (the new point gives you the x and y coordinate and lies on the sweep line, and you could plug the x value from the new point into the formula).

The equation for any parabola given Fortune's sweep line algorithm is given by the equation below were Pjx and Pjy is x and y position to the center point of the cone above the sweep line, Ly is the position of the sweep line, and X is the location of the new intersection point on the sweep line:

The formula must be used on all the points passed by the sweep line in order to find the one that gives you the shortest line from the point to the sweep line. You can now easily generate all the arcs generated from all the points by looping though all the points as the sweep line descend.

However the first intersection point is not that interesting when you are trying to compute the Voronoi diagram. What is interesting is to find the two points in the picture to the left, where the beach line form two intersection points that would form the beginning of a new line in the Voronoi diagram. So we assume that we have to points that form two parabolas, two intersection points and each individual parabola would be governed by the parabolic equation above. We can now set the two parabolic equations equal each other and calculate the values were the two intersections occur by the quadratic formula. The result is two x values that you could plug into any of the two parabolic equations to get the resulting y value, and thereby get the two cross points. The equation is given below were Pi and Pj is the two center site points of the cones:

I do not want to clean up the equation, as it is quite a long, but the idea is to make it into an equation that could be written on the form:

With the two solutions for X being found from the equation below:

But within the equation lies a pretty big problem, because of the site events when are you suppose to make the guess about the line in a Voronoi diagram? If you increment it a little before calculation the line, you will get a large error, and if you wait too long you get a mess as you have to shorten the line between two points, as new intersections would break the calculated lines.

You might be tempted to ask for how long you would have to keep checking the parabolic line for any given point, and the answer is that you have to keep checking the parabola of a point as long as the center point is not enclosed in lines, meaning that no event further down the sweep line would change the boundaries for any of the lines surrounding the point. This is normally not a problem as the checks are event driven by two conditions, and a parabolic arc is terminated when a circle event or sometimes called the vertex event occurs. The length of a parabolic arc (at the beach line) shrinks to zero, the arc disappears and a new Voronoi vertex will be created at this point.

### Circle event

This event could only occur if two parabolas swallow a third parabola, meaning that we only need to check if three (or more) points will form a circle that does not contain any other point.

Witch arc would dissapear you might ask. The answer is that, if you have the three points that form a circle like the picture above, it is the arc beloning to the point with the highest Y value (of the points that make up the circle of course) that would dissapear.

There is a formula for getting the center and the radius of a circle here. If want to know if another point is inside the circle you could simply calculate the distance between the center point of the circle and the potential point, and if the distance is less than the radius the point is inside the circle. If not the center of the circle would give you another vertex point and one parabola would disappear. The calculation is done by using the general equation of a circle:

(X - Xc)*2 + (Y - Yc)*2 = R*2

And substitute the three points for X and Y, you'll have three equations with three unknown variables which then can be solved for Xc, Yc and R.

There is also a third possible "event", which is not discussed so far. This would happen when the intersection will occur outside the predefined size of the rectangle, were an infantry long line would appear. This could be fixed when you are done with the sweep line algorithm, were you could check if some of the lines are too short and have not had any events occur to the parabola.

## Elementary functions used in the code

A heads up first, this code uses the native library in .NET which may have some limited uses for an
arbitrary collection of points. This is done so the celerity of the code is not compromised.

The first thing we need to do with an arbitrary collection of points is to sort them in the sweep line direction. A general algorithem sorting points in either x or y direction,
with lexographically ordered points, meaning that if you sort by y value, and get equal y's, you would sort the two points by the x direction after the y check, is given below:

Private Function SortPoints(ByVal samplepoints As PointCollection, _
ByVal SortByXdirection As Boolean) As PointCollection
Dim copySamplePoints As Point() = New Point(samplepoints.Count - 1) {}
samplepoints.CopyTo(copySamplePoints, 0)
If SortByXdirection Then
Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.X))
Else
Array.Sort(copySamplePoints, New PointSort(PointSort.Mode.Y))
End If
Dim result As New PointCollection
For Each p As Point In copySamplePoints
result.Add(p)
Next
Return result
End Function
Private Class PointSort
Implements IComparer
Public Enum Mode
X
Y
End Enum
Private currentMode As Mode = Mode.X
Public Sub New(ByVal mode As Mode)
currentMode = mode
End Sub
Private Function IComparer_Compare(ByVal a As Object, ByVal b As Object) _
As Integer Implements IComparer.Compare
Dim point1 As Point = CType(a, Point)
Dim point2 As Point = CType(b, Point)
If currentMode = Mode.X Then
If point1.X > point2.X Then
Return 1
ElseIf point1.X < point2.X Then
Return -1
Else
If point1.Y > point2.Y Then
Return 1
ElseIf point1.Y < point2.Y Then
Return -1
Else
Return 0
End If
End If
Else
If point1.Y > point2.Y Then
Return 1
ElseIf point1.Y < point2.Y Then
Return -1
Else
If point1.X > point2.X Then
Return 1
ElseIf point1.X < point2.X Then
Return -1
Else
Return 0
End If
End If
End If
End Function
End Class

A function for finding the intersection parabola of a new site event:

Private Function Intersection(ByVal NewPoint As Point, ByVal OldPoint As Point, _
ByVal SweepLine As Double) As Polyline
Dim result As New Polyline
result.Points.Add(NewPoint)
Dim endpoint As New Point
endpoint = NewPoint
Dim res As Double
res = (1 / (2 * ((OldPoint.Y) - SweepLine)))
res *= (NewPoint.X ^ 2 - 2 * OldPoint.X * NewPoint.X + OldPoint.X ^ 2 + OldPoint.Y ^ 2 - SweepLine ^ 2)
endpoint.Y = res
result.Points.Add(endpoint)
result.Stroke = Brushes.Red
result.StrokeThickness = 2
Return result
End Function

This next function is used just to calculate the same as the function above, but this just returns the y value in the arc created by the oldpoint.

Private Function Intersection2(ByVal NewPoint As Point, _
ByVal OldPoint As Point, ByVal SweepLine As Double) As Double
Dim res As Double
res = (1 / (2 * ((OldPoint.Y) - SweepLine)))
res *= (NewPoint.X ^ 2 - 2 * OldPoint.X * NewPoint.X + OldPoint.X ^ 2 + OldPoint.Y ^ 2 - SweepLine ^ 2)
Return (res)
End Function

The next function we need to handle is the so-called site event, or as I called it a Parabolic cut (the sweep line is set at the location ys + 500 in the code below):

Private Function ParabolicCut(ByVal Point1 As Point, ByVal point2 As Point, ys As Double) As Polyline
Dim x1, y1, x2, y2 As Double
x1 = Point1.X
y1 = Point1.Y
x2 = point2.X
y2 = point2.Y
ys = ys + 500
Dim a1 As Double = 1 / (2 * (y1 - ys))
Dim a2 As Double = 1 / (2 * (y2 - ys))
Dim xs1 As Double = 0.5 / (2 * a1 - 2 * a2) * (4 * a1 * x1 - 4 * a2 * x2 + 2 * _
Math.Sqrt(-8 * a1 * x1 * a2 * x2 - 2 * a1 * y1 + 2 * a1 * y2 + 4 * a1 * a2 _
* x2 * x2 + 2 * a2 * y1 + 4 * a2 * a1 * x1 * x1 - 2 * a2 * y2))
Dim xs2 As Double = 0.5 / (2 * a1 - 2 * a2) * (4 * a1 * x1 - 4 * a2 * x2 - 2 * _
Math.Sqrt(-8 * a1 * x1 * a2 * x2 - 2 * a1 * y1 + 2 * a1 * y2 + 4 * a1 * _
a2 * x2 * x2 + 2 * a2 * y1 + 4 * a2 * a1 * x1 * x1 - 2 * a2 * y2))
Dim p1, p2 As New Point
p1.X = xs1
p2.X = xs2
p1.Y = 0
p2.Y = 0
p1.Y = Intersection2(p1, Point1, ys)
p2.Y = Intersection2(p2, Point1, ys)
Dim result As New Polyline
result.Points.Add(p1)
result.Points.Add(p2)
result.Stroke = Brushes.Black
result.StrokeThickness = 2
Return result
End Function

Circle event with a check to see if another point is inside the circle we create an own class called Circle:

Public Class Circle
Private _CenterPoint As Point
Public Property CenterPoint() As Point
Get
Return _CenterPoint
End Get
Set(ByVal value As Point)
_CenterPoint = value
End Set
End Property
Private _radius As Double
Public Property Radius() As Double
Get
Return _radius
End Get
Set(ByVal value As Double)
_radius = value
End Set
End Property
Private Function Distance(ByVal p1 As Point, ByVal p2 As Point) As Double
Return Math.Sqrt((p1.X - p2.X) ^ 2 + (p1.Y - p2.Y) ^ 2)
End Function
Public Function DoesCircleContainPoint(ByVal p1 As Point) As Boolean
Dim d As Double
d = Distance(p1, CenterPoint)
If d < Radius Then
Return True
Else
Return False
End If
End Function
Public Sub ReturnCenterPoint(ByVal can As Canvas)
Dim p As New Ellipse
p.Fill = Brushes.Blue
p.StrokeThickness = 2
p.Stroke = Brushes.Blue
p.Width = 5
p.Height = 5
can.Children.Add(p)
Canvas.SetLeft(p, CenterPoint.X - 2.5)
Canvas.SetTop(p, CenterPoint.Y - 2.5)
End Sub
Public Sub ReturnEllipse(ByVal can As Canvas)
Dim p As New Ellipse
p.Fill = Brushes.Yellow
p.StrokeThickness = 2
p.Stroke = Brushes.Black
p.Width = Radius * 2
p.Height = Radius * 2
can.Children.Add(p)
Canvas.SetLeft(p, CenterPoint.X - Radius)
Canvas.SetTop(p, CenterPoint.Y - Radius)
End Sub
Public Sub CreatCircleFromThreePoints(ByVal Point1 As Point, _
ByVal Point2 As Point, ByVal Point3 As Point)
Dim Pt(3) As Point
Pt(0) = Point1
Pt(1) = Point2
Pt(2) = Point3
GetCircleRectFromPoints(Pt)
End Sub
Private Sub GetCircleRectFromPoints(ByVal Pts() As Point)
Dim N As Integer = Pts.Length - 1
Dim X(N, N) As Double, Y(N) As Double
For I As Integer = 0 To N
X(I, 0) = 1 : X(I, 1) = Pts(I).X : X(I, 2) = Pts(I).Y
Y(I) = X(I, 1) * X(I, 1) + X(I, 2) * X(I, 2)
Next
Dim MatInv As New Elimination
MatInv.ComputeCoefficents(X, Y)
Dim A As Single = CSng(Y(1) / 2)
Dim B As Single = CSng(Y(2) / 2)
Dim R As Single = CSng(Math.Sqrt(Y(0) + A * A + B * B))
CenterPoint = New Point(A, B)
Radius = R
End Sub
Public Class Elimination
Sub ComputeCoefficents(ByVal X(,) As Double, ByVal Y() As Double)
Dim I, J, K, K1, N As Integer
N = Y.Length
For K = 0 To N - 1
K1 = K + 1
For I = K To N - 1
If X(I, K) <> 0 Then
For J = K1 To N - 1
X(I, J) /= X(I, K)
Next J
Y(I) /= X(I, K)
End If
Next I
For I = K1 To N - 1
If X(I, K) <> 0 Then
For J = K1 To N - 1
X(I, J) -= X(K, J)
Next J
Y(I) -= Y(K)
End If
Next I
Next K
For I = N - 2 To 0 Step -1
For J = N - 1 To I + 1 Step -1
Y(I) -= X(I, J) * Y(J)
Next J
Next I
End Sub
End Class
End Class

## The main code

As I was writing the article I soon realized that it would become a monster if I included all the necessary steps for creating a complete Voronoi diagram.
So I stopped with a version that does not, bt any means, contain the creation of a complete Voronoi diagram, but rather give you a chance to understand the most
necessary functions for creating it. The final complete implementation and creation will have to happen in the next article.

## History

First release: 29.06.2012.

Second release: 01.07.2012, Minor changes in languange, added some references inside the article.

## Referances