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
Submit an article or tip