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

【 tulaoshi.com - 编程语言 】

　　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;