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

核心提示：PRocedure CornerDetect(sWidth, sHeight: longint; Quality: extended);var i, j, fi, fj: longint; a, b, c, sum, MinAccept, MaxEigenvalue: extended;begin Feature

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;