Click here to Skip to main content
Click here to Skip to main content

Create a Voronoi diagram 1 of 2

, 2 Jul 2012 CPOL
Rate this:
Please Sign up or sign in to vote.
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
        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

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

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

    'Add the two calculated poitns
    result.Points.Add(p1)
    result.Points.Add(p2)

    '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

    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

    '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
        If d < Radius Then
            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

        can.Children.Add(p)
        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.
        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) '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)
        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  

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Kenneth Haugland
Engineer
Norway 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[^].

Comments and Discussions


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.
 
GeneralMy vote of 5 PinmemberSam Dolc17-Aug-13 13:19 
QuestionVoronoi.vb PinmemberTeruteru31420-Mar-13 8:19 
AnswerRe: Voronoi.vb PinmemberKenneth Haugland22-Mar-13 2:51 
GeneralRe: Voronoi.vb PinmemberTeruteru31423-Mar-13 4:28 
GeneralRe: Voronoi.vb PinmemberKenneth Haugland23-Mar-13 4:34 
GeneralMy vote of 5 PinmemberShahan Ayyub17-Mar-13 9:09 
GeneralRe: My vote of 5 PinmemberKenneth Haugland18-Mar-13 6:16 
GeneralMy vote of 5 Pinmemberdim1313-Jan-13 13:28 
GeneralRe: My vote of 5 PinmemberKenneth Haugland10-Mar-13 1:26 
GeneralWell done, Kenneth PinmvpEspen Harlinn15-Aug-12 9:03 
GeneralRe: Well done, Kenneth PinmemberKenneth Haugland15-Aug-12 9:06 
Questionvery nice PinmemberCIDev18-Jul-12 4:46 
AnswerRe: very nice PinmemberKenneth Haugland19-Jul-12 7:32 
GeneralMy vote of 5 Pinmembermanoj kumar choubey15-Jul-12 19:18 
GeneralMy vote of 5 Pinmemberalexander chupeev10-Jul-12 16:39 
GeneralMy vote of 5 PinmemberGregory.Gadow10-Jul-12 4:50 
GeneralRe: My vote of 5 PinmemberKenneth Haugland19-Jul-12 8:57 
GeneralMy vote of 5 PinmemberS Houghtelin9-Jul-12 2:55 
QuestionVote of 5 PinmemberGanesanSenthilvel1-Jul-12 7:41 
GeneralMy vote of 5 Pinmembersam.hill29-Jun-12 13:09 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Terms of Use | Mobile
Web01 | 2.8.141223.1 | Last Updated 2 Jul 2012
Article Copyright 2012 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2014
Layout: fixed | fluid