Simba 1.5 simple concave hull

program new;

function DouglasPeucker(TPA: TPointArray; epsilon:Double): TPointArray;
var
  L, i, index: Int32;
  dmax, d: Double;
  Slice1,Slice2: TPointArray;
begin
  L := Length(TPA);
  if L = 0 then Exit;

  for i:=1 to High(TPA) do
  begin
    d := DistToLine(TPA[i], TPA[0], TPA[High(TPA)]);
    if ( d > dmax ) then
    begin
      index := i;
      dmax  := d;
    end;
  end;

  if (dmax > epsilon) then
  begin
    Slice1 := DouglasPeucker(Copy(TPA, 0, index), epsilon);
    Slice2 := DouglasPeucker(Copy(TPA, index), epsilon);

    Result := Slice1;
    Result += Slice2;
  end else
    Result := [TPA[0], TPA[High(TPA)]];
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 ConcaveHullEx(TPA: TPointArray; MaxLeap: Double=-1; Epsilon:Double=2): T2DPointArray;
var
  res, pts, poly: TPointArray;
  tree: TSlackTree;
  i,c: Int32;
  B: TBox;
begin
  B := TPA.Bounds();
  TPA := TPA.PartitionEx([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(PolygonArea(TPA.ConvexHull()) / TPA.Length())*Sqrt(2));

  MaxLeap := Max(MaxLeap, Epsilon*2);
  SetLength(res, 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;
    pts := pts.ConvexHull().Connect();

    if Length(pts)+c >= Length(res) then
      SetLength(res, Max(c*2, c+Length(pts)));
    Move(pts[0], res[c], Length(pts)*SizeOf(TPoint));

    Inc(c, Length(pts));
  end;

  SetLength(pts, c);
  pts := pts.Unique(); //not needed, but might be faster
  for pts in res.Cluster(2) do
    Result += DouglasPeucker(pts.Border(), Epsilon);
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 ConcaveHull(TPA: TPointArray; Epsilon:Double=2.5; kCount:Int32=5): TPointArray;
var
  res, pts, poly: TPointArray;
  tree: TSlackTree;
  i,c: Int32;
  B: TBox;
begin
  B := TPA.Bounds();
  TPA := TPA.PartitionEx([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);

  SetLength(res, 256);
  for i:=0 to High(tree.data) do
  begin
    pts := tree.KNearest(tree.data[i].split, kCount, False);
    //tree.data[i].hidden := True;

    if Length(pts) <= 1 then continue;
    pts := pts.ConvexHull().Connect();

    if Length(pts)+c >= Length(res) then
      SetLength(res, Max(c*2, c+Length(pts)));
    Move(pts[0], res[c], Length(pts)*SizeOf(TPoint));

    Inc(c, Length(pts));
  end;

  SetLength(res, c);

  Result := DouglasPeucker(res.Border(), Max(2,Epsilon/2));
end;

procedure GenTestHullEx(bmp: TMufasaBitmap; TPA: TPointArray);
var
  i,j,color: Int32;
  t: Double;
  polys: T2DPointArray;
begin
  t := PerformanceTimer();
  polys := ConcaveHullEx(TPA);
  WriteLn PerformanceTimer() - t;

  for i:=0 to High(polys) do
  begin
    color := Random($FFFFFF);
    bmp.DrawTPA(polys[i].Connect(), color);
    for j:=0 to High(polys[i]) do BMP.DrawCircleFilled(polys[i][j], 2, color);
  end;
end;

procedure GenTestHull(bmp: TMufasaBitmap; TPA: TPointArray);
var
  i,j,color: Int32;
  t: Double;
  poly: TPointArray;
begin
  t := PerformanceTimer();
  poly := ConcaveHull(TPA);
  WriteLn PerformanceTimer() - t;

  color := Random($FFFFFF);
  bmp.DrawTPA(poly.Connect(), color);
  for j:=0 to High(poly) do BMP.DrawCircleFilled(poly[j], 2, color);
end;

var
  TPA,pts,poly: TPointArray;
  polys: T2DPointArray;
  bmp: TMufasaBitmap;
  i,d,j,color: Int32;
  t,x: Double;

begin
  bmp := TMufasaBitmap.CreateFromFile('images/shapes.png');
  TPA := bmp.Finder.FindColor(0, 90, [0,0,bmp.GetWidth()-1, bmp.GetHeight()-1]);
  WriteLn(Length(TPA));
  TPA := RandomTPA(10000, [50,50,850,850]);
  bmp.Clear();
  bmp.DrawTPA(TPA, $555555);

  GenTestHullEx(bmp, TPA);
  bmp.Show();

  //GenTestHull(bmp, TPA);
  //bmp.Show();
end.