VBA code for Delaunay and ConvexHull

Option Explicit

' VBA code for Delaunay and ConvexHull
' Created by auadmin on http://www.audeser.com/vba-delaunay-convexhull/
'
' The VBA formatting was scrambled so this is my attempt to sort it out.
' Ben Sacherich - 5/17/2023
' https://pastebin.com/BJvEUnD8
'
' Constant EPSILON was not defined.
' Function fQuickSortArrayDbl() was missing.
' Function fDistance2DPoints() was missing.
'

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     ' This is not used and won't compile so Ben commented it out.

Private Const PI_8th As Double = 0.392699081698724
Private Const PI As Double = 3.14159265358979

'# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
Const EPSILON = 1       ' ### Placeholder for missing code ###

Function fDistance2DPoints(x As tXYZ, y As tXYZ)
    ' ### Placeholder for missing function ###
    Stop
End Function

Function fQuickSortArrayDbl(arry As Variant, i, j, k, bAscendent As Boolean) As Variant
    ' ### Placeholder for missing function ###
    Stop
End Function
'# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Private Function fDrawTriangles() As Boolean
' # This function is not used #
    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
' # This function is not used #

    '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
' # This function was missing but is not needed. #
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 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
    Dim T As Single
    T = Timer

    oP1 = NewPoint(50, 1000)
    oP2 = NewPoint(500, 3200)
    oP3 = NewPoint(-2000, 3500)

    'dtDate = VBA.Now()
    For lgT = 1 To 1000000000
        Call xTriangle(oP1, oP2, oP3)
    Next lgT
    Debug.Print "Completed processing " & lgT - 1 & " iterations in " & Round(Timer - T, 2) & " seconds."
    'Debug.Print dtDate & "  --- " & VBA.Now()
    Debug.Print
    Beep

End Sub