Advertisement
Ben_S

VBA code for Delaunay and ConvexHull

May 17th, 2023 (edited)
1,892
0
Never
1
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
VBScript 42.89 KB | Source Code | 0 0
  1. Option Explicit
  2.  
  3. ' VBA code for Delaunay and ConvexHull
  4. ' Created by auadmin on http://www.audeser.com/vba-delaunay-convexhull/
  5. '
  6. ' The VBA formatting was scrambled so this is my attempt to sort it out.
  7. ' Ben Sacherich - 5/17/2023
  8. ' https://pastebin.com/BJvEUnD8
  9. '
  10. ' Constant EPSILON was not defined.
  11. ' Function fQuickSortArrayDbl() was missing.
  12. ' Function fDistance2DPoints() was missing.
  13. '
  14.  
  15. Public Type tXYZ
  16.     x As Double
  17.     y As Double
  18.     Z As Double
  19. End Type
  20.  
  21. Private Type tTriangle
  22.     V1 As Long 'vertex1 pointer
  23.    V2 As Long 'vertex2 pointer
  24.    V3 As Long 'vertex3 pointer
  25.    T1 As Long 'neighbour1 pointer
  26.    T2 As Long 'neighbour2 pointer
  27.    T3 As Long 'neighbour3 pointer
  28.    
  29.     Xmin As Double
  30.     Ymin As Double
  31.     Xmax As Double
  32.     Ymax As Double
  33.    
  34.     Center As tXYZ
  35.     R² As Double
  36. End Type
  37.  
  38. Dim oT() As tTriangle
  39. Dim oP() As tXYZ
  40. Dim oTTmp As tTriangle
  41. Dim aStack() As Long
  42. Dim aBlnk() As Long
  43. ' Dim myFrm As UserForm     ' This is not used and won't compile so Ben commented it out.
  44.  
  45. Private Const PI_8th As Double = 0.392699081698724
  46. Private Const PI As Double = 3.14159265358979
  47.  
  48. '# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
  49. Const EPSILON = 1       ' ### Placeholder for missing code ###
  50.        
  51. Function fDistance2DPoints(x As tXYZ, y As tXYZ)
  52.     ' ### Placeholder for missing function ###
  53.    Stop
  54. End Function
  55.        
  56. Function fQuickSortArrayDbl(arry As Variant, i, j, k, bAscendent As Boolean) As Variant
  57.     ' ### Placeholder for missing function ###
  58.    Stop
  59. End Function
  60. '# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
  61.  
  62. Private Function fDrawTriangles() As Boolean
  63. ' # This function is not used #
  64.    Dim lgT As Long
  65.  
  66.     For lgT = LBound(oT) To UBound(oT)
  67.         'Draw triangle in userform
  68.        Call fDrawTriangle(lgT)
  69.     Next lgT
  70. End Function
  71.  
  72. Private Function fDrawTriangle(ByVal lgT As Long) As Boolean
  73. ' # This function is not used #
  74.  
  75.     'Draw triangle in userform
  76.    'Call fDrawLine(oP(oT(lgT).V1), oP(oT(lgT).V2))
  77.    'Call fDrawLine(oP(oT(lgT).V2), oP(oT(lgT).V3))
  78.    'Call fDrawLine(oP(oT(lgT).V3), oP(oT(lgT).V1))
  79. End Function
  80.  
  81. Private Function fDrawLine(ByVal lgP1 As Long, ByVal lgP2 As Long) As Boolean
  82. ' # This function was missing but is not needed. #
  83. End Function
  84.  
  85. Private Function fPointsSort(ByRef oP() As tXYZ, _
  86.                              Optional ByRef bX As Boolean = True, _
  87.                              Optional ByRef bAscendent As Boolean = True) As Boolean
  88.     Dim oPtTmp() As tXYZ
  89.     Dim arrPtr() As Double 'tPointer
  90.    Dim lgP As Long
  91.  
  92.     ReDim arrPtr(LBound(oP) To UBound(oP), 1 To 2)
  93.     If bX Then
  94.     ' Sort by X
  95.        For lgP = LBound(oP) To UBound(oP)
  96.             arrPtr(lgP, 1) = VBA.CLng(lgP) 'Index
  97.            arrPtr(lgP, 2) = oP(lgP).x     'Value
  98.        Next lgP
  99.     Else
  100.     ' Sort by Y
  101.        For lgP = LBound(oP) To UBound(oP)
  102.             arrPtr(lgP, 1) = VBA.CLng(lgP) 'Index
  103.            arrPtr(lgP, 2) = oP(lgP).y     'Value
  104.        Next lgP
  105.     End If
  106.     Call fQuickSortArrayDbl(arrPtr(), -1, -1, 2, bAscendent)
  107.  
  108.     oPtTmp() = oP() ' copy source
  109.    For lgP = LBound(oP) To UBound(oP)
  110.         oP(lgP) = oPtTmp(arrPtr(lgP, 1))
  111.     Next lgP
  112.     Erase arrPtr()
  113.     Erase oPtTmp()
  114.  
  115. End Function
  116.  
  117. Private Function fPoints(ByRef oP() As tXYZ) As Long
  118.     If Not (Not oP) Then
  119.         fPoints = UBound(oP) - LBound(oP) + 1
  120.     Else
  121.         fPoints = 0
  122.     End If
  123. End Function
  124.  
  125. Private Function fPointsFlip(ByRef oP() As tXYZ) As Boolean
  126.     Dim oPTmp() As tXYZ
  127.     Dim lgP As Long
  128.    
  129.     If Not (Not oP) Then
  130.         oPTmp() = oP()
  131.         ReDim Preserve oP(LBound(oP) To UBound(oP) - 1)
  132.         For lgP = LBound(oP) To UBound(oP)
  133.             oPTmp(UBound(oP) - lgP) = oP(lgP)
  134.         Next lgP
  135.         oP() = oPTmp()
  136.         Erase oPTmp()
  137.     End If
  138. End Function
  139.  
  140. Private Function fPointRemove(ByRef oP() As tXYZ, ByVal lgP As Long) As Boolean
  141.     Dim oPTmp() As tXYZ
  142.     Dim lgRemove As Long
  143.    
  144.     If Not (Not oP) Then
  145.         oPTmp() = oP()
  146.         ReDim Preserve oPTmp(LBound(oP) To UBound(oP) - 1)
  147.         For lgRemove = (UBound(oP) - 1) To lgP Step -1
  148.             oPTmp(lgRemove) = oP(lgRemove + 1)
  149.         Next lgRemove
  150.         oP() = oPTmp()
  151.         Erase oPTmp()
  152.     End If
  153. End Function
  154.  
  155. Private Function fPointAdd(ByRef oP() As tXYZ, ByRef oPt As tXYZ) As Boolean
  156.     If Not (Not oP) Then
  157.         ReDim Preserve oP(LBound(oP) To UBound(oP) + 1)
  158.     Else
  159.         ReDim Preserve oP(0)
  160.     End If
  161.     oP(UBound(oP)) = oPt
  162. End Function
  163.  
  164. Private Function fConvexHull() As tXYZ()
  165. ' Graham-Scan
  166. ' https://www.geeksforgeeks.org/convex-hull-set-2-graham-scan/
  167. ' The algorithm can be adapted to a 3D space, considering these:
  168. ' - the four points A, B, C, D become eight
  169. ' - the quadrilateral becomes a polyhedron with eight vertices
  170. ' - the rectangle becomes a parallelepiped
  171.    
  172.     Dim oP´() As tXYZ
  173.     Dim oP´´() As tXYZ
  174.     Dim oPrun() As tXYZ
  175.     Dim oPrun´() As tXYZ
  176.     Dim oHull() As tXYZ
  177.     Dim lgDecalageRGT As Long
  178.     Dim lgP As Long
  179.     Dim lgPts As Long
  180.     Dim lgTop As Long
  181.     Dim lgPrun´ As Long
  182.     Dim lgPrun´´ As Long
  183.     Dim bPrun As Boolean
  184.     Dim bPrint As Boolean: bPrint = False
  185. bPrint = False
  186.     If fPoints(oP) = 0 Then
  187.     ' Get list of points...
  188.        Call sPointsRnd(100000000, 0, 1000, 0, 1000)
  189. Call fPointsPrint(oP(), 1, 2, 1, -1, -1, bPrint)
  190.     End If
  191.  
  192. Dim dtTimer As Date
  193. dtTimer = VBA.Now()
  194.    
  195.     oPrun() = fPrunePolygon(oP())
  196. Call fPointsPrint(oPrun(), 3, 4, 1, -1, -1, bPrint)
  197.    
  198.     ' Get rectangle R:
  199.    ' X1 = Max(Dx, Ax)
  200.    ' X2 = Min(Bx, Cx)
  201.    ' Y1 = Max(By, Ay)
  202.    ' Y2 = Min(Cy, Dy)
  203.    ' Rectangle R with vertices: (x2, y1), (x2, y2), (x1, y2), (x1, y1)
  204.    
  205.     oPrun´() = oPrun()
  206.     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
  207.    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
  208. Call fPointsPrint(oPrun´(), 5, 6, 1, -1, -1, bPrint)
  209.    
  210.     '#1st pass
  211.    Erase oP´()
  212.     oP´() = oP()
  213.     oP´´() = oP()
  214.     lgPrun´ = LBound(oP´´) - 1
  215.     lgPrun´´ = LBound(oP´´) - 1
  216.     For lgP = LBound(oP) To UBound(oP)
  217.         bPrun = False
  218.         If oPrun´(0).x > oP(lgP).x Then
  219.             lgPrun´´ = lgPrun´´ + 1
  220.             oP´´(lgPrun´´) = oP(lgP)
  221.         ElseIf oPrun´(0).y > oP(lgP).y Then
  222.             lgPrun´´ = lgPrun´´ + 1
  223.             oP´´(lgPrun´´) = oP(lgP)
  224.         ElseIf oPrun´(2).x < oP(lgP).x Then
  225.             lgPrun´´ = lgPrun´´ + 1
  226.             oP´´(lgPrun´´) = oP(lgP)
  227.         ElseIf oPrun´(2).y < oP(lgP).y Then
  228.             lgPrun´´ = lgPrun´´ + 1
  229.             oP´´(lgPrun´´) = oP(lgP)
  230.         Else
  231.             lgPrun´ = lgPrun´ + 1
  232.             oP´(lgPrun´) = oP(lgP)
  233.         End If
  234.     Next lgP
  235.     ReDim Preserve oP´(0 To lgPrun´) ' outside R
  236.    ReDim Preserve oP´´(0 To lgPrun´´) ' inside R
  237.    Call fPointsPrint(oP´(), 7, 8, 1, -1, -1, bPrint)
  238.     Call fPointsPrint(oP´´(), 9, 10, 1, -1, -1, bPrint) '#2nd pass
  239.    lgPrun´´ = LBound(oP´´) - 1
  240.     For lgP = LBound(oP´´) To UBound(oP´´)
  241.     If Not fPointInsidePolygon(oP´´(lgP), oPrun()) Then
  242.         lgPrun´´ = lgPrun´´ + 1
  243.         oP´´(lgPrun´´) = oP´´(lgP)
  244.     End If
  245.     Next lgP
  246.     ReDim Preserve oP´´(0 To lgPrun´´) ' outside Q
  247.    Call fPointsPrint(oP´´(), 11, 12, 1, -1, -1, bPrint)
  248.     Call fPointsSort(oP´´(), False, True)
  249.     ReDim Preserve oHull(LBound(oP´´) To UBound(oP´´)) ' Pick the bottom-most (choose the left-most point in case of tie)
  250.    oHull(LBound(oP)) = oP´´(LBound(oP´´))
  251.     lgPts = LBound(oP´´) ' Right side hull
  252.    
  253.     For lgP = (LBound(oP´´) + 1) To UBound(oP´´)
  254.    
  255.     '!Do While (UBound(oHull) >= 1)
  256.        Do While (lgPts >= 2)
  257.             ' If two or more points make same angle with p0, remove all but the one that is farthest from p0
  258.            ' (in above sorting our criteria was to keep the farthest point at the end)
  259.            '!If CCW(oHull(UBound(oHull) - 1), oHull(UBound(oHull)), oP´´(lgP)) Then Exit Do
  260.            If CCW(oHull(lgPts - 2), oHull(lgPts - 1), oP´´(lgP)) Then Exit Do
  261.            
  262.             'Remove UBound(oHull) point from oHull()
  263.            '!ReDim Preserve oHull(LBound(oHull) To UBound(oHull) - 1)
  264.            lgPts = lgPts - 1
  265.         Loop
  266.        
  267.         'Add Point oP´´(lgP) to oHull()
  268.        '!ReDim Preserve oHull(LBound(oHull) To UBound(oHull) + 1)
  269.        '!oHull(UBound(oHull)) = oP´´(lgP)
  270.        oHull(lgPts) = oP´´(lgP)
  271.         lgPts = lgPts + 1
  272.     Next lgP
  273.    
  274.     ' Left side hull
  275.    lgTop = lgPts
  276.     For lgP = (UBound(oP´´) - 1) To LBound(oP´´) Step -1
  277.         '!Do While (UBound(oHullLFT) >= 1)
  278.        Do While (lgPts > lgTop)
  279.             '!If CCW(oHullLFT(UBound(oHullLFT) - 1), oHullLFT(UBound(oHullLFT)), oP´´(lgP)) Then Exit Do
  280.            If CCW(oHull(lgPts - 2), oHull(lgPts - 1), oP´´(lgP)) Then Exit Do
  281.        
  282.             'Remove UBound(oHull) point from oHull()
  283.            '!ReDim Preserve oHullLFT(LBound(oHullLFT) To UBound(oHullLFT) - 1)
  284.            lgPts = lgPts - 1
  285.         Loop
  286.    
  287.         'Add Point oP´´(lgP) to oHull()
  288.        '!ReDim Preserve oHullLFT(LBound(oHullLFT) To UBound(oHullLFT) + 1)
  289.        '!oHullLFT(UBound(oHullLFT)) = oP´´(lgP)
  290.        oHull(lgPts) = oP´´(lgP)
  291.         lgPts = lgPts + 1
  292.     Next lgP
  293.        
  294.     ReDim Preserve oHull(LBound(oHull) To lgPts - 1)
  295. Call fPointsPrint(oHull(), 13, 14, 1, -1, -1, bPrint)
  296.    
  297.     fConvexHull = oHull()
  298.     'Erase oHull()
  299. Call fPointsPrint(oHull(), 15, 16, 1, -1, -1, bPrint)
  300. Debug.Print dtTimer & " vs " & VBA.Now()
  301. End Function
  302.  
  303. Private Function fPointsPrint(ByRef oP() As tXYZ, _
  304.                               Optional ByVal lC_X As Long = 1, _
  305.                               Optional ByVal lC_Y As Long = 2, _
  306.                               Optional ByVal lR_Start As Long = 1, _
  307.                               Optional ByVal lgP_Start As Long = -1, _
  308.                               Optional ByVal lgP_End As Long = -1, _
  309.                               Optional ByVal bPrint = False) As Boolean
  310.     Dim lgP As Long
  311.     Dim lgR As Long
  312.    
  313.     If Not bPrint Then Exit Function
  314.    
  315.     If lgP_Start < LBound(oP) Then lgP_Start = LBound(oP)
  316.     If lgP_End < LBound(oP) Then lgP_End = UBound(oP)
  317.     lgR = lR_Start
  318.     If lgP_Start <= lgP_End Then
  319.         For lgP = lgP_Start To lgP_End
  320.             Cells(lgR, lC_X).Value2 = oP(lgP).x
  321.             Cells(lgR, lC_Y).Value2 = oP(lgP).y
  322.             lgR = lgR + 1
  323.         Next lgP
  324.     Else 'descending
  325.        For lgP = lgP_Start To lgP_End Step -1
  326.             Cells(lgR, lC_X).Value2 = oP(lgP).x
  327.             Cells(lgR, lC_Y).Value2 = oP(lgP).y
  328.             lgR = lgR + 1
  329.         Next lgP
  330.     End If
  331. End Function
  332.  
  333. Private Function fPointInsidePolygon(ByRef oPoint As tXYZ, ByRef oPolygon() As tXYZ) As Boolean
  334.     Dim lgP As Long
  335.     Dim iInside As Integer
  336.    
  337.     ' Avoid no segment (coincident points)
  338.    For lgP = LBound(oPolygon) To UBound(oPolygon) - 1
  339.         If oPolygon(lgP).x = oPolygon(lgP + 1).x And oPolygon(lgP).y = oPolygon(lgP + 1).y Then
  340.             Call fPointRemove(oPolygon(), lgP)
  341.         End If
  342.     Next lgP
  343.    
  344.     iInside = fLineSide(oPoint, oPolygon(LBound(oPolygon)), oPolygon(LBound(oPolygon) + 1))
  345.     For lgP = LBound(oPolygon) + 1 To UBound(oPolygon) - 1
  346.         If iInside <> fLineSide(oPoint, oPolygon(lgP), oPolygon(lgP + 1)) Then Exit Function
  347.     Next lgP
  348.     If oPolygon(LBound(oPolygon)).x <> oPolygon(UBound(oPolygon)).x Or oPolygon(LBound(oPolygon)).y <> oPolygon(UBound(oPolygon)).y Then
  349.         If iInside <> fLineSide(oPoint, oPolygon(UBound(oPolygon)), oPolygon(LBound(oPolygon))) Then Exit Function
  350.     End If
  351.     fPointInsidePolygon = True
  352. End Function
  353.  
  354. Private Function NewPoint(Optional ByVal x As Double = 0, _
  355.                           Optional ByVal y As Double = 0, _
  356.                           Optional ByVal Z As Double = 0) As tXYZ
  357.     With NewPoint
  358.         .x = x
  359.         .y = y
  360.         .Z = Z
  361.     End With
  362. End Function
  363.  
  364. Private Function fPrunePolygon(ByRef oP() As tXYZ) As tXYZ()
  365. '   A, B, C, D As Points
  366. '   Q As Quadrilateral
  367. '   R As Rectangle
  368. '   DD As Diamond [(Xmid, Ymin)
  369. '
  370. '   for each Point p in oP
  371. '         Update A, B, C, D --> Q
  372. '   # First pass
  373. '   Create R from Q
  374. '   for each Point p in S
  375. '      if p inside R
  376. '         prune p from S
  377. '   # Second pass
  378. '   for each Point p in S
  379. '      if p inside Q
  380. '         prune p from S
  381.  
  382. '   ^         D    / +\
  383. '   |           /....  + \
  384. '   |        /  / Q x....   \ C
  385. '   |     /    /---------....| \
  386. '   |  /      /|   *     *   |    \
  387. '   |  \     /x|*    R  *  * | +  /
  388. '   |     \ /..|.....--------|  /
  389. '   |      A \        ..x...|/ B
  390. '   |           \ +  +   /
  391. '   |              \  /
  392. '   --------------------------->
  393. '
  394. '   * points inside R rectangle are prunable
  395. '   x point inside Q are prunable but not until the end
  396. '   + point outside Q are not prunable
  397. '
  398.    If fPoints(oP) = 0 Then
  399.     ' Get list of points...
  400.        Call sPointsRnd(300, 0, 100, 0, 100)
  401.     End If
  402.    
  403.     Dim PtA As tXYZ, PtB As tXYZ, PtC As tXYZ, PtD As tXYZ
  404.     Dim Q(0 To 3) As tXYZ 'Quadrilateral
  405.    Dim R(0 To 3) As tXYZ 'Rectangle
  406.    Dim DD(0 To 3) As tXYZ  'Diamond
  407.  
  408.     Dim lgP As Long
  409.  
  410.     Q(0) = oP(LBound(oP)): Q(1) = Q(0): Q(2) = Q(0): Q(3) = Q(0)
  411.     For lgP = LBound(oP) To UBound(oP)
  412.         With Q(0) ' Point A
  413.            If (-.x - .y) < (-oP(lgP).x - oP(lgP).y) Then Q(0) = oP(lgP)
  414.         End With
  415.         With Q(1) ' Point B
  416.            If (.x - .y) < (oP(lgP).x - oP(lgP).y) Then Q(1) = oP(lgP)
  417.         End With
  418.         With Q(2) ' Point C
  419.            If (.x + .y) < (oP(lgP).x + oP(lgP).y) Then Q(2) = oP(lgP)
  420.         End With
  421.         With Q(3) ' Point D
  422.            If (-.x + .y) < (-oP(lgP).x + oP(lgP).y) Then Q(3) = oP(lgP)
  423.         End With
  424.     Next lgP
  425.  
  426.     fPrunePolygon = Q()
  427. End Function
  428.  
  429. Private Sub sPointsRnd(Optional ByVal lgPts As Long = 0, _
  430.                        Optional ByVal Xmin As Double = 0, _
  431.                        Optional ByVal Xmax As Double = 0, _
  432.                        Optional ByVal Ymin As Double = 0, _
  433.                        Optional ByVal Ymax As Double = 0)
  434.     Dim lgP As Long
  435.     Dim dbTmp As Long
  436.  
  437.     If Xmin = 0 Then Xmin = (Rnd() * 100) + 0
  438.     If Xmax = 0 Then Xmax = (Rnd() * 100) + 0
  439.     If Ymin = 0 Then Ymin = (Rnd() * 100) + 0
  440.     If Ymax = 0 Then Ymax = (Rnd() * 100) + 0
  441.     If Xmax < Xmin Then
  442.         Xmax = dbTmp
  443.         Xmax = Xmin
  444.         Xmin = dbTmp
  445.     End If
  446.     If Ymax < Ymin Then
  447.         Ymax = dbTmp
  448.         Ymax = Ymin
  449.         Ymin = dbTmp
  450.     End If
  451.     ReDim Preserve oP(0 To lgPts - 1)
  452.     For lgP = LBound(oP) To UBound(oP)
  453.         oP(lgP) = NewPoint((Rnd() * Xmax - Xmin) + Xmin, (Rnd() * Ymax - Ymin) + Ymin)
  454.     Next lgP
  455. End Sub
  456.  
  457. Private Sub sDelaunay()
  458.     Dim oChrt As Excel.ChartObject
  459.     Dim rgData As Excel.Range
  460.     Dim aData As Variant
  461.     Dim lgP As Long
  462.     Dim lgT As Long
  463.     Dim lgTs As Long ' total number of triangles
  464.    Dim Xmin As Double, Ymin As Double
  465.     Dim Xmax As Double, Ymax As Double
  466.     Dim Xmid As Double, Ymid As Double
  467.     Dim HSide As Double, VSide As Double
  468.     Dim i12 As Integer, i23 As Integer, i31 As Integer
  469.     Dim bDalaunay As Boolean
  470.    
  471.     ' Elements are sorted by X
  472.    With ActiveSheet
  473.         Call sPointsCreate(10)
  474.         Set rgData = .Range("A1", .Range("B1", .Range("B1").End(xlDown)))
  475.     End With
  476.     aData = rgData.Value2
  477.     'ReDim oP(0 To 2 + (UBound(aData, 1) - LBound(aData, 1) + 1))
  478.    ReDim aStack(LBound(oP) To UBound(oP)) As Long
  479.     ReDim aBlnk(LBound(oP) To UBound(oP)) As Long
  480.    
  481.     For lgP = LBound(oP) + 3 To UBound(oP)
  482.         With oP(lgP)
  483.             '.X = aData(lgP - 2, 1)
  484.            '.Y = aData(lgP - 2, 2)
  485.            If .x < Xmin Then Xmin = .x
  486.             If .x > Xmax Then Xmax = .x
  487.             If .y < Ymin Then Ymin = .y
  488.             If .y > Ymax Then Ymax = .y
  489.         End With
  490.     Next lgP
  491.    
  492.     lgT = -1
  493.    
  494.     ' Generate the super-triangle
  495.    Xmid = 0.5 * (Xmax + Xmin)
  496.     Ymid = 0.5 * (Ymax + Ymin)
  497.     HSide = 2 * (Xmid - Xmin)
  498.     VSide = 2 * (Ymid - Ymin)
  499.     With oP(LBound(oP) + 0)
  500.         .x = Xmid - (3 * HSide)
  501.         .y = Ymid - (3 * VSide)
  502.     End With
  503.     With oP(LBound(oP) + 1)
  504.         .x = Xmid
  505.         .y = Ymid + (3 * VSide)
  506.     End With
  507.     With oP(LBound(oP) + 2)
  508.         .x = Xmid + (3 * HSide)
  509.         .y = Ymid - (3 * VSide)
  510.     End With
  511.    
  512.     Set oChrt = ActiveSheet.ChartObjects("ChrtPts")
  513.     Call fChrtSeriesDelete(oChrt)
  514.  
  515.     ' Add supertriangle points
  516.    Call fChrtTriangleAdd(oP(LBound(oP) + 0), _
  517.                           oP(LBound(oP) + 1), _
  518.                           oP(LBound(oP) + 2), _
  519.                           oChrt)
  520.    
  521.     lgT = lgT + 1
  522.     lgTs = lgT
  523.     ReDim Preserve oT(0 To lgT)
  524.     With oT(lgT)
  525.         .T1 = -1
  526.         .T2 = -2
  527.         .T3 = -3
  528.    
  529.         .V1 = LBound(oP) + 0
  530.         .V2 = LBound(oP) + 1
  531.         .V3 = LBound(oP) + 2
  532.     End With
  533.    
  534.     ' Get new max-min coordinates...
  535. '    Call fTriangleParameters(lgT)
  536.  
  537.     ' Generate new triangles:
  538.    For lgP = LBound(oP) + 3 To UBound(oP)
  539.         Call fChrtPtHighlight(oChrt, 1, lgP - 3 + 1, True)
  540.         'Call fChrtPtAdd(oP(lgP), oChrt)
  541.  
  542.         For lgT = LBound(oT) To UBound(oT)
  543.             ' Discard point inside boundary
  544.            If oP(lgP).x < oT(lgT).Xmax Then
  545.             If oP(lgP).x > oT(lgT).Xmin Then
  546.             If oP(lgP).y < oT(lgT).Ymax Then
  547.             If oP(lgP).y > oT(lgT).Ymin Then
  548.                 ' Point is inside boundary of triangle, so check more intensively
  549.                i12 = fLineSide(oP(lgP), oP(oT(lgT).V2), oP(oT(lgT).V1))
  550.                 i23 = fLineSide(oP(lgP), oP(oT(lgT).V3), oP(oT(lgT).V2))
  551.                 i31 = fLineSide(oP(lgP), oP(oT(lgT).V1), oP(oT(lgT).V3))
  552.                
  553.                 If (i12 = i23 And i23 = i31) Then
  554.                     ReDim Preserve oT(0 To lgTs + 2)
  555.                     ' Always store points in Counter-clockwise order (last point is new point)
  556.                    With oT(lgTs + 1)
  557.                         .T1 = oT(lgT).T2
  558.                         .T2 = lgTs + 2
  559.                         .T3 = lgT
  560.                    
  561.                         .V1 = oT(lgT).V2
  562.                         .V2 = oT(lgT).V3
  563.                         .V3 = lgP
  564.                    
  565.                         ' Add supertriangle points
  566.                        Call fChrtTriangleAdd(oP(.V1), _
  567.                                               oP(.V2), _
  568.                                               oP(.V3), _
  569.                                               ActiveSheet.ChartObjects("ChrtPts"))
  570.                     End With
  571.                    
  572.                     With oT(lgTs + 2)
  573.                         .T1 = oT(lgT).T3
  574.                         .T2 = lgT
  575.                         .T3 = lgTs + 1
  576.                    
  577.                         .V1 = oT(lgT).V3
  578.                         .V2 = oT(lgT).V1
  579.                         .V3 = lgP
  580.                        
  581.                         ' Add supertriangle points
  582.                        Call fChrtTriangleAdd(oP(.V1), _
  583.                                               oP(.V2), _
  584.                                               oP(.V3), _
  585.                                               ActiveSheet.ChartObjects("ChrtPts"))
  586.                     End With
  587.                    
  588.                     ' Reformat initial triangle
  589.                    With oT(lgT)
  590.                         '.T1 = oT(lgT).T1 ' this neighbour is not rewritten
  591.                        .T2 = lgTs + 1
  592.                         .T3 = lgTs + 2
  593.                    
  594.                         '.V1 = .V1 ' this vertex is not rewritten
  595.                        '.V2 = .V2 ' this vertex is not rewritten
  596.                        .V3 = lgP
  597.                     End With
  598.                    
  599.                     ' Get new max-min coordinates for final triangles...
  600. '                    Call fTriangleParameters(lgT)
  601. '                    Call fTriangleParameters(lgTs + 1)
  602. '                    Call fTriangleParameters(lgTs + 2)
  603.                    
  604.                     lgTs = lgTs + 2
  605.    
  606.                     ' Is the optimum delaunay?
  607.                    ' Each new triangle has three neighbours that have to be checked for Delaunay condition
  608.                    ' For Neighbour triangle 1
  609.                    'call fDelaunay(lgT, oT(lgT).T1)
  610.                    
  611.                     ' For Neighbour triangle 2
  612.                    'call fDelaunay(lgT, oT(lgT).T2)
  613.                    
  614.                     ' For Neighbour triangle 3
  615.                    'Call fDelaunay(lgT, oT(lgT).T3)
  616.                    
  617.                     Exit For ' Goto NextTriangle
  618.                End If
  619.             End If
  620.             End If
  621.             End If
  622.             End If
  623.         Next lgT
  624.        
  625.         Call fChrtPtHighlight(oChrt, 1, lgP - 3, False)
  626.     Next lgP
  627.  
  628.     ' Delete triangles on the supertriangle structure
  629.    Dim oT_() As tTriangle: oT_() = oT()
  630.     For lgT = LBound(oT) To UBound(oT)
  631.         With oT(lgT)
  632.             If .V1 >= 0 And _
  633.                .V2 >= 0 And _
  634.                .V3 >= 0 Then
  635.                'lgT_ = lgT_ + 1
  636.               'oT_(0 to lgT) = oT(lgT)
  637.            End If
  638.         End With
  639.     Next lgT
  640.     'Redim preserve oT_(0 to lgT_)
  641.    'oT() = oT_()
  642.    'Erase oT_()
  643. Stop
  644. End Sub
  645.  
  646. Private Function fDelaunay(ByVal lgT1 As Long, _
  647.                            ByVal lgT2 As Long) As Boolean
  648.     Dim bSwap            As Boolean
  649.     Dim TPtrTmp          As Long
  650.  
  651.     ' Rotate vertices of both triangles so in both triangles vertices 3 are opposed
  652.    If oT(lgT1).V1 = oT(lgT2).V2 And oT(lgT1).V2 = oT(lgT2).V1 Then
  653.         ' Do nothing...
  654.    ElseIf oT(lgT1).V1 = oT(lgT2).V1 And oT(lgT1).V2 = oT(lgT2).V3 Then
  655.         ' rotate 2 clockwise
  656.        Call fTriangleRotate(lgT2, False)
  657.     ElseIf oT(lgT1).V1 = oT(lgT2).V3 And oT(lgT1).V2 = oT(lgT2).V2 Then
  658.         ' rotate 2 counter-clockwise
  659.        Call fTriangleRotate(lgT2, True)
  660.         '--
  661.    ElseIf oT(lgT1).V2 = oT(lgT2).V1 And oT(lgT1).V3 = oT(lgT2).V3 Then
  662.         ' rotate 1 counter-clockwise
  663.        Call fTriangleRotate(lgT1, True)
  664.         ' rotate 2 clockwise
  665.        Call fTriangleRotate(lgT2, False)
  666.     ElseIf oT(lgT1).V2 = oT(lgT2).V2 And oT(lgT1).V3 = oT(lgT2).V1 Then
  667.         ' rotate 1 counter-clockwise
  668.        Call fTriangleRotate(lgT1, True)
  669.     ElseIf oT(lgT1).V2 = oT(lgT2).V3 And oT(lgT1).V3 = oT(lgT2).V2 Then
  670.         ' rotate 1 counter-clockwise
  671.        Call fTriangleRotate(lgT1, True)
  672.         ' rotate 2 counter-clockwise
  673.        Call fTriangleRotate(lgT2, True)
  674.         '--
  675.    ElseIf oT(lgT1).V3 = oT(lgT2).V1 And oT(lgT1).V1 = oT(lgT2).V3 Then
  676.         ' rotate 1 clockwise
  677.        Call fTriangleRotate(lgT1, False)
  678.         ' rotate 2 clockwise
  679.        Call fTriangleRotate(lgT2, False)
  680.     ElseIf oT(lgT1).V3 = oT(lgT2).V2 And oT(lgT1).V1 = oT(lgT2).V1 Then
  681.         ' rotate 1 clockwise
  682.        Call fTriangleRotate(lgT1, False)
  683.     ElseIf oT(lgT1).V3 = oT(lgT2).V3 And oT(lgT1).V1 = oT(lgT2).V2 Then
  684.         ' rotate 1 clockwise
  685.        Call fTriangleRotate(lgT1, False)
  686.         ' rotate 2 counter-clockwise
  687.        Call fTriangleRotate(lgT2, True)
  688.     End If
  689.  
  690.     If fDistance2DPoints(oP(oT(lgT2).V2), oT(lgT1).Center) < EPSILON Then
  691.         bSwap = True
  692.         TPtrTmp = oT(lgT1).T2    'Destroy oT(lgT1).T1 and oT(lgT2).T1
  693.        oT(lgT1).T1 = oT(lgT2).T2
  694.         oT(lgT1).T2 = lgT2 'oT(lgT1).T3 = does not change
  695.        oT(lgT2).T1 = TPtrTmp
  696.         oT(lgT2).T2 = lgT1 'oT(lgT2).T3 = does not change
  697.        ' Swap vertices
  698.        oT(lgT1).V2 = oT(lgT2).V3
  699.         oT(lgT2).V2 = oT(lgT1).V3
  700.     End If
  701.     If bSwap Then
  702.         Call fTriangleParameters(lgT1)
  703.         Call fTriangleParameters(lgT2) '
  704.        'Call fDelaunay(lgT1, oT(lgT1).T1) Then
  705.        'Call fDelaunay(lgT1, oT(lgT1).T2) Then
  706.        'Call fDelaunay(lgT1, oT(lgT1).T3) Then '
  707.        'Call fDelaunay(lgT2, oT(lgT2).T1) Then
  708.        'Call fDelaunay(lgT2, oT(lgT2).T2) Then
  709.        'Call fDelaunay(lgT2, oT(lgT2).T3) Then
  710.    Else
  711.         fDelaunay = True
  712.     End If
  713. End Function
  714.        
  715. Private Function fTriangleRotate(ByRef lgT As Long, Optional ByVal bCCW As Boolean = True) As Boolean
  716. ' Given a triangle, rotate vertices oTTmp = oT(lgT)
  717.        
  718.     With oT(lgT)
  719.         If bCCW Then ' counter-clockwise
  720.            .V1 = oTTmp.V2
  721.             .V2 = oTTmp.V3
  722.             .V3 = oTTmp.V1
  723.             .T1 = oTTmp.T2
  724.             .T2 = oTTmp.T3
  725.             .T3 = oTTmp.T1
  726.         Else
  727.             .V1 = oTTmp.V3
  728.             .V2 = oTTmp.V1
  729.             .V3 = oTTmp.V2
  730.             .T1 = oTTmp.T3
  731.             .T2 = oTTmp.T1
  732.             .T3 = oTTmp.T2
  733.         End If
  734.     End With
  735. End Function
  736.  
  737. Private Function fDelaunayDiagonal(ByRef lgT1 As Long, ByRef lgT2 As Long) As Boolean
  738. ' Given two triangles, find the opposed vertices
  739. ' Rotate triangles so the 3 vertices are opposed...
  740.    Dim dX1 As Double, dY1 As Double
  741.     Dim dX2 As Double, dY2 As Double
  742.     dX1 = oP(oT(lgT1).V3).x - oP(oT(lgT2).V3).x
  743.     dY1 = oP(oT(lgT1).V3).y - oP(oT(lgT2).V3).y
  744.     dX2 = oP(oT(lgT1).V1).x - oP(oT(lgT2).V2).x
  745.     dY2 = oP(oT(lgT1).V1).y - oP(oT(lgT2).V2).y
  746.     If ((dX1 * dX1) + (dY1 * dY1)) > ((dX2 * dX2) + (dY2 * dY2)) Then
  747.         Call fSwapDiagonal(lgT1, lgT2)
  748.     End If
  749. '    End If
  750. End Function
  751.  
  752. Private Function fSwapDiagonal(ByRef lgT1 As Long, _
  753.                                ByRef lgT2 As Long) As Boolean
  754. ' For given diagonal D1 on triangle T1 = diagonal D2 on triangle T2, will swap diagonals with opposed vertices
  755. ' Rotate triangles so the 3 vertices are opposed...
  756.    Dim TPtrTmp As Long
  757.  
  758.     With oT(lgT1)
  759.         TPtrTmp = .T2
  760.        
  761.         'Destroy oT(lgT1).T1 and oT(lgT2).T1
  762.        .T1 = oT(lgT2).T2
  763.         .T2 = lgT2
  764.         '.T3 = does not change
  765.    End With
  766.    
  767.     With oT(lgT2)
  768.         .T1 = TPtrTmp
  769.         .T2 = lgT1
  770.         '.T3 = does not change
  771.    End With
  772.    
  773.     ' Swap vertices
  774.    oT(lgT1).V2 = oT(lgT2).V3
  775.     oT(lgT2).V2 = oT(lgT1).V3
  776.  
  777. End Function
  778.  
  779. Private Function fNeigbourSide(ByRef iSide As Integer, _
  780.                                ByRef lgT1 As Long, _
  781.                                ByRef lgT2 As Long) As Integer
  782. ' For given iSide on triangle1, return neighbour side on triangle2.
  783. ' Both triangles turn the same: 1->2->3
  784.    
  785.     If iSide = 1 Then
  786.         If oT(lgT1).V1 = oT(lgT2).V1 Then
  787. '            fNeigbourSide = 3
  788.        ElseIf oT(lgT1).V1 = oT(lgT2).V2 Then
  789. '            fNeigbourSide = 2
  790.        ElseIf oT(lgT1).V1 = oT(lgT2).V3 Then
  791. '            fNeigbourSide = 1
  792.        End If
  793.     ElseIf iSide = 2 Then
  794.         If oT(lgT1).V2 = oT(lgT2).V1 Then
  795. '            fNeigbourSide = 3
  796.        ElseIf oT(lgT1).V2 = oT(lgT2).V2 Then
  797. '            fNeigbourSide = 2
  798.        ElseIf oT(lgT1).V2 = oT(lgT2).V3 Then
  799. '            fNeigbourSide = 1
  800.        End If
  801.     ElseIf iSide = 3 Then
  802.         If oT(lgT1).V3 = oT(lgT2).V1 Then
  803. '            fNeigbourSide = 3
  804.        ElseIf oT(lgT1).V3 = oT(lgT2).V2 Then
  805. '            fNeigbourSide = 2
  806.        ElseIf oT(lgT1).V3 = oT(lgT2).V3 Then
  807. '            fNeigbourSide = 1
  808.        End If
  809.     End If
  810. End Function
  811.  
  812. Private Function fTriangleFromPoints(ByRef oP1 As tXYZ, _
  813.                                      ByRef oP2 As tXYZ, _
  814.                                      ByRef oP3 As tXYZ) As tTriangle
  815.     With fTriangleFromPoints
  816.         .V1 = 1
  817.         .V2 = 2
  818.         .V3 = 3
  819.        
  820.         .T1 = -1
  821.         .T2 = -2
  822.         .T3 = -3
  823.        
  824.         .Xmin = oP1.x
  825.         .Xmax = oP1.x
  826.         .Ymin = oP1.y
  827.         .Ymax = oP1.y
  828.    
  829.         If .Xmin > oP2.x Then .Xmin = oP2.x
  830.         If .Xmin > oP3.x Then .Xmin = oP3.x
  831.         If .Xmax < oP2.x Then .Xmax = oP2.x
  832.         If .Xmax < oP3.x Then .Xmax = oP3.x
  833.         If .Ymin > oP2.y Then .Ymin = oP2.y
  834.         If .Ymin > oP3.y Then .Ymin = oP3.y
  835.         If .Ymax < oP2.y Then .Ymax = oP2.y
  836.         If .Ymax < oP3.y Then .Ymax = oP3.y
  837.         Dim B As tXYZ
  838.         Dim C As tXYZ
  839.         Dim D As Double
  840.         Dim Bx² As Double
  841.         Dim By² As Double
  842.         Dim Cx² As Double
  843.         Dim Cy² As Double
  844.        
  845.         B.x = oP2.x - oP1.x
  846.         B.y = oP2.y - oP1.y
  847.         C.x = oP3.x - oP1.x
  848.         C.y = oP3.y - oP1.y
  849.         Bx² = B.x * B.x: By² = B.y * B.y
  850.         Cx² = C.x * C.x: Cy² = C.y * C.y
  851.         D = 1 / (2 * (B.x * C.y - B.y * C.x))
  852.         .Center.x = D * (C.y * (Bx² + By²) - B.y * (Cx² + Cy²))
  853.         .Center.y = D * (B.x * (Cx² + Cy²) - C.x * (Bx² + By²))
  854.         .R² = (.Center.x * .Center.x) + (.Center.y * .Center.y)
  855.         .Center.x = .Center.x + oP1.x
  856.         .Center.y = .Center.y + oP1.y
  857.     End With
  858.    
  859. End Function
  860.  
  861. Private Function fTriangleParameters(ByRef lgT As Long) As Boolean
  862.     With oT(lgT)
  863.         .Xmin = oP(.V1).x
  864.         .Xmax = oP(.V1).x
  865.         .Ymin = oP(.V1).y
  866.         .Ymax = oP(.V1).y
  867.         If .Xmin > oP(.V2).x Then .Xmin = oP(.V2).x
  868.         If .Xmin > oP(.V3).x Then .Xmin = oP(.V3).x
  869.         If .Xmax < oP(.V2).x Then .Xmax = oP(.V2).x
  870.         If .Xmax < oP(.V3).x Then .Xmax = oP(.V3).x
  871.         If .Ymin > oP(.V2).y Then .Ymin = oP(.V2).y
  872.         If .Ymin > oP(.V3).y Then .Ymin = oP(.V3).y
  873.         If .Ymax < oP(.V2).y Then .Ymax = oP(.V2).y
  874.         If .Ymax < oP(.V3).y Then .Ymax = oP(.V3).y
  875.        
  876.         .Center = Circumcenter(oP(.V1).x, oP(.V1).y, oP(.V2).x, oP(.V2).y, oP(.V3).x, oP(.V3).y)
  877.         .R² = (.Center.x - oP(.V1).x) ^ 2 + (.Center.y - oP(.V1).y) ^ 2
  878.        
  879.         'Dim B As tXYZ
  880.        'Dim C As tXYZ
  881.        'Dim D As Double
  882.        'Dim Bx² As Double
  883.        'Dim By² As Double
  884.        'Dim Cx² As Double
  885.        'Dim Cy² As Double
  886.        '
  887.        'B.X = oP(.V2).X - oP(.V1).X
  888.        'B.Y = oP(.V2).Y - oP(.V1).Y
  889.        'C.X = oP(.V3).X - oP(.V1).X
  890.        'C.Y = oP(.V3).Y - oP(.V1).Y
  891.        'Bx² = B.X * B.X: By² = B.Y * B.Y
  892.        'Cx² = C.X * C.X: Cy² = C.Y * C.Y
  893.        'D = 1 / (2 * (B.X * C.Y - B.Y * C.X))
  894.        '.Center.X = D * (C.Y * (Bx² + By²) - B.Y * (Cx² + Cy²))
  895.        '.Center.Y = D * (B.X * (Cx² + Cy²) - C.X * (Bx² + By²))
  896.        '.R² = (.Center.X * .Center.X) + (.Center.Y * .Center.Y)
  897.        '.Center.X = .Center.X + oP(.V1).X
  898.        '.Center.Y = .Center.Y + oP(.V1).Y
  899.    End With
  900. End Function
  901.  
  902. Private Function Circumcenter(ByVal x1 As Double, ByVal y1 As Double, _
  903.                               ByVal x2 As Double, ByVal y2 As Double, _
  904.                               ByVal x3 As Double, ByVal y3 As Double) As tXYZ
  905.     Dim A As Double
  906.     Dim B As Double
  907.     Dim C As Double
  908.     Dim D As Double
  909.  
  910.     A = x1 * x1 + y1 * y1
  911.     B = x2 * x2 + y2 * y2
  912.     C = x3 * x3 + y3 * y3
  913.     D = 2 * (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2))
  914.    
  915.     If D <> 0 Then
  916.         With Circumcenter
  917.             .x = (A * (y2 - y3) + B * (y3 - y1) + C * (y1 - y2)) / D
  918.             .y = (A * (x3 - x2) + B * (x1 - x3) + C * (x2 - x1)) / D
  919.         End With
  920.     End If
  921. End Function
  922.  
  923. Private Sub sLineSide()
  924.     Dim oPoint As tXYZ, oPt1 As tXYZ, oPt2 As tXYZ
  925.  
  926.     With oPoint
  927.         .x = 1
  928.         .y = -1
  929.     End With
  930.     With oPt1
  931.         .x = 0
  932.         .y = 0
  933.     End With
  934.     With oPt2
  935.         .x = 10
  936.         .y = 0
  937.     End With
  938.     Debug.Print fLineSide(oPoint, oPt1, oPt2)
  939. End Sub
  940.  
  941. Private Function CCW(ByRef oPt1 As tXYZ, _
  942.                      ByRef oPt2 As tXYZ, _
  943.                      ByRef oPt3 As tXYZ) As Boolean
  944. ' If counter clock-wise, then CCW is true
  945.    CCW = ((oPt2.x - oPt1.x) * (oPt3.y - oPt1.y)) > ((oPt2.y - oPt1.y) * (oPt3.x - oPt1.x))
  946. End Function
  947.  
  948. Private Function fLineSide(ByRef oPoint As tXYZ, _
  949.                            ByRef oPt1 As tXYZ, _
  950.                            ByRef oPt2 As tXYZ) As Integer
  951. ' Use the sign of the determinant of vectors (AB,AM), where M(X,Y) is the query point:
  952. ' Position = Sign((Bx - Ax) * (Y - Ay) - (By - Ay) * (X - Ax))
  953. ' It is 0 on the line, and -1 on right side, +1 on the left side.
  954.    
  955.     fLineSide = VBA.Sgn((oPt2.x - oPt1.x) * (oPoint.y - oPt1.y) - (oPt2.y - oPt1.y) * (oPoint.x - oPt1.x))
  956.  
  957. End Function
  958.  
  959. Private Function fChrtSeriesDelete(ByVal oChrt As Excel.ChartObject)
  960.     Dim lgSeries As Long
  961.    
  962.     With oChrt
  963.         With .Chart.SeriesCollection
  964.             For lgSeries = .Count To 2 Step -1
  965.                 oChrt.Chart.SeriesCollection(lgSeries).Delete
  966.             Next lgSeries
  967.         End With
  968.     End With
  969. End Function
  970.  
  971. Private Function fChrtCircleAdd(ByRef oCenter As tXYZ, _
  972.                                 ByRef Radius As Double, _
  973.                                 Optional ByVal oChrt As Excel.ChartObject) As Boolean
  974.     Dim oSeries As Excel.Series
  975.     Dim lgAngle As Long
  976.     Dim dbAngleRAD As Double
  977.     Dim strX As String
  978.     Dim strY As String
  979.    
  980.     With oChrt
  981.         With .Chart
  982.             Set oSeries = .SeriesCollection.NewSeries
  983.             With oSeries
  984.                 For lgAngle = 0 To 17
  985.                     dbAngleRAD = lgAngle * PI_8th
  986.                     strX = strX & Replace(oCenter.x + Radius * Cos(dbAngleRAD), ",", ".") & ","
  987.                     strY = strY & Replace(oCenter.y + Radius * Sin(dbAngleRAD), ",", ".") & ","
  988.                 Next lgAngle
  989.                 strX = VBA.Left$(strX, Len(strX) - 1)
  990.                 strY = VBA.Left$(strY, Len(strY) - 1)
  991.                 .XValues = "={" & strX & "}"
  992.                 .Values = "={" & strY & "}"
  993.                 With .Format.Line
  994.                     .Visible = msoTrue
  995.                     '.ForeColor.ObjectThemeColor = msoThemeColorAccent1
  996.                    '.ForeColor.TintAndShade = 0
  997.                    '.ForeColor.Brightness = 0
  998.                    .Weight = 0.25
  999.                     .Visible = msoTrue
  1000.                     .DashStyle = msoLineSysDash
  1001.                 End With
  1002.             End With
  1003.         End With
  1004.     End With
  1005. End Function
  1006.  
  1007. Private Function fChrtTriangleAdd(ByRef oPt1 As tXYZ, _
  1008.                                   ByRef oPt2 As tXYZ, _
  1009.                                   ByRef oPt3 As tXYZ, _
  1010.                                   Optional ByVal oChrt As Excel.ChartObject) As Boolean
  1011.     Dim oSeries As Excel.Series
  1012.    
  1013.     With oChrt
  1014.         With .Chart
  1015.             Set oSeries = .SeriesCollection.NewSeries
  1016.             With oSeries
  1017.                 .XValues = "={" & Replace(oPt1.x, ",", ".") & "," & Replace(oPt2.x, ",", ".") & "," & Replace(oPt3.x, ",", ".") & "," & Replace(oPt1.x, ",", ".") & "}"
  1018.                 .Values = "={" & Replace(oPt1.y, ",", ".") & "," & Replace(oPt2.y, ",", ".") & "," & Replace(oPt3.y, ",", ".") & "," & Replace(oPt1.y, ",", ".") & "}"
  1019.                 With .Format.Line
  1020.                     .Visible = msoTrue
  1021.                     .Weight = 0.25
  1022.                     .DashStyle = msoLineSysDash
  1023.                     '.ForeColor.ObjectThemeColor = msoThemeColorAccent1
  1024.                    '.ForeColor.TintAndShade = 0
  1025.                    '.ForeColor.Brightness = 0
  1026.                End With
  1027.             End With
  1028.         End With
  1029.     End With
  1030. End Function
  1031.  
  1032. Private Function fChrtPtAdd(ByRef oPt As tXYZ, _
  1033.                             Optional ByVal oChrt As Excel.ChartObject)
  1034.     Dim oSeries As Excel.Series
  1035.    
  1036.     With oChrt
  1037.         With .Chart
  1038.             Set oSeries = .SeriesCollection.NewSeries
  1039.             With oSeries
  1040.                 .XValues = "={" & Replace(oPt.x, ",", ".") & "}"
  1041.                 .Values = "={" & Replace(oPt.y, ",", ".") & "}"
  1042.                 With .Format.Fill
  1043.                     .Visible = msoTrue
  1044.                     .ForeColor.RGB = RGB(255, 0, 0)
  1045.                     '.ForeColor.ObjectThemeColor = msoThemeColorAccent1
  1046.                    '.ForeColor.TintAndShade = 0
  1047.                    '.ForeColor.Brightness = 0
  1048.                End With
  1049.             End With
  1050.         End With
  1051.     End With
  1052. End Function
  1053.  
  1054. Private Function fChrtPtHighlight(ByVal oChrt As Excel.ChartObject, _
  1055.                                   ByVal lgSeries As Long, _
  1056.                                   ByVal lgData As Long, _
  1057.                                   Optional ByVal bActive As Boolean = False)
  1058.     With oChrt
  1059.         With .Chart.FullSeriesCollection(lgSeries).Points(lgData)
  1060.             With .Format
  1061.                 With .Fill
  1062.                     .Visible = msoTrue
  1063.                     If bActive Then
  1064.                         .ForeColor.RGB = RGB(255, 0, 0)
  1065.                     Else
  1066.                         .ForeColor.ObjectThemeColor = msoThemeColorAccent1
  1067.                         '.ForeColor.TintAndShade = 0
  1068.                        '.ForeColor.Brightness = 0
  1069.                    End If
  1070.                     .Transparency = 0
  1071.                     .Solid
  1072.                 End With
  1073.             End With
  1074.            
  1075.             '.ApplyDataLabels
  1076.        End With
  1077.     End With
  1078. End Function
  1079.  
  1080. Private Sub sPointsCreate(ByVal lgPoints As Long)
  1081.     Dim Xmin As Double, Ymin As Double
  1082.     Dim Xmax As Double, Ymax As Double
  1083.     Dim lgPoint As Long
  1084.     Dim dX As Double, dY As Double
  1085.     Dim lgR As Long
  1086.     Dim rgData As Excel.Range
  1087.     Dim XValue As Excel.Range
  1088.     Dim YValue As Excel.Range
  1089.     Dim oChrt As Excel.ChartObject
  1090.  
  1091.     ' Create points and print out to range
  1092.    With ActiveSheet
  1093.         Xmin = 0
  1094.         Xmax = 1000
  1095.         Ymin = 0
  1096.         Ymax = 1000
  1097.         dX = (Xmax - Xmin)
  1098.         dY = (Ymax - Ymin)
  1099.         ReDim oP(0 To 2 + lgPoints)
  1100.        
  1101.         For lgPoint = 3 To (lgPoints - 1) + 3
  1102.             oP(lgPoint).x = Xmin + (Rnd() * dX)
  1103.             oP(lgPoint).y = Ymin + (Rnd() * dY)
  1104.             lgR = lgPoint - 1 + 3
  1105.             .Cells(lgR, 1).Value2 = oP(lgPoint).x
  1106.             .Cells(lgR, 2).Value2 = oP(lgPoint).y
  1107.         Next lgPoint
  1108.        
  1109.         Set XValue = .Range("A1", .Range("A1", .Range("A1").End(xlDown)))
  1110.         'XValue.Select
  1111.        Set YValue = .Range("B1", .Range("B1", .Range("B1").End(xlDown)))
  1112.         'YValue.Select
  1113.        Set rgData = Union(XValue, YValue)
  1114.    
  1115.         ' Sort points by X
  1116.        Call fPointsSort(oP(), True, True)
  1117.         With .Sort
  1118.         '    .SortFields.Clear
  1119.        '    .SortFields.Add _
  1120.                 Key:=rgData, _
  1121.                 SortOn:=xlSortOnValues, _
  1122.                 Order:=xlAscending, _
  1123.                 DataOption:=xlSortNormal
  1124.            
  1125.         '    .SetRange rgData
  1126.        '    .Header = xlGuess
  1127.        '    .MatchCase = False
  1128.        '    .Orientation = xlTopToBottom
  1129.        '    .SortMethod = xlPinYin
  1130.        '    '.Apply
  1131.        End With
  1132.        
  1133.         ' Create chart
  1134.        'Set oChrt = .Shapes.AddChart2(240, xlXYScatter).ChartObject
  1135.        With .Shapes.AddChart2(240, xlXYScatter)
  1136.             .Name = "ChrtPts"
  1137.             .Left = ActiveSheet.Range("C1").Left
  1138.             .Top = ActiveSheet.Range("C1").Top
  1139.            
  1140.             .Height = 800
  1141.             .Width = 800
  1142.             With .Chart
  1143.                 With .PlotArea
  1144.                     '.Height = 700
  1145.                    '.Width = 700
  1146.                End With
  1147.  
  1148.                 .SetSourceData Source:=rgData
  1149.                 .ChartTitle.Delete
  1150.                 With .Axes(xlValue)
  1151.                     '.MinimumScaleIsAuto = True
  1152.                    '.MaximumScaleIsAuto = True
  1153.                    .MinimumScale = 0
  1154.                     .MaximumScale = 1000
  1155.                     .MajorUnit = 250
  1156.                     .MinorUnit = 50
  1157.                 End With
  1158.                 With .Axes(xlCategory)
  1159.                     '.MinimumScaleIsAuto = True
  1160.                    '.MaximumScaleIsAuto = True
  1161.                    .MinimumScale = 0
  1162.                     .MaximumScale = 1000
  1163.                     .MajorUnit = 250
  1164.                     .MinorUnit = 50
  1165.                 End With
  1166.             End With
  1167.         End With
  1168.     End With
  1169. End Sub
  1170.  
  1171. '-----------------------------------------
  1172. '   T E S T    F U N C T I O N S
  1173. '-----------------------------------------
  1174.  
  1175. Private Sub sTestPoint()
  1176.     Dim oTriangle1 As tTriangle
  1177.     Dim oTriangle2 As tTriangle
  1178.     Dim oP() As tXYZ
  1179.     Dim lgP As Long
  1180.  
  1181.     'Set oChrt = ActiveSheet.ChartObjects("ChrtPts")
  1182.    
  1183.     'Call sPointsCreate(4)
  1184.    
  1185.     ReDim oP(1 To 4)
  1186.     oP(1) = NewPoint(50, 1000)
  1187.     oP(2) = NewPoint(500, 3200)
  1188.     oP(3) = NewPoint(-2000, 3500)
  1189.     oP(4) = NewPoint(-2000, -3500)
  1190.  
  1191.     For lgP = LBound(oP) To UBound(oP)
  1192.         Cells(lgP, 1).Value2 = oP(lgP).x
  1193.         Cells(lgP, 2).Value2 = oP(lgP).y
  1194.     Next lgP
  1195.  
  1196.     oTriangle1 = fTriangleFromPoints(oP(1), oP(2), oP(3))
  1197.     With oTriangle1
  1198.         .V1 = 1
  1199.         .V2 = 2
  1200.         .V3 = 3
  1201.     End With
  1202.    
  1203.     Call fChrtTriangleAdd(oP(1), _
  1204.                           oP(2), _
  1205.                           oP(3), _
  1206.                           ActiveSheet.ChartObjects("ChrtPts"))
  1207.    
  1208.     oTriangle2 = fTriangleFromPoints(oP(3), oP(2), oP(1))
  1209.     With oTriangle2
  1210.         .V1 = 3
  1211.         .V2 = 2
  1212.         .V3 = 4
  1213.     End With
  1214.     Call fChrtTriangleAdd(oP(3), _
  1215.                           oP(2), _
  1216.                           oP(4), _
  1217.                           ActiveSheet.ChartObjects("ChrtPts"))
  1218.    
  1219.     'Call fChrtCircleAdd(oTriangle.Center, _
  1220.                         Sqr(oTriangle.R²), _
  1221.                         ActiveSheet.ChartObjects("ChrtPts"))
  1222. Stop
  1223. End Sub
  1224.  
  1225. '----------------------------------------------
  1226. '    T E S T
  1227. '----------------------------------------------
  1228.  
  1229. Private Function xTriangle(ByRef oP1 As tXYZ, ByRef oP2 As tXYZ, ByRef oP3 As tXYZ) As tXYZ
  1230.     With xTriangle
  1231.         Dim B As tXYZ
  1232.         Dim C As tXYZ
  1233.         Dim D As Double
  1234.         Dim Bx² As Double
  1235.         Dim By² As Double
  1236.         Dim Cx² As Double
  1237.         Dim Cy² As Double
  1238.        
  1239.         B.x = oP2.x - oP1.x
  1240.         B.y = oP2.y - oP1.y
  1241.         C.x = oP3.x - oP1.x
  1242.         C.y = oP3.y - oP1.y
  1243.         Bx² = B.x * B.x: By² = B.y * B.y
  1244.         Cx² = C.x * C.x: Cy² = C.y * C.y
  1245.         D = (2 * (B.x * C.y - B.y * C.x))
  1246.         If D <> 0 Then
  1247.             D = 1 / D
  1248.             .x = D * (C.y * (Bx² + By²) - B.y * (Cx² + Cy²))
  1249.             .y = D * (B.x * (Cx² + Cy²) - C.x * (Bx² + By²))
  1250.             '.R² = (.Center.X * .Center.X) + (.Center.Y * .Center.Y)
  1251.            .x = .x + oP1.x
  1252.             .y = .y + oP1.y
  1253.         End If
  1254.     End With
  1255. End Function
  1256.  
  1257. Sub sTestingPerformance()
  1258.     Dim lgT As Long
  1259.     Dim dtDate As Date
  1260.     Dim oP1 As tXYZ, oP2 As tXYZ, oP3 As tXYZ
  1261.     Dim T As Single
  1262.     T = Timer
  1263.  
  1264.     oP1 = NewPoint(50, 1000)
  1265.     oP2 = NewPoint(500, 3200)
  1266.     oP3 = NewPoint(-2000, 3500)
  1267.    
  1268.     'dtDate = VBA.Now()
  1269.    For lgT = 1 To 1000000000
  1270.         Call xTriangle(oP1, oP2, oP3)
  1271.     Next lgT
  1272.     Debug.Print "Completed processing " & lgT - 1 & " iterations in " & Round(Timer - T, 2) & " seconds."
  1273.     'Debug.Print dtDate & "  --- " & VBA.Now()
  1274.    Debug.Print
  1275.     Beep
  1276.    
  1277. End Sub
  1278.  
Tags: vba
Advertisement
Comments
  • Ben_S
    1 year
    # VBScript 0.60 KB | 0 0
    1. I was searching for a VBA algorithm for Delaunay when I came across this.  Unfortunately, the VBA code here is poorly formatted and takes some guesswork to straighten out.  I have done so and posted it in the link below.  I'm not sure how to implement the functions so you as the reader will have to figure that out unless the author makes some updates.
    2.  
    3. My code compiles in Excel VBA, but it is incomplete because of two missing function and the value for EPSILON.  The missing functions are fQuickSortArrayDbl() and fDistance2DPoints().  I created placeholder function stubs just so the code could compile.
Add Comment
Please, Sign In to add comment
Advertisement