[Simba1300] Canny & Hough

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 CannyThreshold(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 := CannyThreshold(nms, TLow, THigh);
end;

function Hough(Edges: TSingleMatrix; ntx, mry: Int32): TSingleMatrix;
var
  H,W, x,y, ts,iry: Int32;
  rmax, dr,dth, th,r, col: Single;
begin
  W := Edges.Width;
  H := Edges.Height;
  mry := (mry div 2) * 2;

  Result.SetSize(ntx, mry);

  rmax := Hypot(W,H);
  dr  := rmax / (mry / 2);
  dth := PI / ntx;

  for y:=0 to H-1 do
    for x:=0 to W-1 do
    begin
      if Edges[y,x] = 0 then
        continue;

      for ts:=0 to ntx-1 do
      begin
        th  := ts * dth;
        r   := (x - w/2)*Cos(th) + (y - h/2)*Sin(th);
        Result[Trunc(r + mry / 2), ts] += 1;
      end;
    end;
end;


function HoughLine(acc: TSingleMatrix; v: TPoint; iw,ih: Int32): TBox;
var
  t,r,tdeg: Double;
begin
  tdeg := v.x / 2;
  t := Radians(tdeg);
  r := v.y;

  with Result do
  begin
    r := (r - acc.Height / 2);
      if (tdeg >= 45) and (tdeg <= 135) then
    begin
          x1 := 0;
          y1 := (ih div 2) + Trunc((r - (x1 - iw/2) * Cos(t)) / Sin(t));
          x2 := iw;
          y2 := (ih div 2) + Trunc((r - (x2 - iw/2) * Cos(t)) / Sin(t));
      end else
      begin
          y1 := 0;
          x1 := (iw div 2) + Trunc((r - (y1 - ih/2) * Sin(t)) / Cos(t));
          y2 := ih;
          x2 := (iw div 2) + Trunc((r - (y2 - ih/2) * Sin(t)) / Cos(t));
    end;
  end;
end;

function Test(bmp: TMufasaBitmap; thresh: Double): TBoxArray;
var
  i: Int32;
  lp: TPointArray;
  acc,edges: TSingleMatrix;
  t: Double;
begin
  edges := CannyEdge(bmp, 35,50);
  acc   := Hough(edges, 360, Trunc(Sqrt(2) * Max(bmp.GetWidth, bmp.GetHeight)));

  lp := acc.Indices(acc.Max*thresh, __GT__);
  lp := lp.ToATPA(3).Means();

  for i:=0 to High(lp) do
  begin
    Result += HoughLine(acc, lp[i], bmp.GetWidth, bmp.GetHeight);
  end;
end;


var
  bmp, tmp: TMufasaBitmap;
  i,W,H: Int32;
  TPA: TPointArray;
  Lines: TBoxArray;
begin
  bmp.Init(client.GetMBitmaps);
  bmp.LoadFromFile('Images/chess.jpg');

  (*
  tmp.Init(client.GetMBitmaps);
  bmp.Downsample(2,tmp);
  bmp.Free();
  bmp := tmp;
  *)

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

  Lines := Test(bmp, 0.40);
  for i:=0 to High(Lines) do
  begin
    TPA := TPAFromLine(Lines[i].x1, Lines[i].y1, Lines[i].x2, Lines[i].y2);
    TPA := TPA + TPA.OffsetFunc(Point(1,0));
    FilterPointsBox(TPA, 0,0, bmp.GetWidth-1, bmp.GetHeight-1);
    bmp.DrawTPA(TPA, $00FFFF);
  end;

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