Autodesk 2018 Direct Download Links

Autodesk 2018 Direct Download Links (Until Available on Virtual Agent)

Issue:  You are having a download failure error, similar to previous releases, that are causing installation errors of your 2018 product.  However, the Autodesk Virtual Agent does not yet list the products available for 2018.

Work-Around:  Until the Virtual Agent is updated, we collective peers will try to provide Official direct download links from Autodesk.  This list may be updated as needed and links may or may not work as products are released.  You can also try the Browser Download method from your Autodesk Accounts Management page.

If you have an official direct download link, that is tested and working, or have a URL of a direct download link that should work once the products are released, then please feel free to update the links in the comments below.  Please only use official Autodesk links.  Piracy is not tolerated here.

AUTODESK 2018 DIRECT DOWNLOAD LINKS

AutoCAD 2018

English 32 bit

English 64 bit – Part 1
English 64 bit – Part 2

AutoCAD LT 2018English 32 bit

English 64 bit

Inventor Professional 2018 x64Part 1
Part 2
Part 3

Revit Live 2018
Download
Navisworks Manage 2018Part 1
Part 2
Navisworks Simulate 2018Part 1
Part 2

Here are more working links for other 2018 products:

DWG TrueView 2018

32 bit
64 bit

Inventor View 2018Download

Recap 360 ProDownload

AutoCAD Raster Design 2018

32-bit

64-bit

AutoCAD MAP 3D 2018 (x64)

Part 1

Part 2

AutoCAD Electrical 2018

32-bit Part 1

32-bit Part 2

64-bit Part 1

64-bit Part 2

AutoCAD Mechanical 2018

32-bit

64-bit Part 1

64-bit Part 2

AutoCAD Architecture 2018

32-bit Part 1

32-bit Part 2

64-bit Part 1

64-bit Part 2

AutoCAD MEP 2018

32-bit Part 1

32-bit Part 2

32-bit Part 3

64-bit Part 1

64-bit Part 2

64-bit Part 3

Autodesk Sketchbook Pro Enterprise 2018 (x64)

Download

Inventor LT 2018 (x64)

Part 1

Part 2

Vault 2018 Pro

Vault Pro 2018 Server

Vault Pro 2018 Client

Vault 2018 File Server

Vault 2018 File Server

Vault Basic 2018

Vault Basic 2018 Server

Vault Basic 2018 Client

Here are additional links for Vault:

Autodesk Vault 2018 Basic – Client (x64)

Download

Autodesk Vault 2018 Basic – Server (x64)Download

I wonder if the e-fulfillment links expire after a certain amount of time?  Thanks for the updated links Darren.

Advance Steel 2018 (64 bit)

Part 1

Part 2

AutoDesk Advance Steel 2018 (x64)

Part 1

Part 2

AutoDesk Alias Design 2018 (x64)Part 1
Part 2

AutoDesk Alias Surface 2018 (x64)Download

AutoDesk Alias Speedform 2018 (x64)Download

Autodesk Moldflow Insight Ultimate (x64)

Download

Autodesk Moldflow Adviser Ultimate (x64)Part 1
Part 2

AutoCAD Plant 3D English 2018 (x64)Part 1

Part 2

Revit 2018

Part 1

Part 2

Part 3

Civil 3D 2018

Part 1

Part 2

Part 3

Vehicle Tracking 2018

Vehicle Tracking 2018

Revit LT 2018

Part 1

Part 2

Revit Server 2018 (x64)Download

I also wanted to provide an additional link for the Autodesk Network License Manager for 2018

Autodesk Network License Manager 2017/2018

Info Here

Infraworks 2018

Infraworks 2018
Autodesk Building Design Suite Premium 2018 (x64)*Part 1
Part 2
Part 3
Part 4
Part 5

* This Suite is still available if you have maintained your maintenance Subscription for the Building Design Suite

AEC Collection 2018

AutoCAD 2018English 32 bit

English 64 bit – Part 1
English 64 bit – Part 2

Navisworks Manage 2018Part 1
Part 2

AutoCAD MAP 3D 2018 (x64)

Part 1

Part 2

AutoCAD Architecture 2018

32-bit Part 1

32-bit Part 2

64-bit Part 1

64-bit Part 2

AutoCAD MEP 2018

32-bit Part 1

32-bit Part 2

32-bit Part 3

64-bit Part 1

64-bit Part 2

64-bit Part 3

AutoCAD Plant 3D English 2018 (x64)Part 1

Part 2

3ds Max 2018

Part 1

Part 2

Revit 2018

Part 1

Part 2

Part 3

Civil 3D 2018

Part 1

Part 2

Part 3

Vehicle Tracking 2018

Vehicle Tracking 2018

Revit Server 2018 (x64)Download

Recap 360 ProDownload

Infraworks 2018

Infraworks 2018

AutoCAD Raster Design 2018

32-bit

64-bit

AutoCAD Electrical 2018*

32-bit Part 1

32-bit Part 2

64-bit Part 1

64-bit Part 2

*Autodesks website says Autocad Electrical comes with the AEC Collection i don’t think thats correct but adding anyway

Product Design Ultimate 2018

Part 1

Part 2

Part 3

Part 4

Part 5

Here are some new links as well as some corrected links.  Thanks for your patience.

3DS Max 2018 (x64)Part 1
Part 2

AutoCAD Civil 3D English 2018 (x64)Part 1
Part 2
Part 3

Revit 2018 (x64)Part 1
Part 2
Part 3

Revit LT 2018 (x64)Part 1
Part 2

Vehicle tracking 2018Download

Autodesk Building Design Suite Ultimate 2018 (x64)Part 1
Part 2
Part 3
Part 4
Part 5
Part 6

Autodesk Product Design Suite Ultimate 2018 (x64)Part 1
Part 2
Part 3
Part 4
Part 5

Autodesk Infrastructure Design Suite Premium 2018 (x64)Part 1
Part 2
Part 3
Part 4
Part 5

Autodesk InfraWorks 360 Pro 2018 (x64)Download

Autodesk Robot Structural Analysis Professional 2018 (x64)Download

Autodesk Factory Design Suite Ultimate 2018 (x64)Part 1
Part 2
Part 3
Part 4
Part 5
Part 6

Autodesk Infrastructure Design Suite Ultimate 2018 (x64)Part 1
Part 2
Part 3
Part 4
Part 5

AEC Collection 2018 (x64)AutoCAD 2018Part 1
Part 2
Revit 2018Part 1
Part 2
Part 3
Revit Server Download
Civil 3D 2018Part 1
Part 2
Part 3
Infraworks 2018Download
Navisworks Manage 2018Part 1
Part 2
AutoCAD Raster Design 2018Download
3DS Max 2018Part 1
Part 2
Vehicle Tracking 2018Download
AutoCAD MAP 3D 2018Part 1
Part 2
AutoCAD Architecture 2018Part 1
Part 2
AutoCAD Electrical 2018Part 1
Part 2
AutoCAD MEP 2018Part 1
Part 2
Part 3
AutoCAD Plant 3D 2018Part 1
Part 2
Recap 360 ProDownload


AutoCAD Mobile AppDownload
FormIt Pro AppDownload
Insight Plug-in for Revit 2018Download
Structural Analysis for Revit 2018Subscriber Login

Autodesk Nastran In-CAD 2018 (x64)Download

Autodesk Nastran 2018 (x64)Download

Autodesk HSM Ultimate 2018 for Inventor and Solidworks (x64)Part 1
Part 2

Autodesk Simulation CFD 2018 (x64)Part 1
Part 2

A couple of Revit 2018 add-ins were made available last night

Revit 2018 Steel Connections

Revit 2018 Site Designer

Autodesk InfraWorks 360 Pro 2018 (x64)

Link

Vault Workgroup 2018:

Vault Workgroup 2018 (Server)

Vault Workgroup 2018 (Client)

Polyline simplification

I think I have posted something about Douglas-Peucker algorithm, to reduce the number of vertices of a polyline.

Polyline simplification is the process of reducing the resolution of a polyline, achieved by removing vertices and edges, while maintaining a good approximation of the original curve. In the end is a compromise between waste of resources and level of detail-the resolution of the polyline-.  There are some algorithms you can recall to:

Simplification algorithms

  • Nth point – A naive algorithm that keeps only each nth point
  • Distance between points – Removes successive points that are clustered together
  • Perpendicular distance – Removes points based on their distance to the line segment defined by their left and right neighbors
  • Reumann-Witkam – Shifts a strip along the polyline and removes points that fall outside
  • Opheim – Similar to Reumann-Witkam, but constrains the search area using a minimum and maximum tolerance
  • Lang – Similar to the Perpendicular distance routine, but instead of looking only at direct neighbors, an entire search region is processed
  • Douglas-Peucker – A classic simplification algorithm that provides an excellent approximation of the original line

Error algorithms

  • Positional errors – Distance of each point from an original polyline to its simplification

In this articles there is a bit more information:

 

Genetic algorithm with VBA

For more than a couple of years I’ve been dealing with a nightmare-like problem that I could not get solved with Excel.

The solution I was working with implies choosing between nearly infinite combinatorial possibilities. The problem in question has a lot of dependencies in a variable number of variables you can set, is an open problem in its definition, and I was looking for the optimal one. You bet I will never achieve it. I was also stuck wandering if I could face situations where several options were possible; and I couldn’t discard that as not pausible.

I considered GA as an alternative, but never put the enough time on studying the possibilities, considering I have not -or do not know, that is the same- a tool for GA in VBA. Times I spend  a bit more learning GA in YouTube or in blogs did not scratch the surface, and the only XLS file I played with was very basic, and came from this old post from Paras Chopra blog, and another one from Dermont Balson’s InsaneExcel death blog (under Artificial Intelligence).

Continue reading “Genetic algorithm with VBA”

Function fitter VBA – Conjugate Gradient

I found a nice post about applying conjugate gradient for a function fitter. They show how to use the conjugate gradient with POLAK and RIBIÈRE factor to fit a circle to a set of points. The PDF is very explicative, so I learnt a lot about the mathematics behaind this. Now it will be a lot easier to apply this to a clothoid.

The original code is a Java one, so here is translated to VBA. Will come very handy for an “alignment from points” procedure, although you will need to separate the different alignments before applying this.

Option Explicit

' Class fitting a circle to a set of points.
' This class implements the fitting algorithms described in the paper 
' Finding the circle that best fits a set of points
' @author Luc Maisonobe

' Copyright (c) 2005-2007, Luc Maisonobe
' All rights reserved.
'
' Redistribution and use in source and binary forms, with
' or without modification, are permitted provided that
' the following conditions are met:
'
'    Redistributions of source code must retain the
'    above copyright notice, this list of conditions and
'    the following disclaimer.
'    Redistributions in binary form must reproduce the
'    above copyright notice, this list of conditions and
'    the following disclaimer in the documentation
'    and/or other materials provided with the
'    distribution.
'    Neither the names of spaceroots.org, spaceroots.com
'    nor the names of their contributors may be used to
'    endorse or promote products derived from this
'    software without specific prior written permission.
'
' THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
' CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED
' WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
' WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
' PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
' THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY
' DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
' CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
' PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
' USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
' HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
' IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
' NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
' USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
' POSSIBILITY OF SUCH DAMAGE.

Public Type tXYZ
    X As Double
    Y As Double
    Z As Double
End Type

Private Type tCircleFitter
    Pt() As tXYZ
    Center As tXYZ
    Radius As Double
End Type

Private fitter As tCircleFitter
Private r´ As Double       ' optimum circle radius
Private Jacb As Double
Private dJx As Double
Private dJy As Double

Private Const PI As Double = 3.14159265358979
Private Const EPSILON As Double = 0.00000001

Private Sub sCircleFitter_Main()
    On Local Error GoTo CtrlErr

    Dim iter As Long

    ' fit a circle to the test points
    Call fPointsInitialize(fitter.Pt(), 15, 100, 100, 35)

    'fitter = New CircleFitter
    Call fitterInitialize(fitter.Pt())
'Debug.Print "initial circle: " & fitter.Center.X & " " & fitter.Center.Y & " " & fitter.Radius

    ' minimize the residuals
    iter = minimize(100, 0.1, EPSILON)
'Debug.Print "converged after " & iter & " iterations"
'Debug.Print "final circle: " & fitter.Center.X & " " & fitter.Center.Y & " " & fitter.Radius

ExitProc:
    Exit Sub

CtrlErr:
    Stop
End Sub

Public 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 fPointsInitialize(ByRef oPt() As tXYZ, _
                                   ByVal lgPts As Long, _
                                   ByVal Xcenter As Double, _
                                   ByVal Ycenter As Double, _
                                   ByVal Radius As Double) As Boolean
' Initialize an approximate circle based on all triplets.
' @param lgPts number of sample points
    Dim lgPt As Long
    Dim dbAngleRAD As Double

    ReDim oPt(0 To lgPts - 1)
    VBA.Randomize
    For lgPt = 1 To lgPts
        dbAngleRAD = (Rnd() * 2 * PI)
        With oPt(lgPt - 1)
            .X = (Xcenter + Radius * Cos(dbAngleRAD)) + (CLng(IIf(Rnd() > 0.5, 1, -1)) * Rnd() * 0.01)
            .Y = (Ycenter + Radius * Sin(dbAngleRAD)) + (CLng(IIf(Rnd() > 0.5, 1, -1)) * Rnd() * 0.01)
        End With
    Next lgPt
End Function

Public Sub fitterInitialize(ByRef oPt() As tXYZ)
' Initialize an approximate circle based on all triplets.
' @param oPt() circular ring sample points
' @exception Error if all points are aligned
 
    Dim n As Long 'number of filtered points
    Dim lgPt1 As Long
    Dim lgPt2 As Long
    Dim lgPt3 As Long
    Dim cc As tXYZ
    Dim lgRetVal As Long

    With fitter
        With .Center
            .X = 0
            .Y = 0
        End With
        .Radius = 0
    End With

    ' analyze all possible points triplets
    n = 0
    For lgPt1 = LBound(oPt) To UBound(oPt) - 2
        For lgPt2 = lgPt1 + 1 To UBound(oPt) - 1
            For lgPt3 = lgPt2 + 1 To UBound(oPt)
                ' compute the triangle circumcenter
                ' Check points are not aligned
                If Not fAligned(oPt(lgPt1), oPt(lgPt2), oPt(lgPt3)) Then
                    cc = circumcenter(oPt(lgPt1), oPt(lgPt2), oPt(lgPt3))
                    ' the points are not aligned, we have a circumcenter
                    n = (n + 1)
                    fitter.Center.X = (fitter.Center.X + cc.X)
                    fitter.Center.Y = (fitter.Center.Y + cc.Y)
                End If
            Next lgPt3
        Next lgPt2
    Next lgPt1
    
    If (n = 0) Then
        lgRetVal = VBA.MsgBox("all points are aligned")
    Else
        ' initialize using the circumcenters average
        With fitter
            With .Center
                .X = .X / n
                .Y = .Y / n
            End With
            Call updateRadius(.Pt())
        End With
    End If
End Sub

Public Function fAligned(ByRef oPt1 As tXYZ, _
                         ByRef oPt2 As tXYZ, _
                         ByRef oPt3 As tXYZ) As Boolean
    fAligned = (((oPt2.X - oPt1.X) * (oPt3.Y - oPt1.Y)) - ((oPt2.Y - oPt1.Y) * (oPt3.X - oPt1.X)) < EPSILON)
End Function

Private Function updateRadius(ByRef oPt() As tXYZ) As Boolean
' Update the circle radius

    Dim lgPt As Long
    Dim dx As Double
    Dim dy As Double
    
    r´ = 0
    If Not (Not oPt) Then
        For lgPt = LBound(oPt) To UBound(oPt)
            dx = (oPt(lgPt).X - fitter.Center.X)
            dy = (oPt(lgPt).Y - fitter.Center.Y)
            r´ = (r´ + VBA.Sqr(((dx * dx) + (dy * dy))))
        Next lgPt
        
        r´ = r´ / (UBound(oPt) - LBound(oPt) + 1)
        fitter.Radius = r´
    End If
End Function

Private Function circumcenter(ByRef PI As tXYZ, ByRef pJ As tXYZ, ByRef pK As tXYZ) As tXYZ
' Compute the circumcenter of three points.
' @param pI first point
' @param pJ second point
' @param pK third point
' @return circumcenter of pI, pJ and pK or null if the points are aligned

    Dim dIJ As tXYZ
    Dim dJK As tXYZ
    Dim dKI As tXYZ
    Dim sqI As Double
    Dim sqJ As Double
    Dim sqK As Double
    Dim det As Double

    dIJ = NewPoint((pJ.X - PI.X), (pJ.Y - PI.Y))
    dJK = NewPoint((pK.X - pJ.X), (pK.Y - pJ.Y))
    dKI = NewPoint((PI.X - pK.X), (PI.Y - pK.Y))
    sqI = ((PI.X * PI.X) + (PI.Y * PI.Y))
    sqJ = ((pJ.X * pJ.X) + (pJ.Y * pJ.Y))
    sqK = ((pK.X * pK.X) + (pK.Y * pK.Y))

    ' determinant of the linear system: 0 for aligned points
    det = ((dJK.X * dIJ.Y) - (dIJ.X * dJK.Y))

    If (VBA.Abs(det) < EPSILON) Then Exit Function ' points are almost aligned, we cannot compute the circumcenter
    
    ' beware, there is a minus sign on Y coordinate!
    circumcenter = NewPoint(((sqI * dJK.Y) + ((sqJ * dKI.Y) + (sqK * dIJ.Y))) / (2 * det), _
                           -((sqI * dJK.X) + ((sqJ * dKI.X) + (sqK * dIJ.X))) / (2 * det))
End Function

Public Function minimize(ByVal iterMax As Integer, ByVal innerThreshold As Double, ByVal outerThreshold As Double) As Long
' Minimize the distance residuals between the points and the circle.
' We use a non-linear conjugate gradient method with the Polak and
' Ribiere coefficient for the computation of the search direction. The
' inner minimization along the search direction is performed using a
' few Newton steps. It is worthless to spend too much time on this inner
' minimization, so the convergence threshold can be rather large.
'
' @param maxIter maximal iterations number on the inner loop (cumulated across outer loop iterations)
' @param innerThreshold inner loop threshold, as a relative difference on the cost function value between the two last iterations
' @param outerThreshold outer loop threshold, as a relative difference on the cost function value between the two last iterations
' @return number of inner loop iterations performed (cumulated across outer loop iterations)
' @exception LocalException if we come accross a singularity or if we exceed the maximal number of iterations
    
    Dim lgRetVal As Long
    Dim u As Double
    Dim v As Double
    Dim beta As Double
    Dim lambda As Double
    Dim innerJ As Double
    Dim previousJ As Double: previousJ = Jacb
    Dim previousV As Double: previousV = 0
    Dim previousU As Double: previousU = 0
    Dim previousdJy As Double: previousdJy = 0
    Dim previousdJx As Double: previousdJx = 0
    Dim iterations As Integer: iterations = 0
  
    Call computeCost(fitter.Pt())

    If ((Jacb < EPSILON) Or (VBA.Sqr(((dJx * dJx) + (dJy * dJy))) < EPSILON)) Then
        ' we consider we are already at a local minimum
        minimize = 0
        Exit Function
    End If
    
    Do While (iterations < iterMax)
        ' search direction
        u = (dJx * -1)
        v = (dJy * -1)
        If (iterations <> 0) Then
            ' Polak-Ribiere coefficient
            beta = (((dJx * (dJx - previousdJx)) + (dJy * (dJy - previousdJy))) / ((previousdJx * previousdJx) + (previousdJy * previousdJy)))
            u = (u + (beta * previousU))
            v = (v + (beta * previousV))
        End If
        
        previousdJx = dJx
        previousdJy = dJy
        previousU = u
        previousV = v

        ' rough minimization along the search direction
        Do
            With fitter
                innerJ = Jacb
                lambda = newtonStep(.Pt(), u, v)
                .Center.X = (fitter.Center.X + (lambda * u))
                .Center.Y = (fitter.Center.Y + (lambda * v))
                Call updateRadius(.Pt())
                Call computeCost(.Pt())
            
                iterations = iterations + 1
            End With
        Loop While ((iterations < iterMax) And ((VBA.Abs((Jacb - innerJ)) / Jacb) > innerThreshold))
        
        If ((VBA.Abs((Jacb - previousJ)) / Jacb) < outerThreshold) Then
            minimize = iterations
            Exit Function
        End If
        
        previousJ = Jacb
    Loop
    
    lgRetVal = VBA.MsgBox("unable to converge after " & iterMax & " iterations")
End Function

Private Sub computeCost(ByRef oPt() As tXYZ)
' Compute the cost function and its gradient.
    Dim lgRetVal As Long
    Dim lgPt As Long
    Dim dx As Double
    Dim dy As Double
    Dim di As Double
    Dim dr As Double
    Dim ratio As Double

    Jacb = 0
    dJx = 0
    dJy = 0
    For lgPt = LBound(oPt) To UBound(oPt)
        dx = (oPt(lgPt).X - fitter.Center.X)
        dy = (oPt(lgPt).Y - fitter.Center.Y)
        di = VBA.Sqr(((dx * dx) + (dy * dy)))
        If (di < EPSILON) Then lgRetVal = VBA.MsgBox(("cost singularity: point at the circle center"))
        
        dr = (di - r´)
        ratio = (dr / di)
        Jacb = (Jacb + (dr * (di + r´)))
        dJx = (dJx + (dx * ratio))
        dJy = (dJy + (dy * ratio))
    Next lgPt
    
    dJx = (dJx * 2)
    dJy = (dJy * 2)
End Sub

Private Function newtonStep(ByRef oPt() As tXYZ, ByVal u As Double, ByVal v As Double) As Double
' Compute the length of the Newton step in the search direction.
' @param u abscissa of the search direction
' @param v ordinate of the search direction
' @return value of the step along the search direction

    Dim sumFac2R As Double
    Dim sum1 As Double
    Dim sum2 As Double
    Dim sumFac As Double
    Dim lgPt As Long
    Dim dx As Double
    Dim dy As Double
    Dim di As Double
    Dim coeff1 As Double
    Dim coeff2 As Double

    sumFac2R = 0
    sum1 = 0
    sum2 = 0
    sumFac = 0

    ' compute the first and second derivatives of the cost along the specified search direction
    For lgPt = LBound(oPt) To UBound(oPt)
        dx = (fitter.Center.X - oPt(lgPt).X)
        dy = (fitter.Center.Y - oPt(lgPt).Y)
        di = VBA.Sqr(((dx * dx) + (dy * dy)))
        coeff1 = (((dx * u) + (dy * v)) / di)
        coeff2 = (di - r´)
        sum1 = (sum1 + (coeff1 * coeff2))
        sum2 = (sum2 + (coeff2 / di))
        sumFac = (sumFac + coeff1)
        sumFac2R = (sumFac2R + (coeff1 * (coeff1 / di)))
    Next lgPt
    
    ' step length attempting to nullify the first derivative
    newtonStep = ((sum1 / (((((u * u) + (v * v)) * sum2) - (sumFac * (sumFac / (UBound(oPt) - LBound(oPt) + 1)))) + (r´ * sumFac2R))) * -1)
End Function

'Public Function getCenter(ByRef fitter As tCircleFitter) As Boolean
'    fitter.Center = oCenter: getCenter = True
'End Function

'Public Function getRadius(ByRef fitter As tCircleFitter) As Boolean
'    fitter.Radius = r´: getRadius = True
'End Function

VBA Concave Hull

In previous post was shown an algorithm to obtain the convex hull of a set of points.

But the convex hull, beeing extremely fast, has some disadvantages, finding the most important that it acts like a rubber bounding a figurine so, although  it can embrace all the set of points, it will left big spare spaces from that set to the hull.

A more interesting bounding polygon will be a concave hull. Found this paper about an algorithm to do so, and two implementations, one in JavaScript and the other in C++. Trying to port to VBA, here I left them, so they are far from finished and tested, but for sure will return to it since it should be the algorithm on preference to the convex hull.

There is more information on this topic on this paper on this other one,  also on this one, and finally this other.

First the JS version, that I think is more interesting:

'https://github.com/mapbox/concaveman/blob/master/index.js

Option Explicit

Private Function concaveman(ByRef points() As tXYZ, ByVal concavity As Boolean, ByVal lengthThreshold As Double)
    ' a relative measure of concavity; higher value means simpler hull
    concavity = fMax(0, VBA.IIf(concavity = undefined, 2, concavity))
    
    ' when a segment goes below this length threshold, it won't be drilled down further
    lengthThreshold = (lengthThreshold Or 0)
    
    ' start with a convex hull of the points
    Dim hull: hull = fastConvexHull(points)
    
    ' index the points with an R-tree
    Dim tree:tree  = rbush(16, [, Microsoft.VisualBasic.ChrW(91), Microsoft.VisualBasic.ChrW(91), Microsoft.VisualBasic.ChrW(91), Microsoft.VisualBasic.ChrW(91))
    Load (points)
    
    ' turn the convex hull into a linked list and populate the initial edge queue with the nodes
    Dim queue
    Dim last
    
    Do While (i < hull.Length)
        Dim p: p = hull(i)
        Dim i: i = 0
        tree.Remove (p)
        last = insertNode(p, last)
        queue.push (last)
        i = (i + 1)
    Loop
    
    ' index the segments with an R-tree (for intersection checks)
    Dim segTree: segTree = rbush(16)
    
    i = 0
    Do While (i < queue.Length)
        segTree.Insert (updateBBox(queue(i)))
        i = (i + 1)
    Loop
    
    Dim sqConcavity: sqConcavity = (concavity * concavity)
    Dim sqLenThreshold: sqLenThreshold = (lengthThreshold * lengthThreshold)
    
    ' process edges one by one
    While queue.Length
        Dim node: node = queue.Shift
        Dim a: a = node.p
        Dim b: b = node.Next.p
        
        ' skip the edge if it's already short enough
        Dim sqLen: sqLen = getSqDist(a, b)
        If (sqLen < sqLenThreshold) Then
            'TODO: Warning!!! continue If
        End If
        
        Dim maxSqLen: maxSqLen = (sqLen / sqConcavity)
        ' find the best connection point for the current edge to flex inward to
        p = findCandidate(tree, node.prev.p, a, b, node.Next.Next.p, maxSqLen, segTree)
        ' if we found a connection and it satisfies our concavity measure
        If (p And (Math.Min(getSqDist(p, a), getSqDist(p, b)) <= maxSqLen)) Then
            ' connect the edge endpoints through this point and add 2 new edges to the queue
            queue.push (node)
            queue.push (insertNode(p, node))
            ' update point and segment indexes
            tree.Remove (p)
            segTree.Remove (node)
            segTree.Insert (updateBBox(node))
            segTree.Insert (updateBBox(node.Next))
        End If
    Loop
    
    ' convert the resulting hull linked list to an array of points
    node = last
    Dim concave
    
    Do
        concave.push (node.p)
        node = node.Next
    Loop While node <> last
    
    concave.push (node.p)
    concaveman() = concave()
End Function

Private Function findCandidate(ByVal tree, ByVal a, ByVal b, ByVal c, ByVal d, ByVal maxDist, ByVal segTree) As tXYZ
    Dim queue: queue = NewQueue(Nothing, compareDist)
    Dim node: node = tree.Data
    
    ' search through the point R-tree with a depth-first search using a priority queue
    ' in the order of distance to the edge (b, c)
    While node
        Dim i: i = 0
        Do While (i < node.Children.Length) Dim child = node.children(i) Dim dist = node.leaf 'TODO: Warning!!!, inline IF is not supported ? If (dist > maxDist) Then
                'TODO: Warning!!! continue If
            End If
            
            ' skip the node if it's farther than we ever need
            queue.push({, node:= child, dist:= dist)
            i = (i + 1)
        Loop
    End While
    
    
    While (queue.Length And Not queue.peek.node.Children)
        Dim item: item = queue.pop
        Dim p: p = item.node
        ' skip all points that are as close to adjacent edges (a,b) and (c,d),
        ' and points that would introduce self-intersections when connected
        Dim d0: d0 = sqSegDist(p, a, b)
        Dim d1: d1 = sqSegDist(p, c, d)
        If ((item.dist < d0) And ((item.dist < d1) And (noIntersections(b, p, segTree) And noIntersections(c, p, segTree)))) Then Return p End If Loop node = queue.pop If node Then node = node.node End Function Private Function compareDist(ByVal a, ByVal b) As Boolean compareDist = (a.dist - b.dist) End Function ' square distance from a segment bounding box to the given one Private Function sqSegBoxDist(ByVal a, ByVal b, ByVal bbox) As Double If (inside(a, bbox) Or inside(b, bbox)) Then sqSegBoxDist = 0: Exit Function Dim d1 = sqSegSegDist(a(0), a(1), b(0), b(1), bbox.minX, bbox.minY, bbox.maxX, bbox.minY) If d1 = 0 Then sqSegBoxDist = 0: Exit Function Dim d2 = sqSegSegDist(a(0), a(1), b(0), b(1), bbox.minX, bbox.minY, bbox.minX, bbox.maxY) If d2 = 0 Then sqSegBoxDist = 0: Exit Function Dim d3 = sqSegSegDist(a(0), a(1), b(0), b(1), bbox.maxX, bbox.minY, bbox.maxX, bbox.maxY) If d3 = 0 Then sqSegBoxDist = 0: Exit Function Dim d4 = sqSegSegDist(a(0), a(1), b(0), b(1), bbox.minX, bbox.maxY, bbox.maxX, bbox.maxY) If d4 = 0 Then sqSegBoxDist = 0: Exit Function sqSegBoxDist = fMin(d1, d2, d3, d4) End Function Private Function inside(ByVal a, ByVal bbox) As Boolean inside = ((a(0) >= bbox.minX) And _
             ((a(0) <= bbox.maxX) And _ ((a(1) >= bbox.minY) And _
             (a(1) <= bbox.maxY))))
End Function

Private Function noIntersections(ByVal a, ByVal b, ByVal segTree) As Boolean
' check if the edge (a,b) doesn't intersect any other edges
    Dim minX: minX = fMin(a(0), b(0))
    Dim minY: minY = fMin(a(1), b(1))
    Dim maxX: maxX = fMax(a(0), b(0))
    Dim maxY: maxY = fMax(a(1), b(1))
    Dim edges
    
    edges = segTree.Search(minX:=minX, minY:=minY, maxX:=maxX, maxY:=maxY)
    For i = 0 To edges.Length
        If (intersects(edges(i).p, edges(i).Next.p, a, b)) Then
            noIntersections = False: Exit Function
        End If
    Next i
    
    noIntersections = True
End Function

Private Function intersects(ByRef p1, ByRef q1, ByRef p2, ByRef q2) As Boolean
' check if the edges (p1,q1) and (p2,q2) intersect
    intersects = (p1 <> q2) And _
                 (q1 <> p2) And _
                 (orient(p1, q1, p2) > 0) <> (orient(p1, q1, q2) > 0) And _
                 (orient(p2, q2, p1) > 0) <> (orient(p2, q2, q1) > 0)
End Function

Private Function updateBBox(ByVal p As tXYZ) As Boolean
' update the bounding box of a node's edge
    Dim p1 = node.p
    Dim p2 = node.next.p
    node.minX = Math.Min(p1(0), p2(0))
    node.minY = Math.Min(p1(1), p2(1))
    node.maxX = Math.Max(p1(0), p2(0))
    node.maxY = Math.Max(p1(1), p2(1))
    Return node
End Function

Private Function fastConvexHull(ByVal Unknown As points) As tXYZ()
' speed up convex hull by filtering out points inside quadrilateral formed by 4 extreme points
    Dim left = points(0)
    Dim top = points(0)
    Dim right = points(0)
    Dim bottom = points(0)
    
    ' find the leftmost, rightmost, topmost and bottommost points
    Dim i = 0
    Do While (i < points.Length)
        Dim p = points(i)
        If (p(0) < Left(0)) Then Left = p End If If (p(0) > Right(0)) Then
            Right = p
        End If
        
        If (p(1) < Top(1)) Then Top = p End If If (p(1) > Bottom(1)) Then
            Bottom = p
        End If
        
        i = (i + 1)
    Loop
    
    ' filter out points that are inside the resulting quadrilateral
    Dim cull
    Left
    Top
    Right
    Bottom
    Dim filtered = cull.slice
    i = 0
    Do While (i < points.Length)
        If Not pointInPolygon(points(i), cull) Then
            filtered.push (points(i))
        End If
        
        i = (i + 1)
    Loop
    
    ' get convex hull around the filtered points
    Dim indices = convexHull(filtered)
    ' return the hull as array of points (rather than indices)
    Dim hull
    i = 0
    Do While (i < indices.Length) hull.push (filtered(indices(i))) i = (i + 1) Loop Return hull End Function ' create a new node in a doubly linked list Private Function insertNode(ByRef p, ByRef prev) As tXYZ Dim node p: p prev: Nothing Next: Nothing minX: 0 minY: 0 maxX: 0 maxY: 0 If Not prev Then node.prev = node node.Next = node Else node.Next = prev.Next node.prev = prev prev.Next.prev = node prev.Next = node End If Return node End Function Private Function getSqDist(ByRef p1 As tXYZ, ByRef p2 As tXYZ) As Double ' square distance between 2 points Dim dy As Double: dy = (p1(1) - p2(1)) Dim dx As Double: dx = (p1(0) - p2(0)) getSqDist = ((dx * dx) + (dy * dy)) End Function Private Function sqSegDist(ByRef p As tXYZ, ByRef p1 As tXYZ, ByRef p2 As tXYZ) As Double ' square distance from a point to a segment Dim dy As Double: dy = (p2(1) - y) Dim x As Double: x = p1(0) Dim y As Double: y = p1(1) Dim dx As Double: dx = (p2(0) - x) Dim t As Double t = ((((p(0) - x) * dx) + ((p(1) - y) * dy)) _ / ((dx * dx) + (dy * dy))) If (t > 1) Then
        x = p2(0)
        y = p2(1)
    ElseIf (t > 0) Then
        x = (x + (dx * t))
        y = (y + (dy * t))
    End If
    
    dx = (p(0) - x)
    dy = (p(1) - y)
    
    sqSegDist = ((dx * dx) + (dy * dy))
End Function
    
Private Function sqSegSegDist(ByVal x0 As Double, _
                              ByVal y0 As Double, _
                              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 Double
' segment to segment distance, ported from http://geomalgorithms.com/a07-_distance.html by Dan Sunday
    
    Dim ux As Double: ux = (x1 - x0)
    Dim uy As Double: uy = (y1 - y0)
    Dim vx As Double: vx = (x3 - x2)
    Dim vy As Double: vy = (y3 - y2)
    Dim wx As Double: wx = (x0 - x2)
    Dim wy As Double: wy = (y0 - y2)
    Dim a As Double: a = ((ux * ux) + (uy * uy))
    Dim b As Double: b = ((ux * vx) + (uy * vy))
    Dim c As Double: c = ((vx * vx) + (vy * vy))
    Dim d As Double: d = ((ux * wx) + (uy * wy))
    Dim e As Double: e = ((vx * wx) + (vy * wy))
    Dim D_ As Double: D_ = ((a * c) - (b * b))
    Dim tN As Double
    Dim sc As Double
    Dim sN As Double
    Dim tc As Double
    Dim sD As Double: sD = D_
    Dim tD As Double: tD = D_
    
    sN = 0
    sD = 1
    tN = e
    tD = c
    sN = ((b * e) - (c * d))
    tN = ((a * e) - (b * d))
    
    If (sN < 0) Then sN = 0 tN = e tD = c ElseIf (sN > sD) Then
        sN = sD
        tN = (e + b)
        tD = c
    End If
    
    If (tN < 0) Then
        tN = 0
        If ((d * -1) < 0) Then sN = 0 ElseIf ((d * -1) > a) Then
            sN = sD
        Else
            sN = (d * -1)
            sD = a
        End If
        
    ElseIf (tN > tD) Then
        tN = tD
        If (((d * -1) + b) < 0) Then sN = 0 ElseIf ((d * -1) + (b > a)) Then
            sN = sD
        Else
            sN = ((d * -1) + b)
            sD = a
        End If
        
    End If

    sc = VBA.IIf(sN = 0, 0, sN / sD)
    tc = VBA.IIf(tN = 0, 0, tN / tD)
    
    Dim cx As Double: cx = (((1 - sc) * x0) + (sc * x1))
    Dim cy As Double: cy = (((1 - sc) * y0) + (sc * y1))
    Dim cx2 As Double: cx2 = (((1 - tc) * x2) + (tc * x3))
    Dim cy2 As Double: cy2 = (((1 - tc) * y2) + (tc * y3))
    Dim dx As Double: dx = (cx2 - cx)
    Dim dy As Double: dy = (cy2 - cy)
    
    sqSegSegDist = ((dx * dx) + (dy * dy))
End Function


For the C++ version (it seems that the C++ source code gets stuck on some test, so watch out):

'https://www.researchgate.net/publication/220868874_Concave_hull_A_k-nearest_neighbours_approach_for_the_computation_of_the_region_occupied_by_a_set_of_points
'https://www.codeproject.com/Articles/1201438/The-Concave-Hull-of-a-Set-of-Points


'Option Explicit
'
'Private Type PointValue
'    point As tXYZ
'    distance As Double
'    Angle As Double
'End Type
'
'Const stride As Long = 24 ' size in bytes of x, y, id
'
'using PointVector = std::vector;
'using PointValueVector = std::vector;
'using LineSegment = std::pair<Point, Point>;
'
''' Floating point comparisons
''Private Function Equal(double a, double b) As boolean;
''Private Function Zero(double a) As boolean;
''Private Function LessThan(double a, double b) As boolean;
''Private Function LessThanOrEqual(double a, double b) As boolean;
''Private Function GreaterThan(double a, double b) As boolean;
'
''' I/O
''Private Function Usage() As void;
''Private Function FindArgument(int argc, char **argv, const std::string &name) As int;
''Private Function ParseArgument(int argc, char **argv, const std::string &name, std::string &val) As int;
''Private Function ParseArgument(int argc, char **argv, const std::string &name, int &val) As int;
''Private Function HasExtension(const std::string &filename, const std::string &ext) As boolean;
''Private Function ReadFile(const std::string &filename, int fieldX = 1, int fieldY = 2) As PointVector;
''Private Function Print(const std::string &filename, const PointVector &points) As void;
''Private Function Print(FILE *out, const PointVector &points, const char *format = "%.3f  %.3f\n") As void;
''Private Function Split(const std::string &value, const char *delims) As std::vector;
'
''' Algorithm-specific
''Private Function NearestNeighboursFlann(flann::Index<flann::L2> &index, const Point &p, size_t k) As PointValueVector;
''Private Function ConcaveHull(PointVector &dataset, size_t k, bool iterate) As PointVector;
''Private Function ConcaveHull(PointVector &dataset, size_t k, PointVector &hull) As boolean;
''Private Function SortByAngle(PointValueVector &values, const Point &p, double prevAngle) As PointVector;
''Private Function AddPoint(PointVector &points, const Point &p) As void;
'
''' General maths
''Private Function PointsEqual(const Point &a, const Point &b) As boolean;
''Private Function Angle(const Point &a, const Point &b) As double;
''Private Function NormaliseAngle(double radians) As double;
''Private Function PointInPolygon(const Point &p, const PointVector &list) As boolean;
''Private Function Intersects(const LineSegment &a, const LineSegment &b) As boolean;
'
''' Point list utilities
''Private Function FindMinYPoint(const PointVector &points) As Point;
''Private Function RemoveDuplicates(PointVector &points) As void;
''Private Function IdentifyPoints(PointVector &points) As void;
''Private Function RemoveHull(PointVector &points, const PointVector &hull) As PointVector::iterator;
''Private Function MultiplePointInPolygon(PointVector::iterator begin, PointVector::iterator end, const PointVector &hull) As boolean;
'
'Private Function main(argc As Integer, argv() As Byte) As Integer
'
'    std::cout << "Concave hull: A k-nearest neighbours approach.\n";
'
'    ' input filename is the only requirement
'    if (argc == 1)
'        Usage();
'        return EXIT_FAILURE;
'    End If
'
'    std::string filename(argv[1]);
'
'    ' The input field numbers for x and y coordinates
'    int fieldX = 1;
'    int fieldY = 2;
'    if (FindArgument(argc, argv, "-field_for_x") != -1) then         ParseArgument(argc, argv, "-field_for_x", fieldX);
'    if (FindArgument(argc, argv, "-field_for_y") != -1) then         ParseArgument(argc, argv, "-field_for_y", fieldY);
'
'    ' Read input
'    PointVector points = ReadFile(filename, fieldX, fieldY);
'    size_t uncleanCount = points.size();
'
'    ' Remove duplicates and id the points
'    RemoveDuplicates(points);
'    size_t cleanCount = points.size();
'    IdentifyPoints(points);
'
'    ' Starting k-value
'    int k = 0;
'    if (FindArgument(argc, argv, "-k") != -1)
'        ParseArgument(argc, argv, "-k", k);
'    k = std::min(std::max(k, 3), (int)points.size() - 1);
'
'    ' For debug purposes, optionally disable iterating k.
'    bool iterate = true;
'    if (FindArgument(argc, argv, "-no_iterate") != -1) then iterate = false;
'
'    std::cout << "Filename         : " << filename << "\n";
'    std::cout << "Input points     : " << uncleanCount << "\n";
'    std::cout << "Input (cleaned)  : " << cleanCount << "\n";
'    std::cout << "Initial 'k'      : " << k << "\n";
'    std::cout << "Final 'k'        : " << k;
'
'    Private Function startTime = std::chrono::high_resolution_clock::now();
'
'    PointVector hull = ConcaveHull(points, (size_t)k, iterate);
'
'    Private Function endTime = std::chrono::high_resolution_clock::now();
'    Private Function duration = std::chrono::duration_cast(endTime - startTime).count();
'
'    std::cout << "\n";
'    std::cout << "Output points    : " << hull.size() << "\n";
'    std::cout << "Time (excl. i/o) : " << std::fixed << std::setprecision(1) << (double)duration / 1000.0 << "s\n";
'    std::cout << "\n";
'
'    ' Optional no further output
'    if (FindArgument(argc, argv, "-no_out") != -1)
'        if (FindArgument(argc, argv, "-out") != -1)
'            std::cout << "Output to file overridden by switch -no_out.\n";
'        return EXIT_SUCCESS;
'        End If
'
'        ' Output to file or stdout
'        if (FindArgument(argc, argv, "-out") != -1)
'            std::string output;
'            ParseArgument(argc, argv, "-out", output);
'
'            Print(output, hull);
'            std::cout << output << " written.\n";
'        End If
'    Else
'        ' Nothing specified, so output to console
'        Print(stdout, hull);
'    End If
'
'    return EXIT_SUCCESS;
'End Function
'
'Private Function Usage() As void
'' Print program usage info.
'    std::cout << "Usage: concave.exe filename [-out arg] [-k arg] [-field_for_x arg] [-field_for_y arg] [-no_out] [-no_iterate]\n";
'    std::cout << "\n";
'    std::cout << " filename      (required) : file of input coordinates, one row per point.\n";
'    std::cout << " -out          (optional) : output file for the hull polygon coordinates. Default=stdout.\n";
'    std::cout << " -k            (optional) : start iteration K value. Default=3.\n";
'    std::cout << " -field_for_x  (optional) : 1-based column number of input for x-coordinate. Default=1.\n";
'    std::cout << " -field_for_y  (optional) : 1-based column number of input for y-coordinate. Default=2.\n";
'    std::cout << " -no_out       (optional) : disable output of the hull polygon coordinates.\n";
'    std::cout << " -no_iterate   (optional) : stop after only one iteration of K, irrespective of result.\n";
'End Function
'
'Private Function FindArgument(int argc, char **argv, const std::string &name) As int
'' Get command line index of name
'    for (int i = 1; i < argc; ++i) ' { ' if (std::string(argv[i]) == name) ' return i; ' } ' return -1; 'End Function ' 'Private Function ParseArgument(int argc, char **argv, const std::string &name, std::string &val) As int '' Get the command line value (string) for name ' int index = FindArgument(argc, argv, name) + 1; ' if (index > 0 && index < argc) ' val = argv[index]; ' ' return index - 1; 'End Function ' 'Private Function ParseArgument(int argc, char **argv, const std::string &name, int &val) As int '' Get the command line value (int) for name ' int index = FindArgument(argc, argv, name) + 1; ' ' if (index > 0 && index < argc) ' val = atoi(argv[index]); ' ' return (index - 1); 'End Function ' 'Private Function HasExtension(const std::string &filename, const std::string &ext) As boolean '' Check whether a string ends with a specified suffix. ' if (filename.length() >= ext.length())
'        return (0 == filename.compare(filename.length() - ext.length(), ext.length(), ext));
'    return false;
'End Function
'
'Private Function ReadFile(const std::string &filename, int fieldX, int fieldY) As PointVector
'' Read a file of coordinates into a vector. First two fields of comma/tab/space delimited input are used.
'    fieldX--; ' from 1-based index to 0-based
'    fieldY--;
'
'    PointVector list;
'    Point p;
'    std::string line;
'    std::vector tokens;
'
'    std::ifstream fin(filename.c_str());
'    if (fin.is_open())
'        {
'        While (fin.good())
'            {
'            getline(fin, line);
'            if (!line.empty())
'                {
'                tokens = Split(line, " ,\t");
'                if (tokens.size() >= 2)
'                    {
'                    p.x = std::atof(tokens[fieldX].c_str());
'                    p.y = std::atof(tokens[fieldY].c_str());
'                    list.push_back(p);
'                    }
'                }
'            }
'        }
'
'    return list;
'End Function
'
'' Output a point list to a file
'Private Function Print(const std::string &filename, const PointVector &points) As void
'    std::string format;
'    if (HasExtension(filename, ".csv"))
'        format = "%.3f,%.3f\n";
'    Else
'        format = "%.3f  %.3f\n";
'
'    FILE *out = fopen(filename.c_str(), "w+");
'    if (out)
'        {
'        Print(out, points, format.c_str());
'
'        fclose(out);
'        }
'End Function
'
'' Output a point list to a stream with a given format string
'Private Function Print(FILE *out, const PointVector &points, const char *format) As void
'    for (const Private Function &p : points)
'        {
'        fprintf(out, format, p.x, p.y);
'        }
'End Function
'
'' Iteratively call the main algorithm with an increasing k until success
'Private Function ConcaveHull(PointVector &dataset, size_t k, bool iterate) As PointVector
'    While (k < dataset.Size())
'        {
'        PointVector hull;
'        if (ConcaveHull(dataset, k, hull) || !iterate)
'            {
'            return hull;
'            }
'        k++;
'        }
'
'    return{};
'End Function
'
'Private Function ConcaveHull(PointVector &pointList, size_t k, PointVector &hull) As boolean
'' The main algorithm from the Moreira-Santos paper.
'    Erase hull()
'
'    if (pointList.size() < 3) then return true
'    If (pointList.Size() = 3) Then
'        hull = pointList
'        return true
'    End If
'
'    ' construct a randomized kd-tree index using 4 kd-trees
'    ' 2 columns, but stride = 24 bytes in width (x, y, ignoring id)
'    flann::Matrix matrix(&(pointList.front().x), pointList.size(), 2, stride);
'    flann::Index<flann::L2> flannIndex(matrix, flann::KDTreeIndexParams(4));
'    flannIndex.buildIndex();
'
'    std::cout << "\rFinal 'k'        : " << k; ' ' ' Initialise hull with the min-y point ' Point firstPoint = FindMinYPoint(pointList); ' AddPoint(hull, firstPoint); ' ' ' Until the hull is of size > 3 we want to ignore the first point from nearest neighbour searches
'    Point currentPoint = firstPoint;
'    flannIndex.removePoint(firstPoint.id);
'
'    double prevAngle = 0.0;
'    int step = 1;
'
'    ' Iterate until we reach the start, or until there's no points left to process
'    while ((!PointsEqual(currentPoint, firstPoint) || step == 1) && hull.size() != pointList.size())
'        {
'        if (step == 4)
'            {
'            ' Put back the first point into the dataset and into the flann index
'            firstPoint.id = pointList.size();
'            flann::Matrix firstPointMatrix(&firstPoint.x, 1, 2, stride);
'            flannIndex.addPoints(firstPointMatrix);
'            }
'
'        PointValueVector kNearestNeighbours = NearestNeighboursFlann(flannIndex, currentPoint, k);
'        PointVector cPoints = SortByAngle(kNearestNeighbours, currentPoint, prevAngle);
'
'        bool its = true;
'        size_t i = 0;
'
'        while (its && i < cPoints.size())
'            {
'            size_t lastPoint = 0;
'            if (PointsEqual(cPoints[i], firstPoint))
'                lastPoint = 1;
'
'            size_t j = 2;
'            its = false;
'
'            while (!its && j < hull.size() - lastPoint)
'                {
'                Private Function line1 = std::make_pair(hull[step - 1], cPoints[i]);
'                Private Function line2 = std::make_pair(hull[step - j - 1], hull[step - j]);
'                its = Intersects(line1, line2);
'                j++;
'                }
'
'            if (its)
'                i++;
'            }
'
'        if (its)
'            return false;
'
'        currentPoint = cPoints[i];
'
'        AddPoint(hull, currentPoint);
'
'        prevAngle = Angle(hull[step], hull[step - 1]);
'
'        flannIndex.removePoint(currentPoint.id);
'
'        step++;
'        }
'
'    ' The original points less the points belonging to the hull need to be fully enclosed by the hull in order to return true.
'    PointVector dataset = pointList;
'
'    Private Function newEnd = RemoveHull(dataset, hull);
'    bool allEnclosed = MultiplePointInPolygon(begin(dataset), newEnd, hull);
'
'    return allEnclosed;
'End Function
'
'Private Function Equal(a As Double, b As Double) As Boolean
'' Compare a and b for equality
'    Equal = (Abs(a - b) <= EPSILON)
'End Function
'
'Private Function Zero(double a) As boolean
'' Compare value to zero
'    return fabs(a) <= DBL_EPSILON;
'End Function
'
'Private Function LessThan(double a, double b) As boolean
'' Compare for a < b
'    return a < (b - DBL_EPSILON);
'End Function
'
'Private Function LessThanOrEqual(double a, double b) As boolean
'' Compare for a <= b
'    return a <= (b + DBL_EPSILON); 'End Function ' 'Private Function GreaterThan(double a, double b) As boolean '' Compare for a > b
'    return a > (b + DBL_EPSILON);
'End Function
'
'Private Function PointsEqual(const Point &a, const Point &b) As boolean
'' Compare whether two points have the same x and y
'    return Equal(a.x, b.x) && Equal(a.y, b.y);
'End Function
'
'Private Function RemoveDuplicates(PointVector &points) As void
'' Remove duplicates in a list of points
'    sort(begin(points), end(points), [](const Point & a, const Point & b)
'        {
'        if (Equal(a.x, b.x))
'            return LessThan(a.y, b.y);
'        Else
'            return LessThan(a.x, b.x);
'        });
'
'    Private Function newEnd = unique(begin(points), end(points), [](const Point & a, const Point & b)
'        {
'        return PointsEqual(a, b);
'        });
'
'    points.erase(newEnd, end(points));
'End Function
'
'Private Function IdentifyPoints(PointVector &points) As void
'' Uniquely id the points for binary searching
'    uint64_t id = 0;
'
'    for (Private Function itr = begin(points); itr != end(points); ++itr, ++id)
'        {
'        itr->id = id;
'        }
'End Function
'
'' Find the point having the smallest y-value
'Private Function FindMinYPoint(const PointVector &points) As Point
'    assert(!points.empty());
'
'    Private Function itr = min_element(begin(points), end(points), [](const Point & a, const Point & b)
'        {
'        if (Equal(a.y, b.y))
'            return GreaterThan(a.x, b.x);
'        Else
'            return LessThan(a.y, b.y);
'        });
'
'    return *itr;
'End Function
'
'' Lookup by ID and remove a point from a list of points
'Private Function RemovePoint(PointVector &list, const Point &p) As void
'    Private Function itr = std::lower_bound(begin(list), end(list), p, [](const Point & a, const Point & b)
'        {
'        return a.id < b.id; ' }); ' ' assert(itr != end(list) && itr->id == p.id);
'
'    if (itr != end(list))
'        list.erase(itr);
'End Function
'
'' Add a point to a list of points
'Private Function AddPoint(PointVector &points, const Point &p) As void
'    points.push_back(p);
'End Function
'
'' Return the k-nearest points in a list of points from the given point p (uses Flann library).
'Private Function NearestNeighboursFlann(flann::Index<flann::L2> &index, const Point &p, size_t k) As PointValueVector
'    std::vector vIndices(k);
'    std::vector vDists(k);
'    double test[] = { p.x, p.y };
'
'    flann::Matrix query(test, 1, 2);
'    flann::Matrix mIndices(vIndices.data(), 1, static_cast(vIndices.size()));
'    flann::Matrix mDists(vDists.data(), 1, static_cast(vDists.size()));
'
'    int count_ = index.knnSearch(query, mIndices, mDists, k, flann::SearchParams(128));
'    size_t count = static_cast(count_);
'
'    PointValueVector result(count);
'
'    for (size_t i = 0; i < count; ++i)
'        {
'        int id = vIndices[i];
'        const double *point = index.getPoint(id);
'        result[i].point.x = point[0];
'        result[i].point.y = point[1];
'        result[i].point.id = id;
'        result[i].distance = vDists[i];
'        }
'
'    return result;
'End Function
'
'' Returns a list of points sorted in descending order of clockwise angle
'Private Function SortByAngle(PointValueVector &values, const Point &from, double prevAngle) As PointVector
'    for_each(begin(values), end(values), [from, prevAngle](PointValue & to)
'        {
'        to.angle = NormaliseAngle(Angle(from, to.point) - prevAngle);
'        });
'
'    sort(begin(values), end(values), [](const PointValue & a, const PointValue & b)
'        {
'        return GreaterThan(a.angle, b.angle);
'        });
'
'    PointVector angled(values.size());
'
'    transform(begin(values), end(values), begin(angled), [](const PointValue & pv)
'        {
'        return pv.point;
'        });
'
'    return angled;
'End Function
'
'' Get the angle in radians measured clockwise from +'ve x-axis
'Private Function Angle(const Point &a, const Point &b) As double
'    double angle = -atan2(b.y - a.y, b.x - a.x);
'
'    return NormaliseAngle(angle);
'End Function
'
'' Return angle in range: 0 <= angle < 2PI
'Private Function NormaliseAngle(double radians) As double
'    if (radians < 0.0)
'        return radians + M_PI + M_PI;
'    Else
'        return radians;
'End Function
'
'' Return the new logical end after removing points from dataset having ids belonging to hull
'Private Function RemoveHull(PointVector &points, const PointVector &hull) As PointVector::iterator
'    std::vector ids(hull.size());
'
'    transform(begin(hull), end(hull), begin(ids), [](const Point & p)
'        {
'        return p.id;
'        });
'
'    sort(begin(ids), end(ids));
'
'    return remove_if(begin(points), end(points), [&ids](const Point & p)
'        {
'        return binary_search(begin(ids), end(ids), p.id);
'        });
'End Function
'
'' Uses OpenMP to determine whether a condition exists in the specified range of elements. https://msdn.microsoft.com/en-us/library/ff521445.aspx
'template 
'bool omp_parallel_any_of(InIt first, InIt last, const Predicate &pr)
'{
'    typedef typename std::iterator_traits::value_type item_type;
'
'    ' A flag that indicates that the condition exists.
'    bool found = false;
'
'    #pragma omp parallel for
'    for (int i = 0; i < static_cast(last - first); ++i)
'        {
'        if (!found)
'            {
'            item_type &cur = *(first + i);
'
'            ' If the element satisfies the condition, set the flag to cancel the operation.
'            if (pr(cur))
'                {
'                found = true;
'                }
'            }
'        }
'
'    return found;
'End Function
'
'' Check whether all points in a begin/end range are inside hull.
'Private Function MultiplePointInPolygon(PointVector::iterator begin, PointVector::iterator end, const PointVector &hull) As boolean
'    Private Function test = [&hull](const Point & p)
'        {
'        return !PointInPolygon(p, hull);
'        };
'
'    bool anyOutside = true;
'
'#if defined USE_OPENMP
'
'    anyOutside = omp_parallel_any_of(begin, end, test); ' multi-threaded
'
'#Else
'
'    anyOutside = std::any_of(begin, end, test); ' single-threaded
'
'#End If
'
'    return !anyOutside;
'End Function
'
'' Point-in-polygon test
'Private Function PointInPolygon(p As tXYZ, list As PointVector) As Boolean
'    If (list.Size() <= 2) Then PointInPolygon = False ' ' const double &x = p.x; ' const double &y = p.y; ' ' int inout = 0; ' Private Function v0 = list.begin(); ' Private Function v1 = v0 + 1; ' ' while (v1 != list.end()) ' { ' if ((LessThanOrEqual(v0->y, y) && LessThan(y, v1->y)) || (LessThanOrEqual(v1->y, y) && LessThan(y, v0->y)))
'            {
'            if (!Zero(v1->y - v0->y))
'                {
'                double tdbl1 = (y - v0->y) / (v1->y - v0->y);
'                double tdbl2 = v1->x - v0->x;
'
'                if (LessThan(x, v0->x + (tdbl2 * tdbl1)))
'                    inout++;
'                }
'            }
'
'        v0 = v1;
'        v1++;
'        }
'
'    If (inout = 0) Then
'        PointInPolygon = False
'    ElseIf (inout Mod 2 = 0) Then
'        PointInPolygon = False
'    Else
'        PointInPolygon = True
'    End If
'End Function
'
'' Test whether two line segments intersect each other
'Private Function Intersects(a As LineSegment, b As LineSegment) As Boolean
'' https://www.topcoder.com/community/data-science/data-science-tutorials/geometry-concepts-line-intersection-and-its-applications/
'
'    const double &ax1 = a.first.x;
'    const double &ay1 = a.first.y;
'    const double &ax2 = a.second.x;
'    const double &ay2 = a.second.y;
'    const double &bx1 = b.first.x;
'    const double &by1 = b.first.y;
'    const double &bx2 = b.second.x;
'    const double &by2 = b.second.y;
'
'    double a1 = ay2 - ay1;
'    double b1 = ax1 - ax2;
'    double c1 = a1 * ax1 + b1 * ay1;
'    double a2 = by2 - by1;
'    double b2 = bx1 - bx2;
'    double c2 = a2 * bx1 + b2 * by1;
'    double det = a1 * b2 - a2 * b1;
'
'    If (Zero(det)) Then
'        Intersects = False
'    Else
'        double x = (b2 * c1 - b1 * c2) / det;
'        double y = (a1 * c2 - a2 * c1) / det;
'
'        bool on_both = true;
'        on_both = on_both && LessThanOrEqual(std::min(ax1, ax2), x) && LessThanOrEqual(x, std::max(ax1, ax2));
'        on_both = on_both && LessThanOrEqual(std::min(ay1, ay2), y) && LessThanOrEqual(y, std::max(ay1, ay2));
'        on_both = on_both && LessThanOrEqual(std::min(bx1, bx2), x) && LessThanOrEqual(x, std::max(bx1, bx2));
'        on_both = on_both && LessThanOrEqual(std::min(by1, by2), y) && LessThanOrEqual(y, std::max(by1, by2));
'        return on_both;
'    End If
'End Function
'
'Private Function ToDegrees(ByVal radians As Double) As Double
'    ToDegrees = radians * 180 / PI
'End Function