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