12,453,354 members (60,557 online)
Article
alternative version

52.7K views
59 bookmarked
Posted

# Create a Voronoi diagram 1 of 2

, 2 Jul 2012 CPOL
 Rate this:
A detailed description for creating Voronoi diagrams, review of basic functions
This is an old version of the currently published article.

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

```''' <summary>
''' Returns a sorted pointcollection based on Lexografically sorting
''' </summary>
''' <param name="samplepoints"></param>
''' <param name="SortByXdirection"></param>
''' <returns></returns>
''' <remarks></remarks>
Private Function SortPoints(ByVal samplepoints As PointCollection, _
ByVal SortByXdirection As Boolean) As PointCollection
'Create another array so we can keep the original array out of order
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
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

'Comparing function
'Returns one of three values - 0 (equal), 1 (greater than), -1 (less than)
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
'Compare X values
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 'Identical points
Return 0
End If
End If
Else
If point1.Y > point2.Y Then
'Compare Y Values
Return 1
ElseIf point1.Y < point2.Y Then
Return -1
Else 'Y1 = Y2
If point1.X > point2.X Then
'Compare Y Values
Return 1
ElseIf point1.X < point2.X Then
Return -1
Else 'Identical points
Return 0
End If
End If
End If
End Function
End Class```

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

```''' <summary>
''' Returns a polyline between the newpoint and the contructed arc of the oldpoint
''' </summary>
''' <param name="NewPoint">New occuring point on the sweep line</param>
''' <param name="OldPoint">Point (above the sweep line) that generates the prarabolic arc</param>
''' <param name="SweepLine">Position of the sweep line (Y) </param>
''' <returns>A straight line between the NewPoint to the Intersection point on the parabolic arc</returns>
''' <remarks></remarks>
Private Function Intersection(ByVal NewPoint As Point, ByVal OldPoint As Point, _
ByVal SweepLine As Double) As Polyline
Dim result As New Polyline

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

```''' <summary>
''' Returns the y value construcked by the event from the newpoint and the contructed arc of the oldpoint
''' </summary>
''' <param name="NewPoint">This is the point on the sweep line</param>
''' <param name="OldPoint">This is one point above the sweep line</param>
''' <param name="SweepLine">Y position of the sweep line</param>
''' <returns>Calculates the distance fromn a point (NewPoint) to the Parabolic
''' arc generated by the OldPoint and the Sweep line</returns>
''' <remarks>The Function only gives you the Y value of the position
''' on the parabolic intersection given the X location</remarks>
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):

```''' <summary>
''' Calculates the line between two Parabolic intersections
''' </summary>
''' <param name="Point1">A point in the Voronoi diagram</param>
''' <param name="point2">A Point in the Voronoi diagram (different from point 1)</param>
''' <param name="ys">The position of the sweep line</param>
''' <returns>A staight line between the two intersection poitns</returns>
''' <remarks>It would give wrong values if the two points have the same or
''' close to the same Y coordinate, as the double has limited significant storage</remarks>
Private Function ParabolicCut(ByVal Point1 As Point, ByVal point2 As Point, ys As Double) As Polyline

'Stores the values in double format, as I didnt bother to rewrite Benjamin Dittes quadratic equation.
Dim x1, y1, x2, y2 As Double

'Inizialize Point 1
x1 = Point1.X
y1 = Point1.Y

'Inizialize Point 2
x2 = point2.X
y2 = point2.Y

'Sweep line, added 500 to make sure the two calculated points are well of the boundaries
ys = ys + 500

'Setting ut calculation of the quadratic equation
Dim a1 As Double = 1 / (2 * (y1 - ys))
Dim a2 As Double = 1 / (2 * (y2 - ys))

'The two resulting values of x from the quadratic equation
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))

'Generate two points to store the two intersection points
Dim p1, p2 As New Point
p1.X = xs1
p2.X = xs2
p1.Y = 0
p2.Y = 0

'Find the resulting Y values calculated from the quadratic equation
' (It dosent matter that the Y is 0 as the function Intersection2
' only uses the X value for computation)
p1.Y = Intersection2(p1, Point1, ys)
p2.Y = Intersection2(p2, Point1, ys)

'Make a polyline to store the result
Dim result As New Polyline

'Making the line visible
result.Stroke = Brushes.Black
result.StrokeThickness = 2

'Return the polyline
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:

```''' <summary>
''' A class for preforming circle tests
''' </summary>
''' <remarks></remarks>
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

Get
End Get
Set(ByVal value As Double)
End Set
End Property

'Calculates the Eucladian distence between two poitns
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

''' <summary>
''' Checks if a Point1 is inside the circle
''' </summary>
''' <param name="p1">Test point</param>
''' <returns>If point is on the circle it is considered to be outside,
''' as three or more points is checked for</returns>
''' <remarks></remarks>
Public Function DoesCircleContainPoint(ByVal p1 As Point) As Boolean
Dim d As Double
d = Distance(p1, CenterPoint)
'Replace this with <= if you just want three points to form a circle
Return True
Else
Return False
End If
End Function

''' <summary>
''' Returns the centerpoint of the created circle
''' </summary>
''' <param name="can"></param>
''' <remarks></remarks>
Public Sub ReturnCenterPoint(ByVal can As Canvas)
Dim p As New Ellipse
p.Fill = Brushes.Blue
p.StrokeThickness = 2
p.Stroke = Brushes.Blue

' Set the width and height of the Ellipse.
p.Width = 5
p.Height = 5

Canvas.SetLeft(p, CenterPoint.X - 2.5)
Canvas.SetTop(p, CenterPoint.Y - 2.5)
End Sub

''' <summary>
''' Returns the complete circle created from the three points
''' </summary>
''' <param name="can"></param>
''' <remarks></remarks>
Public Sub ReturnEllipse(ByVal can As Canvas)
Dim p As New Ellipse
p.Fill = Brushes.Yellow
p.StrokeThickness = 2
p.Stroke = Brushes.Black

' Set the width and height of the Ellipse.

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) 'As RectangleF
'(X - Cx)^2 + (Y - Cy)^2 = R^2
'(R^2 - Cx^2 - Cy^2) + 2X*Cx + 2Y*Cy = X^2 +  Y^2
'Cr + A*X + B*Y = X^2 + Y^2
'Solve matrix using gaussian elimination.
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)
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.

## Share

 Engineer Norway
I hope that you like the stuff I have created and if you do wish to say thank you then a donation is always appreciated.
You can donate here[^].

## You may also be interested in...

 Pro Pro

Discussions posted for the Published version of this article. Posting a message here will take you to the publicly available article in order to continue your conversation in public.

 First Prev Next
 My vote of 5 ISpliter5-Mar-15 14:24 ISpliter 5-Mar-15 14:24
 Re: My vote of 5 Kenneth Haugland12-Mar-15 21:22 Kenneth Haugland 12-Mar-15 21:22
 Re: My vote of 5 ISpliter13-Mar-15 12:38 ISpliter 13-Mar-15 12:38
 My vote of 5 Sam Dolc17-Aug-13 12:19 Sam Dolc 17-Aug-13 12:19
 Voronoi.vb Teruteru31420-Mar-13 7:19 Teruteru314 20-Mar-13 7:19
 Re: Voronoi.vb Kenneth Haugland22-Mar-13 1:51 Kenneth Haugland 22-Mar-13 1:51
 Re: Voronoi.vb Teruteru31423-Mar-13 3:28 Teruteru314 23-Mar-13 3:28
 Re: Voronoi.vb Kenneth Haugland23-Mar-13 3:34 Kenneth Haugland 23-Mar-13 3:34
 My vote of 5 Shahan Ayyub17-Mar-13 8:09 Shahan Ayyub 17-Mar-13 8:09
 Re: My vote of 5 Kenneth Haugland18-Mar-13 5:16 Kenneth Haugland 18-Mar-13 5:16
 My vote of 5 dim1313-Jan-13 12:28 dim131 3-Jan-13 12:28
 Re: My vote of 5 Kenneth Haugland10-Mar-13 0:26 Kenneth Haugland 10-Mar-13 0:26
 Well done, Kenneth Espen Harlinn15-Aug-12 8:03 Espen Harlinn 15-Aug-12 8:03
 Re: Well done, Kenneth Kenneth Haugland15-Aug-12 8:06 Kenneth Haugland 15-Aug-12 8:06
 very nice CIDev18-Jul-12 3:46 CIDev 18-Jul-12 3:46
 Re: very nice Kenneth Haugland19-Jul-12 6:32 Kenneth Haugland 19-Jul-12 6:32
 My vote of 5 manoj kumar choubey15-Jul-12 18:18 manoj kumar choubey 15-Jul-12 18:18
 My vote of 5 Member 1152532514-Mar-15 10:49 Member 11525325 14-Mar-15 10:49
 My vote of 5 alexander chupeev10-Jul-12 15:39 alexander chupeev 10-Jul-12 15:39
 Re: My vote of 5 Kenneth Haugland19-Jul-12 7:57 Kenneth Haugland 19-Jul-12 7:57
 My vote of 5 S Houghtelin9-Jul-12 1:55 S Houghtelin 9-Jul-12 1:55
 Vote of 5 GanesanSenthilvel1-Jul-12 6:41 GanesanSenthilvel 1-Jul-12 6:41
 My vote of 5 sam.hill29-Jun-12 12:09 sam.hill 29-Jun-12 12:09
 Last Visit: 31-Dec-99 18:00     Last Update: 29-Aug-16 1:26 Refresh 1