DistanceTransform

{*
 Distance transform
*}

type EDistUnit = (duEuclidean, duManhatten, duChebyshev);

function dtEucDist(x1,x2:Int32): Int32;
begin
  Result := Sqr(x1) + Sqr(x2);
end;

function dtEucSep(i,j, ii,jj:Int32): Int32;
begin
  Result := Round((sqr(j) - sqr(i) + sqr(jj) - sqr(ii))/(2*(j-i)));
end;


function dtMtnDist(x1, x2: Int32): Int32;
begin
  Result := Abs(x1) + x2
end;

function dtMtnSep(i,j, ii, jj: Int32): Int32;
begin
  if (jj >= (ii + j - i)) then
    Exit($FFFFFF);
  if (ii > (jj + j - i)) then
    Exit(-$FFFFFF);
  Result := (jj-ii+j+i) shr 1;
end;

function dtChbDist(x1, x2: Int32): Int32;
begin
  Result := Max(Abs(x1), x2);
end;


function dtChbSep(i,j, ii,jj: Int32): Int32;
begin
  if (ii <= jj) then
    Result := Max(i+jj, Trunc((i+j) / 2))
  else
    Result := Min(j-ii, Trunc((i+j) / 2));
end;


function DistanceTransform(const binIm:TIntegerArray; m,n:Int32; distanceUnit:EDistUnit): TSingleMatrix;
type
  TSepFunc = function(i,j, ii,jj: Int32): Int32;
  TDistFunc = function(x1,x2: Int32): Int32;
var
  x,y,h,w,i,wid:Int32;
  dist: single;
  tmp,s,t:TIntegerArray;
  sf: TSepFunc;
  df: TDistFunc;
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;


  case distanceUnit of
    duEuclidean:
    begin
      df := @dtEucDist;
      sf := @dtEucSep;
    end;
    duManhatten:
    begin
      df := @dtMtnDist;
      sf := @dtMtnSep;
    end;
    duChebyshev:
    begin
      df := @dtChbDist;
      sf := @dtChbSep;
    end;
  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 (df(t[i]-s[i], tmp[y*m+s[i]]) > df(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 + sf(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
      dist := df(x-s[i], tmp[y*m+s[i]]);
      if distanceUnit = duEuclidean then dist := Sqrt(dist);

      Result[y,x] := dist;
      if (x = t[i]) then Dec(i);
    end;
  end;
end;


function DistanceTransform(const TPA:TPointArray; distanceUnit:EDistUnit): TSingleMatrix; overload;
var
  data:TIntegerArray;
  w,h,n,i:Int32;
  area:TBox;
begin
  n := Length(TPA);
  if (n = 0) then Exit;
  area := GetTPABounds(TPA);
  area.y1 -= 1;
  area.x1 -= 1;
  w := (area.x2 - area.x1) + 2;
  h := (area.y2 - area.y1) + 2;
  SetLength(data, h*w);
  for i:=0 to n-1 do
    data[(TPA[i].y-area.y1)*w+(TPA[i].x-area.x1)] := 1;

  Result := DistanceTransform(data,w,h,distanceUnit);
end;