Add TPointArray.DistanceTransform

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

@@ -709,6 +709,16 @@ procedure _LapeTPASortCircular(const Params: PParamArray; const Result: Pointer)
   PPointArray(Result)^ := PPointArray(Params^[0])^.SortCircular(PPoint(Params^[1])^, PInteger(Params^[2])^, PBoolean(Params^[3])^);
 end;
 
+(*
+TPointArray.DistanceTransform
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+> function TPointArray.DistanceTransform: TSingleMatrix;
+*)
+procedure _LapeTPADistanceTransform(const Params: PParamArray; const Result: Pointer); LAPE_WRAPPER_CALLING_CONV
+begin
+  PSingleMatrix(Result)^ := PPointArray(Params^[0])^.DistanceTransform();
+end;
+
 procedure ImportTPA(Compiler: TSimbaScript_Compiler);
 begin
   with Compiler do
@@ -798,8 +808,10 @@ procedure ImportTPA(Compiler: TSimbaScript_Compiler);
     addGlobalFunc('function TPointArray.Difference(Other: TPointArray): TPointArray', @_Lape_Point_Difference);
     addGlobalFunc('function TPointArray.SymmetricDifference(Other: TPointArray): TPointArray', @_Lape_Point_SymmetricDifference);
 
+    addGlobalFunc('function TPointArray.DistanceTransform: TSingleMatrix;', @_LapeTPADistanceTransform);
+
     ImportingSection := '';
   end;
 end;
 
-end.
+end.

Source/simba.tpa.pas

@@ -18,6 +18,7 @@
   - Grow
   - RotateEx
   - PartitionEx
+  - DistanceTransform
 }
 
 {
@@ -128,6 +129,8 @@ interface
     function PartitionEx(BoxWidth, BoxHeight: Integer): T2DPointArray; overload;
 
     function Intersection(Other: TPointArray): TPointArray;
+
+    function DistanceTransform: TSingleMatrix;
   end;
 
 implementation
@@ -2125,5 +2128,105 @@ function TPointArrayHelper.Intersection(Other: TPointArray): TPointArray;
   Result := Algo_Point_Intersection(Self, Other);
 end;
 
+function TPointArrayHelper.DistanceTransform: TSingleMatrix;
+
+  function EucDist(const x1,x2:Int32): Int32; inline;
+  begin
+    Result := Sqr(x1) + Sqr(x2);
+  end;
+
+  function EucSep(const i,j, ii,jj:Int32): Int32; inline;
+  begin
+    Result := Round((sqr(j) - sqr(i) + sqr(jj) - sqr(ii))/(2*(j-i)));
+  end;
+
+  function Transform(const binIm:TIntegerArray; m,n:Int32): TSingleMatrix;
+  var
+    x,y,h,w,i,wid:Int32;
+    tmp,s,t:TIntegerArray;
+  begin
+    // first pass
+    SetLength(tmp, m*n);
+    h := n-1;
+    w := m-1;
+    for x:=0 to w do
+    begin
+      if binIm[x] = 0 then
+        tmp[x] := 0
+      else
+        tmp[x] := m+n;
+
+      for y:=1 to h do
+        if (binIm[y*m+x] = 0) then
+          tmp[y*m+x] := 0
+        else
+          tmp[y*m+x] := 1 + tmp[(y-1)*m+x];
+
+      for y:=h-1 downto 0 do
+        if (tmp[(y+1)*m+x] < tmp[y*m+x]) then
+          tmp[y*m+x] := 1 + tmp[(y+1)*m+x]
+    end;
+
+    // second pass
+    SetLength(Result,n,m);
+    SetLength(s,m);
+    SetLength(t,m);
+    wid := 0;
+    for y:=0 to h do
+    begin
+      i := 0;
+      s[0] := 0;
+      t[0] := 0;
+
+      for x:=1 to W do
+      begin
+        while (i >= 0) and (EucDist(t[i]-s[i], tmp[y*m+s[i]]) > EucDist(t[i]-x, tmp[y*m+x])) do
+          Dec(i);
+        if (i < 0) then
+        begin
+          i := 0;
+          s[0] := x;
+        end else
+        begin
+          wid := 1 + EucSep(s[i], x, tmp[y*m+s[i]], tmp[y*m+x]);
+          if (wid < m) then
+          begin
+            Inc(i);
+            s[i] := x;
+            t[i] := wid;
+          end;
+        end;
+      end;
+
+      for x:=W downto 0 do
+      begin
+        Result[y,x] := Sqrt(EucDist(x-s[i], tmp[y*m+s[i]]));
+        if (x = t[i]) then
+          Dec(i);
+      end;
+    end;
+  end;
+
+var
+  Data:TIntegerArray;
+  w,h,n,i:Int32;
+  B:TBox;
+begin
+  Result := nil;
+  if (Length(Self) = 0) then
+    Exit;
+
+  B := Self.Bounds();
+  B.Y1 -= 1;
+  B.X1 -= 1;
+  w := (B.x2 - B.X1) + 2;
+  h := (B.y2 - B.Y1) + 2;
+  SetLength(Data, h*w);
+  for i:=0 to High(Self) do
+    Data[(Self[i].y-B.Y1)*w+(Self[i].x-B.X1)] := 1;
+
+  Result := Transform(data,w,h);
+end;
+
 end.