Canny Edge

program Canny;
{$I SRL1300/osr.simba}
{$R-}

procedure SplitRGB(a: TMufasaBitmap; out R,G,B: TSingleMatrix);
var
  W,H,x,y: Int32;
begin
  W := a.GetWidth();
  H := a.GetHeight();
  R.SetSize(W,H);
  G.SetSize(W,H);
  B.SetSize(W,H);
  for y:=0 to H-1 do
    for x:=0 to W-1 do begin
      R[y,x] := (a.GetPixel(x,y){shr 00}and $FF);
      G[y,x] := (a.GetPixel(x,y) shr 08 and $FF);
      B[y,x] := (a.GetPixel(x,y) shr 16 and $FF);
    end;
end;

procedure Gradients(a: TMufasaBitmap; out Mag, Dir: TSingleMatrix);
var
  W,H,x,y: Int32;
  dx,dy,theta: Single;
  r,g,b: TSingleMatrix;
const
  Deg22_5:  Single = 0.3926990816987;
  Deg45_0:  Single = 0.7853981633975;
  Deg67_5:  Single = Deg45_0 + Deg22_5;
  Deg90_0:  Single = 1.5707963267949;
  Deg112_5: Single = Deg90_0 + Deg22_5;
  Deg135_0: Single = 2.35619449019234;
  Deg157_5: Single = Deg135_0 + Deg22_5;
begin
  W := a.GetWidth();
  H := a.GetHeight();

  Mag.SetSize(W, H);
  Dir.SetSize(W, H);

  SplitRGB(a, r,g,b);

  for y:=1 to H-2 do
    for x:=1 to W-2 do
    begin
      dx := (-1 * r[y-1,x-1]) + (-2 * r[y-1,x+0]) + (-1 * r[y-1,x+1]) + (+1 * r[y+1,x-1]) + (+2 * r[y+1,x+0]) + (+1 * r[y+1,x+1]);
      dy := (-1 * r[y-1,x-1]) + (+1 * r[y-1,x+1]) + (-2 * r[y+0,x-1]) + (+2 * r[y+0,x+1]) + (-1 * r[y+1,x-1]) + (+1 * r[y+1,x+1]);

      dx += (-1 * g[y-1,x-1]) + (-2 * g[y-1,x+0]) + (-1 * g[y-1,x+1]) + (+1 * g[y+1,x-1]) + (+2 * g[y+1,x+0]) + (+1 * g[y+1,x+1]);
      dy += (-1 * g[y-1,x-1]) + (+1 * g[y-1,x+1]) + (-2 * g[y+0,x-1]) + (+2 * g[y+0,x+1]) + (-1 * g[y+1,x-1]) + (+1 * g[y+1,x+1]);

      dx += (-1 * b[y-1,x-1]) + (-2 * b[y-1,x+0]) + (-1 * b[y-1,x+1]) + (+1 * b[y+1,x-1]) + (+2 * b[y+1,x+0]) + (+1 * b[y+1,x+1]);
      dy += (-1 * b[y-1,x-1]) + (+1 * b[y-1,x+1]) + (-2 * b[y+0,x-1]) + (+2 * b[y+0,x+1]) + (-1 * b[y+1,x-1]) + (+1 * b[y+1,x+1]);

      Mag[y,x] := Sqrt(Sqr(dx * 0.333333334) + Sqr(dy * 0.333333334));
      theta := ArcTan2(dy,dx); // approximate me please
      Dir[y,x] := -1;

      if      ((theta >= -Deg22_5) and (theta <= Deg22_5))  or ((theta <= -Deg157_5) or  (theta >= Deg157_5)) then
        Dir[y,x] := 0
      else if ((theta >= Deg22_5)  and (theta <= Deg67_5))  or ((theta <= -Deg112_5)  and (theta >= -Deg157_5)) then
        Dir[y,x] := 1
      else if ((theta >= Deg67_5)  and (theta <= Deg112_5)) or ((theta <= -Deg67_5)  and (theta >= -Deg112_5)) then
        Dir[y,x] := 2
      else if ((theta >= Deg112_5) and (theta <= Deg157_5)) or ((theta <= -Deg22_5) and (theta >= -Deg67_5)) then
        Dir[y,x] := 3;
    end;
end;

function NonMaximaSuppress(Gm, Gd: TSingleMatrix; low:Single = 1.0): TSingleMatrix;
var
  x,y,dx,dy: Int32;
  mag, theta: Single;
begin
  Result.SetSize(Gm.Width, Gm.Height);

  for y:=1 to Gm.Height()-2 do
    for x:=1 to Gm.Width()-2 do
    begin
      mag := Gm[y,x];
      if mag < low then Continue;
      dx := 0;
      dy := 0;
      case Round(Gd[y,x]) of
        0: if (mag > Gm[y+1,x])   and (mag > Gm[y-1,x]) then
             Result[y,x] := mag;
        1: if (mag > Gm[y-1,x-1]) and (mag > Gm[y+1,x+1]) then
             Result[y,x] := mag;
        2: if (mag > Gm[y,x+1])   and (mag > Gm[y,x-1]) then
             Result[y,x] := mag;
        3: if (mag > Gm[y-1,x+1]) and (mag > Gm[y+1,x-1]) then
             Result[y,x] := mag;
      end;
    end;
end;

function Threshold(Image: TSingleMatrix; low, high: Int32): TSingleMatrix;
var
  i,j,x,y: Int32;
  isValuable, isUnsure: Boolean;
begin
  Result.SetSize(Image.Width, Image.Height);

  for i:=1 to Image.Height-2 do
  begin
    for j:=1 to Image.Width-2 do
    begin
      if(Image[i,j] > high) then
        Result[i,j] := 255
      else if(Image[i,j] < Low) then
        Result[i,j] := 0
      else
      begin
        // figure out if it's low or high, must have high neighbor
        isValuable := False;
        isUnsure   := False;
        for y:=i-1 to i+1 do
        begin
          for x:=j-1 to j+1 do
          begin
            if(Image[y,x] > high) then
            begin
              Result[i,j] := 255;
              isValuable := True;
              break;
            end
            else if (Image[y,x] <= high) and (Image[y,x] >= low) then
              isUnsure := True;
          end;
          if (isValuable) then
            break;
        end;

        if(not isValuable) then
          Result[i,j] := 0;
      end;
    end;
  end;
end;

function CannyEdge(Image: TMufasaBitmap; TLow, THigh: Int32): TSingleMatrix;
var
  Gm, Gd, nms: TSingleMatrix;
begin
  Gradients(Image, Gm, Gd);
  nms := NonMaximaSuppress(Gm, Gd);
  Result := Threshold(nms, TLow, THigh);
end;

var
  bmp, tmp: TMufasaBitmap;
  t: Double;
begin
  bmp.Init(client.GetMBitmaps);
  bmp.LoadFromFile('Images/lenna.jpg');
  tmp.Init(client.GetMBitmaps);

  //bmp.Downsample(2,tmp);
  //bmp.Free();
  //bmp := tmp;

  bmp.Blur(3); bmp.Blur(3);

  t := PerformanceTimer();
  bmp.DrawMatrix(CannyEdge(bmp, 40,70).ToInt);
  WriteLn(PerformanceTimer() - t);

  bmp.Debug();
  bmp.Free();
end.