Click here to Skip to main content
11,640,907 members (69,028 online)
Click here to Skip to main content
Add your own
alternative version

Create a Voronoi diagram 1 of 2

, 26 Jul 2012 CPOL 43.8K 1.5K 59
A detailed description for creating Voronoi diagrams, review of basic functions
WPF_VORONOI_TEST.zip
WPF_VORONOI_TEST
WPF_VORONOI_TEST.suo
WPF_VORONOI_TEST
bin
Debug
WPF_VORONOI_TEST.exe
WPF_VORONOI_TEST.pdb
WPF_VORONOI_TEST.vshost.exe
WPF_VORONOI_TEST.vshost.exe.manifest
Release
My Project
MyExtensions
Settings.settings
obj
x86
Debug
DesignTimeResolveAssemblyReferences.cache
DesignTimeResolveAssemblyReferencesInput.cache
MainWindow.baml
TempPE
My Project.Resources.Designer.vb.dll
WPF_VORONOI_TEST.exe
WPF_VORONOI_TEST.g.resources
WPF_VORONOI_TEST.pdb
WPF_VORONOI_TEST.Resources.resources
WPF_VORONOI_TEST.vbproj.GenerateResource.Cache
WPF_VORONOI_TEST_MarkupCompile.cache
WPF_VORONOI_TEST_MarkupCompile.i.cache
Release
WPF_VORONOI_TEST.vbproj.user
Class MainWindow

    'Stores all the points that you can use to create a Voronoi diagram
    Private points As New PointCollection

    Private Sub Canvas1_MouseDown(sender As System.Object, e As System.Windows.Input.MouseButtonEventArgs)
        'Get the location of the mouse and stor it
        Dim LocationOfNewPoint As New Point
        LocationOfNewPoint = e.GetPosition(Canvas1)

        'Stor the point for future use
        points.Add(LocationOfNewPoint)

        'Creat a reprensentation of the point
        Dim NewPoint As New Ellipse
        NewPoint.Fill = Brushes.Black
        NewPoint.StrokeThickness = 2
        NewPoint.Stroke = Brushes.Black

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

        'Add the ellipse to the canvas
        Canvas1.Children.Add(NewPoint)

        'Set the x and y coordinates to get the ellipse centered around the point
        'If you dont substract Radius/2 (with is 2.5 in this case) your center point would be wrong
        Canvas.SetLeft(NewPoint, LocationOfNewPoint.X - 2.5)
        Canvas.SetTop(NewPoint, LocationOfNewPoint.Y - 2.5)

        'Uncomment these lines to see the circle event
        'If points.Count = 3 Then
        '    Dim g As New Circle
        '    g.CreatCircleFromThreePoints(points(0), points(1), points(2))
        '    g.ReturnEllipse(Canvas1)
        'End If

    End Sub

    Private Sub Button1_Click(sender As System.Object, e As System.Windows.RoutedEventArgs) Handles Button1.Click
        'We need to sort the points according to the Y direction
        Dim SortedPointCollection_Y As New PointCollection
        SortedPointCollection_Y = SortPoints(points, False)


        Dim index As Double = -1
        Dim min_distance As Double = 50000

        For i As Integer = 1 To SortedPointCollection_Y.Count - 1

            'Draw the intersection cut
            Canvas1.Children.Add(Intersection(SortedPointCollection_Y(i), SortedPointCollection_Y(i - 1), SortedPointCollection_Y(i).Y))

            'Draw the line from the parabolic cut, this is a new line in the Voronoi diagram
            'although you would have to cut it in most cases
            Canvas1.Children.Add(ParabolicCut(SortedPointCollection_Y(i), SortedPointCollection_Y(i - 1), SortedPointCollection_Y(i).Y))

            'Looping through all the previous points above the sweep line for the circle check
            If i > 1 And i < SortedPointCollection_Y.Count - 1 Then

                Dim TemporaryCircle As New Circle
                'Create a circle of three points
                TemporaryCircle.CreatCircleFromThreePoints(SortedPointCollection_Y(i), SortedPointCollection_Y(i - 1), SortedPointCollection_Y(i - 2))

                'See if the circle does not contain any other points
                '(this check is incomplete as it currently stands
                If Not TemporaryCircle.DoesCircleContainPoint(SortedPointCollection_Y(i + 1)) Then

                    'Draw a green point if an intersection vertex is found
                    TemporaryCircle.ReturnCenterPoint(Canvas1)
                End If
            End If

        Next
    End Sub

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

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

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


#Region "Sorting algorithm"
    ''' <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

#End Region


End Class

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

By viewing downloads associated with this article you agree to the Terms of Service and the article's licence.

If a file you wish to view isn't highlighted, and is a text file (not binary), please let us know and we'll add colourisation support for it.

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[^].

You may also be interested in...

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.150731.1 | Last Updated 26 Jul 2012
Article Copyright 2012 by Kenneth Haugland
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid