Introduction
Voronoi diagrams are quite useful tools in computational geometry and have a wide range of uses such as, calculating the area per tree in the forest, or figuring out where the poisoned wells were in a city (based on victims' addresses), and so on. In general it is useful for finding
"who is closest to whom." A collection of problems where Voronoi diagrams are used is shown 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 a Voronoi diagram was
done by Descartes in the middle if the 17th century. The general idea was first conceived by Dirichlet in 1850, but was not given a rigid mathematical
treatment before Voronoi's article in 1908. For a complete history behind it 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. (Fun fact: Delaunay actually did his Ph.d. thesis under the supervision of Voronoi.)
In the decades following
the initial discovery there were several approaches for constructing a
Voronoi diagram, but most of them were usually quite slow O(n^2). There was an
algorithm that could do it in O(n*log n) time, but it was 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 method for constructing a Voronoi diagram was
implemented. In this article I am going to try explain the basic functions for this algorithm. However, a brief overview of different algorithms for constructing
the Voronoi diagram are shown 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. This is mostly to clarify
things for me as I try to implement the algorithm, with the hope that it is
interesting for others as well and will hopefully increase your
understanding of the basics. And yes, I am familiar with Ben's article on CodeProject. Some of the code used in this article came from his
implementation.
Background
What is a Voronoi diagram anyway? Well, there are several different types of Voroni diagrams, but they all share some common attributes that are quite general in style and can be explained by using simple metaphors. It is normally described as the area around a point, where the center point is the closest one in the surrounding polygon, as illustrated below:
There is a general underlying pattern, that can be illustrated by finding the halfwayline between just two points. We find the line that is equally far from the two points 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 two points A and B:
If we were just interested in finding the halfway point on the line between points A and B, there is one more possible way to calculate the point. If you take a look at the illustration below where a rope is tied to the points A and B, you immediately see that the halfway point is where the rope is at its lowest point. Therefore, this a complementary way of finding the halway distance between the two points as in the illustration of the two circles above.
If we replace the rope by a blanket and set up wooden sticks that all have the same height, you could imagine that all the valleys would give you all the Voronoi lines in the diagram, assuming that the sticks are high enough (or the tension in the blanket is high enough) so
the blanket doesn’t touch the ground.
Hopefully the aforementioned image is simple enough for everyone to grasp, but it is not convenient to construct the Voronoi diagram from this analogy, or in Steven Fortune's words: "It is notoriously difficult to obtain a practical implementation of an abstractly described geometric algorithm."
A more mathematically sound way of looking at the Voronoi diagram is to place cones (or a
pointy hat if you’d like - see picture below) with the slopes at the side tilted in a
45 degree angle, with the center point as the highest point. This would give
the same overlapping sides that generates the straight lines between the points in
the Voronoi diagram as seen in the first picture.
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 cone is a collection of many circles with different radiuses so the analogy to the halfway line is still valid. The result from the example picture above gives the intersection lines below:
Fortune's sweep line algorithm is based on these principles, but instead of placing the cones with the point at the top, the cones are turned upside down. The sides of the cone are tilted at 45 degrees and an intersection sweep line (the plane is named pi in the picture below) is tilted at 45 degrees opposite to the tilted cones. You will also, I hope, notice that there is no height limitation on the cones by turning them upside down. All angles are given relative to the plain at z = 0. This approach is illustrated below where the sweep line constantly moves from 0 to x = infinite (or rather the boundary rectangle that is certain to complete all the lines in the Voronoi diagram for the area that we are interested in).
For anyone familiar with different geometrical curves (see the article about Conic slides from Wikipedia) it is clear that the intersection line between the sweep plain pi and the cone would form an parabolic equation. This intersection line is dependent on the relative position at the center of the cone P(x,y) and the position of the sweep line L. We will call the curved intersection line generated form the different arcs the beach line from now on.
The algorithm must at all times 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) the order of the set of sites (center point on the cones) 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 only maintain a sweep-line status
and we are interested in simulating the discrete event points where
there is a "significant event" occurs. That is, any event that changes the
topological structure of the Voronoi diagram and the beach line.
Two kinds of events appear that could change the beach line. The first is a site event and the second is a circle event.
Site Events
When the sweep line passes over a new site a new arc will be inserted into the beach line.
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 from all the points, by calculation beta j for 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 to find are the two points in the above picture to the left, where the beach line forms 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, with 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 where the two intersections occur by using the quadratic formula. The result is two different x values, that you could plug into any of the two parabolic equations to get the resulting y value, and thereby the two cross points. The equation is given below where 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 there is a pretty big problem in the equation because of the site events, when are you supposed to make the guess about the line in a Voronoi diagram? If you increment it a little before the 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 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 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 (also called vertex event) occurs.
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.
"Which 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 middle arc that would dissapear.
There is a formula for getting the center and the radius of a circle here and if you 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. 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 are other ways of calculating (which are faster) the center point, but I will leave that investigation up to you.
Elementary Functions Used in the Code
Please note, 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. The code 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
I have also created an own circle class with all the nessecary functions inside:
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, by 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
