{
  FEDT.pas: 高速ユークリッド距離変換
                               Written by maruguu 2005/01/18
  2005/01/24 : RefreshQの場合分けにもれがあったので修正
}
unit FEDT;

interface

uses
  Windows,
  SysUtils,
  Graphics,
  Math;

function FastEuclideanDistanceTransformation(A: TBitmap): TBitmap;

var
  BackGroundColor: Cardinal = $00000000; // 背景色

implementation

type
  TVector2D = record
    x: Integer;
    y: Integer;
  end;

var
  q: array[1..8] of Integer;      // 近傍画素値
  se: array of array of Integer;  // Squared Euclidean Distanceを格納する
  e: array of array of Double;    // 結果を格納する
  r: array of array of TVector2D; // Distance of relative coordinats vector
  q_pos: array[1..8] of TVector2D = ((x: -1; y:  0), // 近傍画素の位置 234
                                     (x: -1; y: -1), //                1 5
                                     (x:  0; y: -1), //                876
                                     (x:  1; y: -1),
                                     (x:  1; y:  0),
                                     (x:  1; y:  1),
                                     (x:  0; y:  1),
                                     (x: -1; y:  1));
  g: array[1..8] of TVector2D = ((x: 1; y: 0), // G(p, q): The difference of relative coordinates
                                 (x: 1; y: 1),
                                 (x: 0; y: 1),
                                 (x: 1; y: 1),
                                 (x: 1; y: 0),
                                 (x: 1; y: 1),
                                 (x: 0; y: 1),
                                 (x: 1; y: 1));

// 背景かそうでないかを判定する。背景でないときはTrueを返す
function IsObject(scanline: pByteArray; x: Integer): Boolean;
var
  red, green, blue: byte;
begin
  red   := ($FF and (BackGroundColor shr 8));
  green := ($FF and (BackGroundColor shr 16));
  blue  := ($FF and (BackGroundColor shr 24));
  if((scanline[x * 4    ] <> blue) or
     (scanline[x * 4 + 1] <> green) or
     (scanline[x * 4 + 2] <> red)) then
    result := True
  else
    result := False;
end;

// 端を考慮してrを返す
function GetRelativeCoordinatesVector(x, y, width, height, qi: Integer): TVector2D;
var
  xx, yy: Integer;
begin
  xx := x + q_pos[qi].x;
  yy := y + q_pos[qi].y;

  if((xx <= 0) or (yy <= 0) or (xx >= width - 1) or (yy >= height - 1))then
  begin
    result.x := 0;
    result.y := 0;
  end
  else
  begin
    result.x := r[xx][yy].x;
    result.y := r[xx][yy].y;
  end;
end;

// h(p, q): The Difference of Squared Euclidean Distances
function SquaredEuclideanDistances(x, y, width, height, qi:Integer): Integer;
var
  v: TVector2D;
begin
  Assert((1 <= qi) and (qi <= 8), 'q[]の引数は1-8の間');

  v := GetRelativeCoordinatesVector(x, y, width, height, qi);
  if((qi = 1) or (qi = 5))then
  begin
    result := 2 * v.x + 1;
  end
  else if((qi = 3) or (qi = 7))then
  begin
    result := 2 * v.y + 1;
  end
  else
  begin
    result := 2 * (v.x + v.y + 1);
  end
end;

// 近傍画素の更新
procedure RefreshQ(x, y, width, height: Integer);
begin
  if((x = 0) and (y = 0))then
  begin
    q[1] := 0;
    q[2] := 0;
    q[3] := 0;
    q[4] := 0;
    q[5] := se[x + 1][y    ];
    q[6] := se[x + 1][y + 1];
    q[7] := se[x    ][y + 1];
    q[8] := 0;
  end
  else if((x = 0) and (y = height - 1))then // 2005/01/24 追加
  begin
    q[1] := 0;
    q[2] := 0;
    q[3] := se[x    ][y - 1];
    q[4] := se[x + 1][y - 1];
    q[5] := se[x + 1][y    ];
    q[6] := 0;
    q[7] := 0;
    q[8] := 0;
  end
  else if((x = width - 1) and (y = 0))then // 2005/01/24 追加
  begin
    q[1] := se[x - 1][y    ];
    q[2] := 0;
    q[3] := 0;
    q[4] := 0;
    q[5] := 0;
    q[6] := 0;
    q[7] := se[x    ][y + 1];
    q[8] := se[x - 1][y + 1];
  end
  else if((x = width - 1) and (y = height - 1))then
  begin
    q[1] := se[x - 1][y    ];
    q[2] := se[x - 1][y - 1];
    q[3] := se[x    ][y - 1];
    q[4] := 0;
    q[5] := 0;
    q[6] := 0;
    q[7] := 0;
    q[8] := 0;
  end
  else if(x = 0)then
  begin
    q[1] := 0;
    q[2] := 0;
    q[3] := se[x    ][y - 1];
    q[4] := se[x + 1][y - 1];
    q[5] := se[x + 1][y    ];
    q[6] := se[x + 1][y + 1];
    q[7] := se[x    ][y + 1];
    q[8] := 0;
  end
  else if(y = 0)then
  begin
    q[1] := se[x - 1][y    ];
    q[2] := 0;
    q[3] := 0;
    q[4] := 0;
    q[5] := se[x + 1][y    ];
    q[6] := se[x + 1][y + 1];
    q[7] := se[x    ][y + 1];
    q[8] := se[x - 1][y + 1];
  end
  else if(x = width - 1)then
  begin
    q[1] := se[x - 1][y    ];
    q[2] := se[x - 1][y - 1];
    q[3] := se[x    ][y - 1];
    q[4] := 0;
    q[5] := 0;
    q[6] := 0;
    q[7] := se[x    ][y + 1];
    q[8] := se[x - 1][y + 1];
  end
  else if(y = height - 1)then
  begin
    q[1] := se[x - 1][y    ];
    q[2] := se[x - 1][y - 1];
    q[3] := se[x    ][y - 1];
    q[4] := se[x + 1][y - 1];
    q[5] := se[x + 1][y    ];
    q[6] := 0;
    q[7] := 0;
    q[8] := 0;
  end
  else
  begin
    q[1] := se[x - 1][y    ];
    q[2] := se[x - 1][y - 1];
    q[3] := se[x    ][y - 1];
    q[4] := se[x + 1][y - 1];
    q[5] := se[x + 1][y    ];
    q[6] := se[x + 1][y + 1];
    q[7] := se[x    ][y + 1];
    q[8] := se[x - 1][y + 1];
  end;
end;

procedure ForwardRasterScan(A: TBitmap);
var
  i, x, y, tmp, minq: Integer;
  v: TVector2D;
begin
  for y := 0 to A.Height - 1 do
    for x := 0 to A.Width - 1 do
    begin

      if(IsObject(A.ScanLine[y], x))then
      begin
        RefreshQ(x, y, A.Width, A.Height);
        se[x][y] := $7FFFFFFF;
        minq := 0;
        for i := 1 to 4 do
        begin
          tmp := q[i] + SquaredEuclideanDistances(x, y, A.Width, A.Height, i);
          if(tmp < se[x][y])then
          begin
            minq := i;
            se[x][y] := tmp;
          end;
        end;

        if(minq > 0)then
        begin
          v := GetRelativeCoordinatesVector(x, y, A.Width, A.Height, minq);
          r[x][y].x := v.x + g[minq].x;
          r[x][y].y := v.y + g[minq].y;
        end;
      end;
    end;
end;

procedure BackwardRasterScan(A: TBitmap);
var
  i, x, y, tmp, minq: Integer;
  v: TVector2D;
begin
  for y := A.Height - 1 downto 0 do
    for x := A.Width - 1 downto 0 do
    begin

      if(IsObject(A.ScanLine[y], x))then
      begin
        RefreshQ(x, y, A.Width, A.Height);
        minq := 0;
        for i := 5 to 8 do
        begin
          tmp := q[i] + SquaredEuclideanDistances(x, y, A.Width, A.Height, i);
          if(tmp < se[x][y])then
          begin
            minq := i;
            se[x][y] := tmp;
          end;
        end;

        if(minq > 0)then
        begin
          v := GetRelativeCoordinatesVector(x, y, A.Width, A.Height, minq);
          r[x][y].x := v.x + g[minq].x;
          r[x][y].y := v.y + g[minq].y;
        end;
        e[x][y] := sqrt(se[x][y]);
      end;
    end;
end;

procedure MakeBitmap(var A: TBitmap);
var
  i, x, y: Integer;
  b: byte;
  Max: Double;
  pA: pByteArray;
begin
  Max := 0;
  for y := 0 to A.Height - 1 do
    for x := 0 to A.Width - 1 do
      if(Max < e[x][y])then
        Max := e[x][y];
  for y := 0 to A.Height - 1 do
    for x := 0 to A.Width - 1 do
    begin
      b := Floor((e[x][y] / Max) * 255);
      pA := A.ScanLine[y];
      for i := 0 to 3 do
        pA^[x * 4 + i] := b;
    end;
end;

function FastEuclideanDistanceTransformation(A: TBitmap): TBitmap;
var
  i, j: Integer;
begin
  Assert((A.PixelFormat = PF32BIT), '引数のTBitmapはPF32BIT形式で');

  // 初期化
  SetLength(r, A.Width);
  SetLength(e, A.Width);
  SetLength(se, A.Width);

  for i := Low(r) to High(r) do
  begin
    SetLength(r[i], A.Height);
    SetLength(e[i], A.Height);
    SetLength(se[i], A.Height);
    ZeroMemory(r[i], Sizeof(r[i]));
    ZeroMemory(e[i], Sizeof(e[i]));
    ZeroMemory(se[i], Sizeof(se[i]));
  end;

  Result := TBitmap.Create;
  Result.Width := A.Width;
  Result.Height := A.Height;
  Result.PixelFormat := PF32BIT;

  ForwardRasterScan(A);
  BackwardRasterScan(A);
  MakeBitmap(Result);
end;

end.