I’ve tried to get a Delaunay code for VBA. I have to state it, I don’t know why the different approximations I had found of this algorithm perform so badly on VBA. In the web there is a lot of resources to compute a mesh, beeing this post a compilation of a lot of them.
I was not to code my implementation from zero, so started cracking a code that was for CAD software; but it behave erratic and with the mentioned bad performance that I had to modify it hard (creating triangles and points as objects type). Either way, the implementation was very complex for the task, running profusively several costy functions (like sqr) and relying on a not well conceived stack system for the Delaunay condition. My understanding of the pseudocode help was very basic -cause I didn’t want to start from that point-. As it would be better to rework the entire implementation, I left it there to rest.
Some time after I hired some freelancer to convert a fantastic JS Voronoi-Fortuny implementation -that on the browser perform ashtonishing well- to VBA. Beeing a 700 line code I thought it would be easy for an experienced coder. It was not. The freelancer did the job (well nearly done, as it breaks with negative coordinates as he/she tried to draw inside Excel and used structures I had specified not to be used). In the end the code was a perfect JS translation, so perfect that I did learn some ways to bypass differences VBA-JS I thought were impossible to atain. But then came the handicaps… to do the bypass it had to rely in Collections and classes, so performance plummeted even below the CAD implementation. And to break the Collection-Class issue and get the code working was again a hard task. Even more, it returns a Voronoi mesh, when I was looking for a Delaunay one.
Finally, I’m trying my own ideas. From zero. Not the best option, but in the middle I’m working on some points I’ll need in other codes, so did try.
Here is an unfinished code, but you can see it draws the triangles and other related geometries inside an Excel Chart, so I can follow what is doing in every step. There is still some work to do on the stack block for the Delaunay diagonal swap, but things are on the right way:
Option Explicit
Public Type tXYZ
X As Double
Y As Double
Z As Double
End Type
Private Type tTriangle
V1 As Long 'vertex1 pointer
V2 As Long 'vertex2 pointer
V3 As Long 'vertex3 pointer
T1 As Long 'neighbour1 pointer
T2 As Long 'neighbour2 pointer
T3 As Long 'neighbour3 pointer
Xmin As Double
Ymin As Double
Xmax As Double
Ymax As Double
Center As tXYZ
R² As Double
End Type
Dim oT() As tTriangle
Dim oP() As tXYZ
Dim oTTmp As tTriangle
Dim aStack() As Long
Dim aBlnk() As Long
Dim myFrm As UserForm
Private Const PI_8th As Double = 0.392699081698724
Private Const PI As Double = 3.14159265358979
Private Function fDrawTriangles() As Boolean
Dim lgT As Long
For lgT = LBound(oT) To UBound(oT)
'Draw triangle in userform
Call fDrawTriangle(lgT)
Next lgT
End Function
Private Function fDrawTriangle(ByVal lgT As Long) As Boolean
'Draw triangle in userform
'Call fDrawLine(oP(oT(lgT).V1), oP(oT(lgT).V2))
'Call fDrawLine(oP(oT(lgT).V2), oP(oT(lgT).V3))
'Call fDrawLine(oP(oT(lgT).V3), oP(oT(lgT).V1))
End Function
Private Function fDrawLine(ByVal lgP1 As Long, ByVal lgP2 As Long) As Boolean
End Function
Private Function fPointsSort(ByRef oP() As tXYZ, _
Optional ByRef bX As Boolean = True, _
Optional ByRef bAscendent As Boolean = True) As Boolean
Dim oPtTmp() As tXYZ
Dim arrPtr() As Double 'tPointer
Dim lgP As Long
ReDim arrPtr(LBound(oP) To UBound(oP), 1 To 2)
If bX Then
' Sort by X
For lgP = LBound(oP) To UBound(oP)
arrPtr(lgP, 1) = VBA.CLng(lgP) 'Index
arrPtr(lgP, 2) = oP(lgP).X 'Value
Next lgP
Else
' Sort by Y
For lgP = LBound(oP) To UBound(oP)
arrPtr(lgP, 1) = VBA.CLng(lgP) 'Index
arrPtr(lgP, 2) = oP(lgP).Y 'Value
Next lgP
End If
Call fQuickSortArrayDbl(arrPtr(), -1, -1, 2, bAscendent)
oPtTmp() = oP() ' copy source
For lgP = LBound(oP) To UBound(oP)
oP(lgP) = oPtTmp(arrPtr(lgP, 1))
Next lgP
Erase arrPtr()
Erase oPtTmp()
End Function
Private Function fPoints(ByRef oP() As tXYZ) As Long
If Not (Not oP) Then
fPoints = UBound(oP) - LBound(oP) + 1
Else
fPoints = 0
End If
End Function
Private Function fPointsFlip(ByRef oP() As tXYZ) As Boolean
Dim oPTmp() As tXYZ
Dim lgP As Long
If Not (Not oP) Then
oPTmp() = oP()
ReDim Preserve oP(LBound(oP) To UBound(oP) - 1)
For lgP = LBound(oP) To UBound(oP)
oPTmp(UBound(oP) - lgP) = oP(lgP)
Next lgP
oP() = oPTmp()
Erase oPTmp()
End If
End Function
Private Function fPointRemove(ByRef oP() As tXYZ, ByVal lgP As Long) As Boolean
Dim oPTmp() As tXYZ
Dim lgRemove As Long
If Not (Not oP) Then
oPTmp() = oP()
ReDim Preserve oPTmp(LBound(oP) To UBound(oP) - 1)
For lgRemove = (UBound(oP) - 1) To lgP Step -1
oPTmp(lgRemove) = oP(lgRemove + 1)
Next lgRemove
oP() = oPTmp()
Erase oPTmp()
End If
End Function
Private Function fPointAdd(ByRef oP() As tXYZ, ByRef oPt As tXYZ) As Boolean
If Not (Not oP) Then
ReDim Preserve oP(LBound(oP) To UBound(oP) + 1)
Else
ReDim Preserve oP(0)
End If
oP(UBound(oP)) = oPt
End Function
Private Function fConvexHull() As tXYZ()
' Graham-Scan
' https://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/
' The algorithm can be adapted to a 3D space, considering these:
' - the four points A, B, C, D become eight
' - the quadrilateral becomes a polyhedron with eight vertices
' - the rectangle becomes a parallelepiped
Dim oP´() As tXYZ
Dim oP´´() As tXYZ
Dim oPrun() As tXYZ
Dim oPrun´() As tXYZ
Dim oHull() As tXYZ
Dim lgDecalageRGT As Long
Dim lgP As Long
Dim lgPts As Long
Dim lgTop As Long
Dim lgPrun´ As Long
Dim lgPrun´´ As Long
Dim bPrun As Boolean
Dim bPrint As Boolean: bPrint = False
bPrint = False
If fPoints(oP) = 0 Then
' Get list of points...
Call sPointsRnd(100000000, 0, 1000, 0, 1000)
Call fPointsPrint(oP(), 1, 2, 1, -1, -1, bPrint)
End If
Dim dtTimer As Date
dtTimer = VBA.Now()
oPrun() = fPrunePolygon(oP())
Call fPointsPrint(oPrun(), 3, 4, 1, -1, -1, bPrint)
' Get rectangle R:
' X1 = Max(Dx, Ax)
' X2 = Min(Bx, Cx)
' Y1 = Max(By, Ay)
' Y2 = Min(Cy, Dy)
' Rectangle R with vertices: (x2, y1), (x2, y2), (x1, y2), (x1, y1)
oPrun´() = oPrun()
If oPrun(0).X < oPrun(3).X Then oPrun´(0).X = oPrun(3).X Else oPrun´(3).X = oPrun(0).X 'maxX A-D If oPrun(1).X > oPrun(2).X Then oPrun´(1).X = oPrun(2).X Else oPrun´(2).X = oPrun(1).X 'minX B-C
If oPrun(0).Y < oPrun(1).Y Then oPrun´(0).Y = oPrun(1).Y Else oPrun´(1).Y = oPrun(0).Y 'maxY A-B If oPrun(2).Y > oPrun(3).Y Then oPrun´(2).Y = oPrun(3).Y Else oPrun´(3).Y = oPrun(2).Y 'minY C-D
Call fPointsPrint(oPrun´(), 5, 6, 1, -1, -1, bPrint)
'#1st pass
Erase oP´()
oP´() = oP()
oP´´() = oP()
lgPrun´ = LBound(oP´´) - 1
lgPrun´´ = LBound(oP´´) - 1
For lgP = LBound(oP) To UBound(oP)
bPrun = False
If oPrun´(0).X > oP(lgP).X Then
lgPrun´´ = lgPrun´´ + 1
oP´´(lgPrun´´) = oP(lgP)
ElseIf oPrun´(0).Y > oP(lgP).Y Then
lgPrun´´ = lgPrun´´ + 1
oP´´(lgPrun´´) = oP(lgP)
ElseIf oPrun´(2).X < oP(lgP).X Then
lgPrun´´ = lgPrun´´ + 1
oP´´(lgPrun´´) = oP(lgP)
ElseIf oPrun´(2).Y < oP(lgP).Y Then lgPrun´´ = lgPrun´´ + 1 oP´´(lgPrun´´) = oP(lgP) Else lgPrun´ = lgPrun´ + 1 oP´(lgPrun´) = oP(lgP) End If Next lgP ReDim Preserve oP´(0 To lgPrun´) ' outside R ReDim Preserve oP´´(0 To lgPrun´´) ' inside R Call fPointsPrint(oP´(), 7, 8, 1, -1, -1, bPrint) Call fPointsPrint(oP´´(), 9, 10, 1, -1, -1, bPrint) '#2nd pass lgPrun´´ = LBound(oP´´) - 1 For lgP = LBound(oP´´) To UBound(oP´´) If Not fPointInsidePolygon(oP´´(lgP), oPrun()) Then lgPrun´´ = lgPrun´´ + 1 oP´´(lgPrun´´) = oP´´(lgP) End If Next lgP ReDim Preserve oP´´(0 To lgPrun´´) ' outside Q Call fPointsPrint(oP´´(), 11, 12, 1, -1, -1, bPrint) Call fPointsSort(oP´´(), False, True) ReDim Preserve oHull(LBound(oP´´) To UBound(oP´´)) ' Pick the bottom-most (choose the left-most point in case of tie) oHull(LBound(oP)) = oP´´(LBound(oP´´)) lgPts = LBound(oP´´) ' Right side hull For lgP = (LBound(oP´´) + 1) To UBound(oP´´) '!Do While (UBound(oHull) >= 1)
Do While (lgPts >= 2)
' If two or more points make same angle with p0, remove all but the one that is farthest from p0
' (in above sorting our criteria was to keep the farthest point at the end)
'!If CCW(oHull(UBound(oHull) - 1), oHull(UBound(oHull)), oP´´(lgP)) Then Exit Do
If CCW(oHull(lgPts - 2), oHull(lgPts - 1), oP´´(lgP)) Then Exit Do
'Remove UBound(oHull) point from oHull()
'!ReDim Preserve oHull(LBound(oHull) To UBound(oHull) - 1)
lgPts = lgPts - 1
Loop
'Add Point oP´´(lgP) to oHull()
'!ReDim Preserve oHull(LBound(oHull) To UBound(oHull) + 1)
'!oHull(UBound(oHull)) = oP´´(lgP)
oHull(lgPts) = oP´´(lgP)
lgPts = lgPts + 1
Next lgP
' Left side hull
lgTop = lgPts
For lgP = (UBound(oP´´) - 1) To LBound(oP´´) Step -1
'!Do While (UBound(oHullLFT) >= 1)
Do While (lgPts > lgTop)
'!If CCW(oHullLFT(UBound(oHullLFT) - 1), oHullLFT(UBound(oHullLFT)), oP´´(lgP)) Then Exit Do
If CCW(oHull(lgPts - 2), oHull(lgPts - 1), oP´´(lgP)) Then Exit Do
'Remove UBound(oHull) point from oHull()
'!ReDim Preserve oHullLFT(LBound(oHullLFT) To UBound(oHullLFT) - 1)
lgPts = lgPts - 1
Loop
'Add Point oP´´(lgP) to oHull()
'!ReDim Preserve oHullLFT(LBound(oHullLFT) To UBound(oHullLFT) + 1)
'!oHullLFT(UBound(oHullLFT)) = oP´´(lgP)
oHull(lgPts) = oP´´(lgP)
lgPts = lgPts + 1
Next lgP
ReDim Preserve oHull(LBound(oHull) To lgPts - 1)
Call fPointsPrint(oHull(), 13, 14, 1, -1, -1, bPrint)
fConvexHull = oHull()
'Erase oHull()
Call fPointsPrint(oHull(), 15, 16, 1, -1, -1, bPrint)
Debug.Print dtTimer & " vs " & VBA.Now()
End Function
Private Function fPointsPrint(ByRef oP() As tXYZ, _
Optional ByVal lC_X As Long = 1, _
Optional ByVal lC_Y As Long = 2, _
Optional ByVal lR_Start As Long = 1, _
Optional ByVal lgP_Start As Long = -1, _
Optional ByVal lgP_End As Long = -1, _
Optional ByVal bPrint = False) As Boolean
Dim lgP As Long
Dim lgR As Long
If Not bPrint Then Exit Function
If lgP_Start < LBound(oP) Then lgP_Start = LBound(oP)
If lgP_End < LBound(oP) Then lgP_End = UBound(oP)
lgR = lR_Start
If lgP_Start <= lgP_End Then
For lgP = lgP_Start To lgP_End
Cells(lgR, lC_X).Value2 = oP(lgP).X
Cells(lgR, lC_Y).Value2 = oP(lgP).Y
lgR = lgR + 1
Next lgP
Else 'descending
For lgP = lgP_Start To lgP_End Step -1
Cells(lgR, lC_X).Value2 = oP(lgP).X
Cells(lgR, lC_Y).Value2 = oP(lgP).Y
lgR = lgR + 1
Next lgP
End If
End Function
Private Function fPointInsidePolygon(ByRef oPoint As tXYZ, ByRef oPolygon() As tXYZ) As Boolean
Dim lgP As Long
Dim iInside As Integer
' Avoid no segment (coincident points)
For lgP = LBound(oPolygon) To UBound(oPolygon) - 1
If oPolygon(lgP).X = oPolygon(lgP + 1).X And oPolygon(lgP).Y = oPolygon(lgP + 1).Y Then
Call fPointRemove(oPolygon(), lgP)
End If
Next lgP
iInside = fLineSide(oPoint, oPolygon(LBound(oPolygon)), oPolygon(LBound(oPolygon) + 1))
For lgP = LBound(oPolygon) + 1 To UBound(oPolygon) - 1
If iInside <> fLineSide(oPoint, oPolygon(lgP), oPolygon(lgP + 1)) Then Exit Function
Next lgP
If oPolygon(LBound(oPolygon)).X <> oPolygon(UBound(oPolygon)).X Or oPolygon(LBound(oPolygon)).Y <> oPolygon(UBound(oPolygon)).Y Then
If iInside <> fLineSide(oPoint, oPolygon(UBound(oPolygon)), oPolygon(LBound(oPolygon))) Then Exit Function
End If
fPointInsidePolygon = True
End Function
Private Function NewPoint(Optional ByVal X As Double = 0, _
Optional ByVal Y As Double = 0, _
Optional ByVal Z As Double = 0) As tXYZ
With NewPoint
.X = X
.Y = Y
.Z = Z
End With
End Function
Private Function fPrunePolygon(ByRef oP() As tXYZ) As tXYZ()
' A, B, C, D As Points
' Q As Quadrilateral
' R As Rectangle
' DD As Diamond [(Xmid, Ymin)
'
' for each Point p in oP
' Update A, B, C, D --> Q
' # First pass
' Create R from Q
' for each Point p in S
' if p inside R
' prune p from S
' # Second pass
' for each Point p in S
' if p inside Q
' prune p from S
' ^ D / +\
' | /.... + \
' | / / Q x.... \ C
' | / /---------....| \
' | / /| * * | \
' | \ /x|* R * * | + /
' | \ /..|.....--------| /
' | A \ ..x...|/ B
' | \ + + /
' | \ /
' --------------------------->
'
' * points inside R rectangle are prunable
' x point inside Q are prunable but not until the end
' + point outside Q are not prunable
'
If fPoints(oP) = 0 Then
' Get list of points...
Call sPointsRnd(300, 0, 100, 0, 100)
End If
Dim PtA As tXYZ, PtB As tXYZ, PtC As tXYZ, PtD As tXYZ
Dim Q(0 To 3) As tXYZ 'Quadrilateral
Dim R(0 To 3) As tXYZ 'Rectangle
Dim DD(0 To 3) As tXYZ 'Diamond
Dim lgP As Long
Q(0) = oP(LBound(oP)): Q(1) = Q(0): Q(2) = Q(0): Q(3) = Q(0)
For lgP = LBound(oP) To UBound(oP)
With Q(0) ' Point A
If (-.X - .Y) < (-oP(lgP).X - oP(lgP).Y) Then Q(0) = oP(lgP)
End With
With Q(1) ' Point B
If (.X - .Y) < (oP(lgP).X - oP(lgP).Y) Then Q(1) = oP(lgP)
End With
With Q(2) ' Point C
If (.X + .Y) < (oP(lgP).X + oP(lgP).Y) Then Q(2) = oP(lgP)
End With
With Q(3) ' Point D
If (-.X + .Y) < (-oP(lgP).X + oP(lgP).Y) Then Q(3) = oP(lgP)
End With
Next lgP
fPrunePolygon = Q()
End Function
Private Sub sPointsRnd(Optional ByVal lgPts As Long = 0, _
Optional ByVal Xmin As Double = 0, _
Optional ByVal Xmax As Double = 0, _
Optional ByVal Ymin As Double = 0, _
Optional ByVal Ymax As Double = 0)
Dim lgP As Long
Dim dbTmp As Long
If Xmin = 0 Then Xmin = (Rnd() * 100) + 0
If Xmax = 0 Then Xmax = (Rnd() * 100) + 0
If Ymin = 0 Then Ymin = (Rnd() * 100) + 0
If Ymax = 0 Then Ymax = (Rnd() * 100) + 0
If Xmax < Xmin Then
Xmax = dbTmp
Xmax = Xmin
Xmin = dbTmp
End If
If Ymax < Ymin Then
Ymax = dbTmp
Ymax = Ymin
Ymin = dbTmp
End If
ReDim Preserve oP(0 To lgPts - 1)
For lgP = LBound(oP) To UBound(oP)
oP(lgP) = NewPoint((Rnd() * Xmax - Xmin) + Xmin, (Rnd() * Ymax - Ymin) + Ymin)
Next lgP
End Sub
Private Sub sDelaunay()
Dim oChrt As Excel.ChartObject
Dim rgData As Excel.Range
Dim aData As Variant
Dim lgP As Long
Dim lgT As Long
Dim lgTs As Long ' total number of triangles
Dim Xmin As Double, Ymin As Double
Dim Xmax As Double, Ymax As Double
Dim Xmid As Double, Ymid As Double
Dim HSide As Double, VSide As Double
Dim i12 As Integer, i23 As Integer, i31 As Integer
Dim bDalaunay As Boolean
' Elements are sorted by X
With ActiveSheet
Call sPointsCreate(10)
Set rgData = .Range("A1", .Range("B1", .Range("B1").End(xlDown)))
End With
aData = rgData.Value2
'ReDim oP(0 To 2 + (UBound(aData, 1) - LBound(aData, 1) + 1))
ReDim aStack(LBound(oP) To UBound(oP)) As Long
ReDim aBlnk(LBound(oP) To UBound(oP)) As Long
For lgP = LBound(oP) + 3 To UBound(oP)
With oP(lgP)
'.X = aData(lgP - 2, 1)
'.Y = aData(lgP - 2, 2)
If .X < Xmin Then Xmin = .X If .X > Xmax Then Xmax = .X
If .Y < Ymin Then Ymin = .Y If .Y > Ymax Then Ymax = .Y
End With
Next lgP
lgT = -1
' Generate the super-triangle
Xmid = 0.5 * (Xmax + Xmin)
Ymid = 0.5 * (Ymax + Ymin)
HSide = 2 * (Xmid - Xmin)
VSide = 2 * (Ymid - Ymin)
With oP(LBound(oP) + 0)
.X = Xmid - (3 * HSide)
.Y = Ymid - (3 * VSide)
End With
With oP(LBound(oP) + 1)
.X = Xmid
.Y = Ymid + (3 * VSide)
End With
With oP(LBound(oP) + 2)
.X = Xmid + (3 * HSide)
.Y = Ymid - (3 * VSide)
End With
Set oChrt = ActiveSheet.ChartObjects("ChrtPts")
Call fChrtSeriesDelete(oChrt)
' Add supertriangle points
Call fChrtTriangleAdd(oP(LBound(oP) + 0), _
oP(LBound(oP) + 1), _
oP(LBound(oP) + 2), _
oChrt)
lgT = lgT + 1
lgTs = lgT
ReDim Preserve oT(0 To lgT)
With oT(lgT)
.T1 = -1
.T2 = -2
.T3 = -3
.V1 = LBound(oP) + 0
.V2 = LBound(oP) + 1
.V3 = LBound(oP) + 2
End With
' Get new max-min coordinates...
Call fTriangleParameters(lgT)
' Generate new triangles:
For lgP = LBound(oP) + 3 To UBound(oP)
Call fChrtPtHighlight(oChrt, 1, lgP - 3 + 1, True)
'Call fChrtPtAdd(oP(lgP), oChrt)
For lgT = LBound(oT) To UBound(oT)
' Discard point inside boundary
If oP(lgP).X < oT(lgT).Xmax Then If oP(lgP).X > oT(lgT).Xmin Then
If oP(lgP).Y < oT(lgT).Ymax Then If oP(lgP).Y > oT(lgT).Ymin Then
' Point is inside boundary of triangle, so check more intensively
i12 = fLineSide(oP(lgP), oP(oT(lgT).V2), oP(oT(lgT).V1))
i23 = fLineSide(oP(lgP), oP(oT(lgT).V3), oP(oT(lgT).V2))
i31 = fLineSide(oP(lgP), oP(oT(lgT).V1), oP(oT(lgT).V3))
If (i12 = i23 And i23 = i31) Then
ReDim Preserve oT(0 To lgTs + 2)
' Always store points in Counter-clockwise order (last point is new point)
With oT(lgTs + 1)
.T1 = oT(lgT).T2
.T2 = lgTs + 2
.T3 = lgT
.V1 = oT(lgT).V2
.V2 = oT(lgT).V3
.V3 = lgP
' Add supertriangle points
Call fChrtTriangleAdd(oP(.V1), _
oP(.V2), _
oP(.V3), _
ActiveSheet.ChartObjects("ChrtPts"))
End With
With oT(lgTs + 2)
.T1 = oT(lgT).T3
.T2 = lgT
.T3 = lgTs + 1
.V1 = oT(lgT).V3
.V2 = oT(lgT).V1
.V3 = lgP
' Add supertriangle points
Call fChrtTriangleAdd(oP(.V1), _
oP(.V2), _
oP(.V3), _
ActiveSheet.ChartObjects("ChrtPts"))
End With
' Reformat initial triangle
With oT(lgT)
'.T1 = oT(lgT).T1 ' this neighbour is not rewritten
.T2 = lgTs + 1
.T3 = lgTs + 2
'.V1 = .V1 ' this vertex is not rewritten
'.V2 = .V2 ' this vertex is not rewritten
.V3 = lgP
End With
' Get new max-min coordinates for final triangles...
Call fTriangleParameters(lgT)
Call fTriangleParameters(lgTs + 1)
Call fTriangleParameters(lgTs + 2)
lgTs = lgTs + 2
' Is the optimum delaunay?
' Each new triangle has three neighbours that have to be checked for Delaunay condition
' For Neighbour triangle 1
'call fDelaunay(lgT, oT(lgT).T1)
' For Neighbour triangle 2
'call fDelaunay(lgT, oT(lgT).T2)
' For Neighbour triangle 3
'Call fDelaunay(lgT, oT(lgT).T3)
Exit For ' Goto NextTriangle
End If
End If
End If
End If
End If
Next lgT
Call fChrtPtHighlight(oChrt, 1, lgP - 3, False)
Next lgP
' Delete triangles on the supertriangle structure
Dim oT_() As tTriangle: oT_() = oT()
For lgT = LBound(oT) To UBound(oT)
With oT(lgT)
If .V1 >= 0 And _
.V2 >= 0 And _
.V3 >= 0 Then
'lgT_ = lgT_ + 1
'oT_(0 to lgT) = oT(lgT)
End If
End With
Next lgT
'Redim preserve oT_(0 to lgT_)
'oT() = oT_()
'Erase oT_()
Stop
End Sub
Private Function fDelaunay(ByVal lgT1 As Long, _
ByVal lgT2 As Long) As Boolean
Dim bSwap As Boolean
Dim TPtrTmp As Long
' Rotate vertices of both triangles so in both triangles vertices 3 are opposed
If oT(lgT1).V1 = oT(lgT2).V2 And oT(lgT1).V2 = oT(lgT2).V1 Then
' Do nothing...
ElseIf oT(lgT1).V1 = oT(lgT2).V1 And oT(lgT1).V2 = oT(lgT2).V3 Then
' rotate 2 clockwise
Call fTriangleRotate(lgT2, False)
ElseIf oT(lgT1).V1 = oT(lgT2).V3 And oT(lgT1).V2 = oT(lgT2).V2 Then
' rotate 2 counter-clockwise
Call fTriangleRotate(lgT2, True)
'--
ElseIf oT(lgT1).V2 = oT(lgT2).V1 And oT(lgT1).V3 = oT(lgT2).V3 Then
' rotate 1 counter-clockwise
Call fTriangleRotate(lgT1, True)
' rotate 2 clockwise
Call fTriangleRotate(lgT2, False)
ElseIf oT(lgT1).V2 = oT(lgT2).V2 And oT(lgT1).V3 = oT(lgT2).V1 Then
' rotate 1 counter-clockwise
Call fTriangleRotate(lgT1, True)
ElseIf oT(lgT1).V2 = oT(lgT2).V3 And oT(lgT1).V3 = oT(lgT2).V2 Then
' rotate 1 counter-clockwise
Call fTriangleRotate(lgT1, True)
' rotate 2 counter-clockwise
Call fTriangleRotate(lgT2, True)
'--
ElseIf oT(lgT1).V3 = oT(lgT2).V1 And oT(lgT1).V1 = oT(lgT2).V3 Then
' rotate 1 clockwise
Call fTriangleRotate(lgT1, False)
' rotate 2 clockwise
Call fTriangleRotate(lgT2, False)
ElseIf oT(lgT1).V3 = oT(lgT2).V2 And oT(lgT1).V1 = oT(lgT2).V1 Then
' rotate 1 clockwise
Call fTriangleRotate(lgT1, False)
ElseIf oT(lgT1).V3 = oT(lgT2).V3 And oT(lgT1).V1 = oT(lgT2).V2 Then
' rotate 1 clockwise
Call fTriangleRotate(lgT1, False)
' rotate 2 counter-clockwise
Call fTriangleRotate(lgT2, True)
End If
If fDistance2DPoints(oP(oT(lgT2).V2), oT(lgT1).Center) < EPSILON Then bSwap = True TPtrTmp = oT(lgT1).T2 'Destroy oT(lgT1).T1 and oT(lgT2).T1 oT(lgT1).T1 = oT(lgT2).T2 oT(lgT1).T2 = lgT2 'oT(lgT1).T3 = does not change oT(lgT2).T1 = TPtrTmp oT(lgT2).T2 = lgT1 'oT(lgT2).T3 = does not change ' Swap vertices oT(lgT1).V2 = oT(lgT2).V3 oT(lgT2).V2 = oT(lgT1).V3 End If If bSwap Then Call fTriangleParameters(lgT1) Call fTriangleParameters(lgT2) ' 'Call fDelaunay(lgT1, oT(lgT1).T1) Then 'Call fDelaunay(lgT1, oT(lgT1).T2) Then 'Call fDelaunay(lgT1, oT(lgT1).T3) Then ' 'Call fDelaunay(lgT2, oT(lgT2).T1) Then 'Call fDelaunay(lgT2, oT(lgT2).T2) Then 'Call fDelaunay(lgT2, oT(lgT2).T3) Then Else fDelaunay = True End If End Function Private Function fTriangleRotate(ByRef lgT As Long, _ Optional ByVal bCCW As Boolean = True) As Boolean ' Given a triangle, rotate vertices oTTmp = oT(lgT) With oT(lgT) If bCCW Then ' counter-clockwise .V1 = oTTmp.V2 .V2 = oTTmp.V3 .V3 = oTTmp.V1 .T1 = oTTmp.T2 .T2 = oTTmp.T3 .T3 = oTTmp.T1 Else .V1 = oTTmp.V3 .V2 = oTTmp.V1 .V3 = oTTmp.V2 .T1 = oTTmp.T3 .T2 = oTTmp.T1 .T3 = oTTmp.T2 End If End With End Function Private Function fDelaunayDiagonal(ByRef lgT1 As Long, _ ByRef lgT2 As Long) As Boolean ' Given two triangles, find the opposed vertices ' Rotate triangles so the 3 vertices are opposed... Dim dX1 As Double, dY1 As Double Dim dX2 As Double, dY2 As Double dX1 = oP(oT(lgT1).V3).X - oP(oT(lgT2).V3).X dY1 = oP(oT(lgT1).V3).Y - oP(oT(lgT2).V3).Y dX2 = oP(oT(lgT1).V1).X - oP(oT(lgT2).V2).X dY2 = oP(oT(lgT1).V1).Y - oP(oT(lgT2).V2).Y If ((dX1 * dX1) + (dY1 * dY1)) > ((dX2 * dX2) + (dY2 * dY2)) Then
Call fSwapDiagonal(lgT1, lgT2)
End If
End Function
Private Function fSwapDiagonal(ByRef lgT1 As Long, _
ByRef lgT2 As Long) As Boolean
' For given diagonal D1 on triangle T1 = diagonal D2 on triangle T2, will swap diagonals with opposed vertices
' Rotate triangles so the 3 vertices are opposed...
Dim TPtrTmp As Long
With oT(lgT1)
TPtrTmp = .T2
'Destroy oT(lgT1).T1 and oT(lgT2).T1
.T1 = oT(lgT2).T2
.T2 = lgT2
'.T3 = does not change
End With
With oT(lgT2)
.T1 = TPtrTmp
.T2 = lgT1
'.T3 = does not change
End With
' Swap vertices
oT(lgT1).V2 = oT(lgT2).V3
oT(lgT2).V2 = oT(lgT1).V3
End Function
Private Function fNeigbourSide(ByRef iSide As Integer, _
ByRef lgT1 As Long, _
ByRef lgT2 As Long) As Integer
' For given iSide on triangle1, return neighbour side on triangle2.
' Both triangles turn the same: 1->2->3
If iSide = 1 Then
If oT(lgT1).V1 = oT(lgT2).V1 Then
' fNeigbourSide = 3
ElseIf oT(lgT1).V1 = oT(lgT2).V2 Then
' fNeigbourSide = 2
ElseIf oT(lgT1).V1 = oT(lgT2).V3 Then
' fNeigbourSide = 1
End If
ElseIf iSide = 2 Then
If oT(lgT1).V2 = oT(lgT2).V1 Then
' fNeigbourSide = 3
ElseIf oT(lgT1).V2 = oT(lgT2).V2 Then
' fNeigbourSide = 2
ElseIf oT(lgT1).V2 = oT(lgT2).V3 Then
' fNeigbourSide = 1
End If
ElseIf iSide = 3 Then
If oT(lgT1).V3 = oT(lgT2).V1 Then
' fNeigbourSide = 3
ElseIf oT(lgT1).V3 = oT(lgT2).V2 Then
' fNeigbourSide = 2
ElseIf oT(lgT1).V3 = oT(lgT2).V3 Then
' fNeigbourSide = 1
End If
End If
End Function
Private Function fTriangleFromPoints(ByRef oP1 As tXYZ, _
ByRef oP2 As tXYZ, _
ByRef oP3 As tXYZ) As tTriangle
With fTriangleFromPoints
.V1 = 1
.V2 = 2
.V3 = 3
.T1 = -1
.T2 = -2
.T3 = -3
.Xmin = oP1.X
.Xmax = oP1.X
.Ymin = oP1.Y
.Ymax = oP1.Y
If .Xmin > oP2.X Then .Xmin = oP2.X
If .Xmin > oP3.X Then .Xmin = oP3.X
If .Xmax < oP2.X Then .Xmax = oP2.X
If .Xmax < oP3.X Then .Xmax = oP3.X If .Ymin > oP2.Y Then .Ymin = oP2.Y
If .Ymin > oP3.Y Then .Ymin = oP3.Y
If .Ymax < oP2.Y Then .Ymax = oP2.Y
If .Ymax < oP3.Y Then .Ymax = oP3.Y Dim B As tXYZ Dim C As tXYZ Dim D As Double Dim Bx² As Double Dim By² As Double Dim Cx² As Double Dim Cy² As Double B.X = oP2.X - oP1.X B.Y = oP2.Y - oP1.Y C.X = oP3.X - oP1.X C.Y = oP3.Y - oP1.Y Bx² = B.X * B.X: By² = B.Y * B.Y Cx² = C.X * C.X: Cy² = C.Y * C.Y D = 1 / (2 * (B.X * C.Y - B.Y * C.X)) .Center.X = D * (C.Y * (Bx² + By²) - B.Y * (Cx² + Cy²)) .Center.Y = D * (B.X * (Cx² + Cy²) - C.X * (Bx² + By²)) .R² = (.Center.X * .Center.X) + (.Center.Y * .Center.Y) .Center.X = .Center.X + oP1.X .Center.Y = .Center.Y + oP1.Y End With End Function Private Function fTriangleParameters(ByRef lgT As Long) As Boolean With oT(lgT) .Xmin = oP(.V1).X .Xmax = oP(.V1).X .Ymin = oP(.V1).Y .Ymax = oP(.V1).Y If .Xmin > oP(.V2).X Then .Xmin = oP(.V2).X
If .Xmin > oP(.V3).X Then .Xmin = oP(.V3).X
If .Xmax < oP(.V2).X Then .Xmax = oP(.V2).X
If .Xmax < oP(.V3).X Then .Xmax = oP(.V3).X If .Ymin > oP(.V2).Y Then .Ymin = oP(.V2).Y
If .Ymin > oP(.V3).Y Then .Ymin = oP(.V3).Y
If .Ymax < oP(.V2).Y Then .Ymax = oP(.V2).Y
If .Ymax < oP(.V3).Y Then .Ymax = oP(.V3).Y
.Center = Circumcenter(oP(.V1).X, oP(.V1).Y, oP(.V2).X, oP(.V2).Y, oP(.V3).X, oP(.V3).Y)
.R² = (.Center.X - oP(.V1).X) ^ 2 + (.Center.Y - oP(.V1).Y) ^ 2
'Dim B As tXYZ
'Dim C As tXYZ
'Dim D As Double
'Dim Bx² As Double
'Dim By² As Double
'Dim Cx² As Double
'Dim Cy² As Double
'
'B.X = oP(.V2).X - oP(.V1).X
'B.Y = oP(.V2).Y - oP(.V1).Y
'C.X = oP(.V3).X - oP(.V1).X
'C.Y = oP(.V3).Y - oP(.V1).Y
'Bx² = B.X * B.X: By² = B.Y * B.Y
'Cx² = C.X * C.X: Cy² = C.Y * C.Y
'D = 1 / (2 * (B.X * C.Y - B.Y * C.X))
'.Center.X = D * (C.Y * (Bx² + By²) - B.Y * (Cx² + Cy²))
'.Center.Y = D * (B.X * (Cx² + Cy²) - C.X * (Bx² + By²))
'.R² = (.Center.X * .Center.X) + (.Center.Y * .Center.Y)
'.Center.X = .Center.X + oP(.V1).X
'.Center.Y = .Center.Y + oP(.V1).Y
End With
End Function
Private Function Circumcenter(ByVal x1 As Double, ByVal y1 As Double, _
ByVal x2 As Double, ByVal y2 As Double, _
ByVal x3 As Double, ByVal y3 As Double) As tXYZ
Dim A As Double
Dim B As Double
Dim C As Double
Dim D As Double
A = x1 * x1 + y1 * y1
B = x2 * x2 + y2 * y2
C = x3 * x3 + y3 * y3
D = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
If D <> 0 Then
With Circumcenter
.X = (A * (y2 - y3) + B * (y3 - y1) + C * (y1 - y2)) / D
.Y = (A * (x3 - x2) + B * (x1 - x3) + C * (x2 - x1)) / D
End With
End If
End Function
Private Sub sLineSide()
Dim oPoint As tXYZ, oPt1 As tXYZ, oPt2 As tXYZ
With oPoint
.X = 1
.Y = -1
End With
With oPt1
.X = 0
.Y = 0
End With
With oPt2
.X = 10
.Y = 0
End With
Debug.Print fLineSide(oPoint, oPt1, oPt2)
End Sub
Private Function CCW(ByRef oPt1 As tXYZ, _
ByRef oPt2 As tXYZ, _
ByRef oPt3 As tXYZ) As Boolean
' If counter clock-wise, then CCW is true
CCW = ((oPt2.X - oPt1.X) * (oPt3.Y - oPt1.Y)) > ((oPt2.Y - oPt1.Y) * (oPt3.X - oPt1.X))
End Function
Private Function fLineSide(ByRef oPoint As tXYZ, _
ByRef oPt1 As tXYZ, _
ByRef oPt2 As tXYZ) As Integer
' Use the sign of the determinant of vectors (AB,AM), where M(X,Y) is the query point:
' Position = Sign((Bx - Ax) * (Y - Ay) - (By - Ay) * (X - Ax))
' It is 0 on the line, and -1 on right side, +1 on the left side.
fLineSide = VBA.Sgn((oPt2.X - oPt1.X) * (oPoint.Y - oPt1.Y) - (oPt2.Y - oPt1.Y) * (oPoint.X - oPt1.X))
End Function
Private Function fChrtSeriesDelete(ByVal oChrt As Excel.ChartObject)
Dim lgSeries As Long
With oChrt
With .Chart.SeriesCollection
For lgSeries = .Count To 2 Step -1
oChrt.Chart.SeriesCollection(lgSeries).Delete
Next lgSeries
End With
End With
End Function
Private Function fChrtCircleAdd(ByRef oCenter As tXYZ, _
ByRef Radius As Double, _
Optional ByVal oChrt As Excel.ChartObject) As Boolean
Dim oSeries As Excel.Series
Dim lgAngle As Long
Dim dbAngleRAD As Double
Dim strX As String
Dim strY As String
With oChrt
With .Chart
Set oSeries = .SeriesCollection.NewSeries
With oSeries
For lgAngle = 0 To 17
dbAngleRAD = lgAngle * PI_8th
strX = strX & Replace(oCenter.X + Radius * Cos(dbAngleRAD), ",", ".") & ","
strY = strY & Replace(oCenter.Y + Radius * Sin(dbAngleRAD), ",", ".") & ","
Next lgAngle
strX = VBA.Left$(strX, Len(strX) - 1)
strY = VBA.Left$(strY, Len(strY) - 1)
.XValues = "={" & strX & "}"
.Values = "={" & strY & "}"
With .Format.Line
.Visible = msoTrue
'.ForeColor.ObjectThemeColor = msoThemeColorAccent1
'.ForeColor.TintAndShade = 0
'.ForeColor.Brightness = 0
.Weight = 0.25
.Visible = msoTrue
.DashStyle = msoLineSysDash
End With
End With
End With
End With
End Function
Private Function fChrtTriangleAdd(ByRef oPt1 As tXYZ, _
ByRef oPt2 As tXYZ, _
ByRef oPt3 As tXYZ, _
Optional ByVal oChrt As Excel.ChartObject) As Boolean
Dim oSeries As Excel.Series
With oChrt
With .Chart
Set oSeries = .SeriesCollection.NewSeries
With oSeries
.XValues = "={" & Replace(oPt1.X, ",", ".") & "," & Replace(oPt2.X, ",", ".") & "," & Replace(oPt3.X, ",", ".") & "," & Replace(oPt1.X, ",", ".") & "}"
.Values = "={" & Replace(oPt1.Y, ",", ".") & "," & Replace(oPt2.Y, ",", ".") & "," & Replace(oPt3.Y, ",", ".") & "," & Replace(oPt1.Y, ",", ".") & "}"
With .Format.Line
.Visible = msoTrue
.Weight = 0.25
.DashStyle = msoLineSysDash
'.ForeColor.ObjectThemeColor = msoThemeColorAccent1
'.ForeColor.TintAndShade = 0
'.ForeColor.Brightness = 0
End With
End With
End With
End With
End Function
Private Function fChrtPtAdd(ByRef oPt As tXYZ, _
Optional ByVal oChrt As Excel.ChartObject)
Dim oSeries As Excel.Series
With oChrt
With .Chart
Set oSeries = .SeriesCollection.NewSeries
With oSeries
.XValues = "={" & Replace(oPt.X, ",", ".") & "}"
.Values = "={" & Replace(oPt.Y, ",", ".") & "}"
With .Format.Fill
.Visible = msoTrue
.ForeColor.RGB = RGB(255, 0, 0)
'.ForeColor.ObjectThemeColor = msoThemeColorAccent1
'.ForeColor.TintAndShade = 0
'.ForeColor.Brightness = 0
End With
End With
End With
End With
End Function
Private Function fChrtPtHighlight(ByVal oChrt As Excel.ChartObject, _
ByVal lgSeries As Long, _
ByVal lgData As Long, _
Optional ByVal bActive As Boolean = False)
With oChrt
With .Chart.FullSeriesCollection(lgSeries).Points(lgData)
With .Format
With .Fill
.Visible = msoTrue
If bActive Then
.ForeColor.RGB = RGB(255, 0, 0)
Else
.ForeColor.ObjectThemeColor = msoThemeColorAccent1
'.ForeColor.TintAndShade = 0
'.ForeColor.Brightness = 0
End If
.Transparency = 0
.Solid
End With
End With
'.ApplyDataLabels
End With
End With
End Function
Private Sub sPointsCreate(ByVal lgPoints As Long)
Dim Xmin As Double, Ymin As Double
Dim Xmax As Double, Ymax As Double
Dim lgPoint As Long
Dim dX As Double, dY As Double
Dim lgR As Long
Dim rgData As Excel.Range
Dim XValue As Excel.Range
Dim YValue As Excel.Range
Dim oChrt As Excel.ChartObject
' Create points and print out to range
With ActiveSheet
Xmin = 0
Xmax = 1000
Ymin = 0
Ymax = 1000
dX = (Xmax - Xmin)
dY = (Ymax - Ymin)
ReDim oP(0 To 2 + lgPoints)
For lgPoint = 3 To (lgPoints - 1) + 3
oP(lgPoint).X = Xmin + (Rnd() * dX)
oP(lgPoint).Y = Ymin + (Rnd() * dY)
lgR = lgPoint - 1 + 3
.Cells(lgR, 1).Value2 = oP(lgPoint).X
.Cells(lgR, 2).Value2 = oP(lgPoint).Y
Next lgPoint
Set XValue = .Range("A1", .Range("A1", .Range("A1").End(xlDown)))
'XValue.Select
Set YValue = .Range("B1", .Range("B1", .Range("B1").End(xlDown)))
'YValue.Select
Set rgData = Union(XValue, YValue)
' Sort points by X
Call fPointsSort(oP(), True, True)
With .Sort
' .SortFields.Clear
' .SortFields.Add _
Key:=rgData, _
SortOn:=xlSortOnValues, _
Order:=xlAscending, _
DataOption:=xlSortNormal
' .SetRange rgData
' .Header = xlGuess
' .MatchCase = False
' .Orientation = xlTopToBottom
' .SortMethod = xlPinYin
' '.Apply
End With
' Create chart
'Set oChrt = .Shapes.AddChart2(240, xlXYScatter).ChartObject
With .Shapes.AddChart2(240, xlXYScatter)
.Name = "ChrtPts"
.Left = ActiveSheet.Range("C1").Left
.Top = ActiveSheet.Range("C1").Top
.Height = 800
.Width = 800
With .Chart
With .PlotArea
'.Height = 700
'.Width = 700
End With
.SetSourceData Source:=rgData
.ChartTitle.Delete
With .Axes(xlValue)
'.MinimumScaleIsAuto = True
'.MaximumScaleIsAuto = True
.MinimumScale = 0
.MaximumScale = 1000
.MajorUnit = 250
.MinorUnit = 50
End With
With .Axes(xlCategory)
'.MinimumScaleIsAuto = True
'.MaximumScaleIsAuto = True
.MinimumScale = 0
.MaximumScale = 1000
.MajorUnit = 250
.MinorUnit = 50
End With
End With
End With
End With
End Sub
'-----------------------------------------
' T E S T F U N C T I O N S
'-----------------------------------------
Private Sub sTestPoint()
Dim oTriangle1 As tTriangle
Dim oTriangle2 As tTriangle
Dim oP() As tXYZ
Dim lgP As Long
'Set oChrt = ActiveSheet.ChartObjects("ChrtPts")
'Call sPointsCreate(4)
ReDim oP(1 To 4)
oP(1) = NewPoint(50, 1000)
oP(2) = NewPoint(500, 3200)
oP(3) = NewPoint(-2000, 3500)
oP(4) = NewPoint(-2000, -3500)
For lgP = LBound(oP) To UBound(oP)
Cells(lgP, 1).Value2 = oP(lgP).X
Cells(lgP, 2).Value2 = oP(lgP).Y
Next lgP
oTriangle1 = fTriangleFromPoints(oP(1), oP(2), oP(3))
With oTriangle1
.V1 = 1
.V2 = 2
.V3 = 3
End With
Call fChrtTriangleAdd(oP(1), _
oP(2), _
oP(3), _
ActiveSheet.ChartObjects("ChrtPts"))
oTriangle2 = fTriangleFromPoints(oP(3), oP(2), oP(1))
With oTriangle2
.V1 = 3
.V2 = 2
.V3 = 4
End With
Call fChrtTriangleAdd(oP(3), _
oP(2), _
oP(4), _
ActiveSheet.ChartObjects("ChrtPts"))
'Call fChrtCircleAdd(oTriangle.Center, _
Sqr(oTriangle.R²), _
ActiveSheet.ChartObjects("ChrtPts"))
Stop
End Sub
'----------------------------------------------
' T E S T
'----------------------------------------------
Private Function xTriangle(ByRef oP1 As tXYZ, ByRef oP2 As tXYZ, ByRef oP3 As tXYZ) As tXYZ
With xTriangle
Dim B As tXYZ
Dim C As tXYZ
Dim D As Double
Dim Bx² As Double
Dim By² As Double
Dim Cx² As Double
Dim Cy² As Double
B.X = oP2.X - oP1.X
B.Y = oP2.Y - oP1.Y
C.X = oP3.X - oP1.X
C.Y = oP3.Y - oP1.Y
Bx² = B.X * B.X: By² = B.Y * B.Y
Cx² = C.X * C.X: Cy² = C.Y * C.Y
D = (2 * (B.X * C.Y - B.Y * C.X))
If D <> 0 Then
D = 1 / D
.X = D * (C.Y * (Bx² + By²) - B.Y * (Cx² + Cy²))
.Y = D * (B.X * (Cx² + Cy²) - C.X * (Bx² + By²))
'.R² = (.Center.X * .Center.X) + (.Center.Y * .Center.Y)
.X = .X + oP1.X
.Y = .Y + oP1.Y
End If
End With
End Function
Sub sTestingPerformance()
Dim lgT As Long
Dim dtDate As Date
Dim oP1 As tXYZ, oP2 As tXYZ, oP3 As tXYZ
oP1 = NewPoint(50, 1000)
oP2 = NewPoint(500, 3200)
oP3 = NewPoint(-2000, 3500)
dtDate = VBA.Now()
For lgT = 1 To 100000000
Call xTriangle(oP1, oP2, oP3)
Next lgT
Debug.Print dtDate & " --- " & VBA.Now()
End Sub
In this code there is a convex-hull algorithm implementation, based on the Graham-Scan Rosetta-Stone VB.Net version, which was not working as expected -not at all I would say-, and finally achieved with the help of this post.
The Graham-Scam algorithm is a method of computing the convex hull of a finite set of points in the plane -with time complexity O(n log n)-. The algorithm finds all vertices of the convex hull ordered along its boundary.
- The first step in this algorithm is to find the point with the lowest y-coordinate.
- Next, the set of points should be sorted in increasing order of the angle and the point P made with the x-axis. As this is computationally prone to errors and with bad performance, is a better option to sort points by Y value, and divide the hull construction in left and right halves.
- The algorithm proceeds by considering each of the points in the sorted array in sequence. For each point, it is determined whether moving from the two previously considered points to this point is a “left turn” or a “right turn”. If it is a “right turn”, this means that the second-to-last point is not part of the convex hull and should be removed from consideration. This process is continued for as long as the set of the last three points is a “right turn”. As soon as a “left turn” is encountered, the algorithm moves on to the next point in the sorted array.