强角点检测算子的Pascal实现代码

PRocedure CornerDetect(sWidth, sHeight: longint; Quality: extended);
var
  i, j, fi, fj: longint;
  a, b, c, sum, MinAccept, MaxEigenvalue: extended;
begin
  FeatureCount := 0;
  {
  下面采用Good Feature To Track介绍的方法
  J. Shi and C. Tomasi "Good Features to Track", CVPR 94
  }
  for i := 1 to sWidth - 2 do
    for j := 1 to sHeight - 2 do begin
      dx[i, j] := ImageGray[i - 1, j - 1] + 2 * ImageGray[i - 1, j] + ImageGray[i - 1, j + 1]
        - (ImageGray[i + 1, j - 1] + 2 * ImageGray[i + 1, j] + ImageGray[i + 1, j + 1]);
      dy[i, j] := ImageGray[i - 1, j + 1] + 2 * ImageGray[i, j + 1] + ImageGray[i + 1, j + 1]
        - (ImageGray[i - 1, j - 1] + 2 * ImageGray[i, j - 1] + ImageGray[i + 1, j - 1]);
      dxy[i, j] := ImageGray[i + 1, j - 1] + ImageGray[i - 1, j + 1]
        - (ImageGray[i - 1, j - 1] + ImageGray[i + 1, j + 1]);
    end;
  {求取Sobel算子的Dx, Dy, Dxy
  Dx:
  |1 0 -1|
  |2 0 -2|
  |1 0 -1|
  Dy:
  |-1 -2 -1|
  | 0  0  0|
  | 1  2  1|
  Dxy
  |-1  0  1|
  | 0  0  0|
  | 1  0 -1|}
  MaxEigenvalue := 0;
  for i := 2 to sWidth - 3 do
    for j := 2 to sHeight - 3 do begin
      a := 0; b := 0; c := 0;
      for fi := i - 1 to i + 1 do
        for fj := j - 1 to j + 1 do begin
          a := a + sqr(dx[fi, fj]);
          b := b + dxy[fi, fj];
          c := c + sqr(dy[fi, fj]);
        end;
      a := a / 2; c := c / 2;
      Eigenvalues[i, j] := (a + c - sqrt((a - c) * (a - c) + b * b));
      if Eigenvalues[i, j] > MaxEigenvalue then MaxEigenvalue := Eigenvalues[i, j];
    end;
  {求取矩阵
    |∑Dx*Dx   ∑Dxy|
  M=|               |
    |∑Dxy   ∑Dy*Dy|
  的特征值
  λ= ∑Dx*Dx + ∑Dy*Dy - ((∑Dx*Dx+∑Dy*Dy)^2-4*(∑Dx*Dx * ∑Dy*Dy - ∑Dxy * ∑Dxy))^1/2}
  MinAccept := MaxEigenvalue * Quality;
  {设置最小允许阀值}

  for i := 8 to sWidth - 9 do
    for j := 8 to sHeight - 9 do
      if Eigenvalues[i, j] > MinAccept then begin
        WBPoint[i, j] := true;
        Inc(FeatureCount);
      end else
        WBPoint[i, j] := false;

  for i := 8 to sWidth - 9 do
    for j := 8 to sHeight - 9 do
      if WBPoint[i, j] then begin
        sum := Eigenvalues[i, j];
        for fi := i - 8 to i + 8 do begin
          for fj := j - 8 to j + 8 do
            if sqr(fi - i) + sqr(fj - j) <= 64 then
              if (Eigenvalues[fi, fj] >= sum) and ((fi <> i) or (fj <> j)) and (WBPoint[fi, fj]) then begin
                WBPoint[i, j] := false;
                Dec(FeatureCount);
                break;
              end;
          if not WBPoint[i, j] then break;
        end;
      end;
  {用非最大化抑制来抑制假角点}

  setlength(Features, FeatureCount); fi := 0;
  for i := 8 to sWidth - 9 do
    for j := 8 to sHeight - 9 do
      if WBPoint[i, j] then begin
        Features[fi].Info.X := i;
        Features[fi].Info.Y := j;
        Features[fi].Index := 0;
        Inc(fi);
      end;
  {输出最终的点序列}
end;