add @slackydev's ConcaveHull

Source/script/imports/simba/simba.import_tpa.pas

@@ -754,6 +754,36 @@ procedure _LapeTPACircularity(const Params: PParamArray; const Result: Pointer);
   PDouble(Result)^ := PPointArray(Params^[0])^.Circularity();
 end;
 
+(*
+TPointArray.DouglasPeucker
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+> function TPointArray.DouglasPeucker(epsilon: Double): TPointArray;
+*)
+procedure _LapeTPADouglasPeucker(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV
+begin
+  PPointArray(Result)^ := PPointArray(Params^[0])^.DouglasPeucker(PDouble(Params^[1])^);
+end;
+
+(*
+TPointArray.ConcaveHull
+~~~~~~~~~~~~~~~~~~~~~~~
+> function TPointArray.ConcaveHull(Epsilon:Double=2.5; kCount:Int32=5): TPointArray;
+*)
+procedure _LapeTPAConcaveHull(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV
+begin
+  PPointArray(Result)^ := PPointArray(Params^[0])^.ConcaveHull(PDouble(Params^[1])^, PInteger(Params^[2])^);
+end;
+
+(*
+TPointArray.ConcaveHullEx
+~~~~~~~~~~~~~~~~~~~~~~~~~
+> function TPointArray.ConcaveHullEx(MaxLeap: Double=-1; Epsilon:Double=2): T2DPointArray;
+*)
+procedure _LapeTPAConcaveHullEx(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV
+begin
+  P2DPointArray(Result)^ := PPointArray(Params^[0])^.ConcaveHullEx(PDouble(Params^[1])^, PDouble(Params^[2])^);
+end;
+
 procedure ImportTPA(Compiler: TSimbaScript_Compiler);
 begin
   with Compiler do
@@ -849,6 +879,10 @@ procedure ImportTPA(Compiler: TSimbaScript_Compiler);
 
     addGlobalFunc('function TPointArray.Circularity: Double;', @_LapeTPACircularity);
 
+    addGlobalFunc('function TPointArray.DouglasPeucker(Epsilon: Double): TPointArray;', @_LapeTPADouglasPeucker);
+    addGlobalFunc('function TPointArray.ConcaveHull(Epsilon: Double = 2.5; kCount: Integer = 5): TPointArray;', @_LapeTPAConcaveHull);
+    addGlobalFunc('function TPointArray.ConcaveHullEx(MaxLeap: Double = -1; Epsilon: Double = 2): T2DPointArray;', @_LapeTPAConcaveHullEx);
+
     ImportingSection := '';
   end;
 end;

Source/simba.tpa.pas

@@ -11,6 +11,7 @@
   - CreateFromEllipse
   - CreateFromCircle
   - ConvexHull
+  - ConcaveHull
   - Border
   - Skeleton
   - MinAreaRect
@@ -136,13 +137,17 @@ interface
     function DistanceTransform: TSingleMatrix;
 
     function Circularity: Double;
+
+    function DouglasPeucker(epsilon: Double): TPointArray;
+    function ConcaveHull(Epsilon:Double=2.5; kCount:Int32=5): TPointArray;
+    function ConcaveHullEx(MaxLeap: Double=-1; Epsilon:Double=2): T2DPointArray;
   end;
 
 implementation
 
 uses
   Math,
-  simba.arraybuffer, simba.geometry, simba.math,
+  simba.atpa, simba.arraybuffer, simba.geometry, simba.math,
   simba.algo_sort, simba.algo_intersection, simba.slacktree, simba.algo_unique;
 
 procedure GetAdjacent4(var Adj: TPointArray; const P: TPoint); inline;
@@ -2417,5 +2422,125 @@ function TPointArrayHelper.MinAreaCircle: TCircle;
     Result.Radius := Ceil(Distance(p1,p2) / 2);
 end;
 
+function TPointArrayHelper.DouglasPeucker(epsilon: Double): TPointArray;
+var
+  H, i, index: Int32;
+  dmax, d: Double;
+  Slice1,Slice2: TPointArray;
+begin
+  H := High(Self);
+  if (H = -1) then
+    Exit(nil);
+
+  dmax := 0;
+  index := 0;
+  for i:=1 to H do
+  begin
+    d := TSimbaGeometry.DistToLine(Self[i], Self[0], Self[H]);
+    if (d > dmax) then
+    begin
+      index := i;
+      dmax  := d;
+    end;
+  end;
+
+  if (dmax > epsilon) then
+  begin
+    Slice1 := Copy(Self, 0, index).DouglasPeucker(epsilon);
+    Slice2 := Copy(Self, index).DouglasPeucker(epsilon);
+
+    Result := Slice1;
+    Result += Slice2;
+  end else
+    Result := [Self[0], Self[High(Self)]];
+end;
+
+(*
+  Concave hull approximation using k nearest neighbors
+  Instead of describing a specific max distance we assume that the boundary points are evenly spread out
+  so we can simply extract a number of neighbors and connect the hull of those.
+  Worst case it cuts off points.
+
+  Will reduce the TPA to a simpler shape if it's dense, defined by epsilon.
+
+  If areas are cut off, you have two options based on your needs:
+  1. Increase "Epsilon", this will reduce accurate.. But it's faster.
+  2. Increase "kCount", this will maintain accuracy.. But it's slower.
+*)
+function TPointArrayHelper.ConcaveHull(Epsilon:Double=2.5; kCount:Int32=5): TPointArray;
+var
+  TPA, pts: TPointArray;
+  Buffer: TSimbaPointBuffer;
+  tree: TSlackTree;
+  i: Int32;
+  B: TBox;
+begin
+  B := Self.Bounds();
+  TPA := Self.PartitionEx(TPoint.Create(B.X1-Round(Epsilon), B.Y1-Round(Epsilon)), Round(Epsilon*2-1), Round(Epsilon*2-1)).Means();
+  if Length(TPA) <= 2 then
+    Exit(TPA);
+
+  tree.Init(TPA);
+  Buffer.Init(256);
+  for i:=0 to High(tree.data) do
+  begin
+    pts := tree.KNearest(tree.data[i].split, kCount, False);
+    if Length(pts) <= 1 then
+      Continue;
+    pts := pts.ConvexHull().Connect();
+
+    Buffer.Add(pts);
+  end;
+
+  TPA := Buffer.ToArray(False);
+
+  Result := TPA.Border().DouglasPeucker(Max(2, Epsilon/2));
+end;
+
+(*
+  Concave hull approximation using range query based on given distance "MaxLeap".
+  if maxleap doesn't cover all of the input then several output polygons will be created.
+  MaxLeap is by default automatically calulcated by the density of the polygon
+  described by convexhull. But can be changed.
+
+  Higher maxleap is slower.
+  Epsilon describes how accurate you want your output, and have some impact on speed.
+*)
+function TPointArrayHelper.ConcaveHullEx(MaxLeap: Double=-1; Epsilon:Double=2): T2DPointArray;
+var
+  TPA, pts: TPointArray;
+  tree: TSlackTree;
+  i: Int32;
+  B: TBox;
+  Buffer: TSimbaPointBuffer;
+  BufferResult: TSimbaPointArrayBuffer;
+begin
+  B := Self.Bounds();
+  TPA := Self.PartitionEx(TPoint.Create(B.X1-Round(Epsilon), B.Y1-Round(Epsilon)), Round(Epsilon*2-1), Round(Epsilon*2-1)).Means();
+  if Length(TPA) <= 2 then
+    Exit([TPA]);
+  tree.Init(TPA);
+
+  if MaxLeap = -1 then
+    MaxLeap := Ceil(Sqrt(TSimbaGeometry.PolygonArea(TPA.ConvexHull()) / Length(TPA)) * Sqrt(2));
+
+  MaxLeap := Max(MaxLeap, Epsilon*2);
+  Buffer.Init(256);
+  for i:=0 to High(tree.data) do
+  begin
+    pts := tree.RangeQueryEx(tree.data[i].split, MaxLeap,MaxLeap, False);
+    if Length(pts) <= 1 then
+      Continue;
+
+    Buffer.Add(pts.ConvexHull().Connect());
+  end;
+
+  pts := Buffer.ToArray(False);
+  for pts in pts.Cluster(2) do
+    BufferResult.Add(pts.Border().DouglasPeucker(Epsilon));
+
+  Result := BufferResult.ToArray();
+end;
+
 end.