Create a Voronoi diagram 1 of 2

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