procedure TOpticalFlowLK.InitFeatures(Frames: TBitmap);
var
  i, j: longint;
  Line: pRGBTriple;
begin
  SetStretchBltMode(Frame.Canvas.Handle, Stretch_Halftone);
  StretchBlt(Frame.Canvas.Handle, 0, 0, Width, Height, Frames.Canvas.Handle, 0, 0, Frames.Width, Frames.Height, SrcCopy);
  for i := 0 to Height - 1 do begin
    Line := Frame.ScanLine[i];
    for j := 0 to Width - 1 do begin
      ImageOld[j, i, 0] := (Line^.rgbtBlue * 11 + Line^.rgbtGreen * 59 + Line^.rgbtRed * 30) div 100;
      ImageGray[j, i] := ImageOld[j, i, 0];
      Inc(Line);
    end;
  end;
  {初始化金字塔图像第一层ImageOld[x,y,0]}
  MakePyramid(ImageOld, Width, Height, L);
  {生成金字塔图像}
  CornerDetect(Width, Height, 0.01);
  {进行强角点检测}
end;

procedure TOpticalFlowLK.CalOpticalFlowLK(Frames: TBitmap);
var
  i, j, fi, fj, k, ll, m, dx, dy, gx, gy, px, py, kx, ky, ed, edc, nWidth, nHeight: longint;
  nx, ny, vx, vy, A, B, C, D, E, F, Ik: extended;
  Ix, Iy: TDoubleExtendedArray;
  Line: pRGBTriple;
  Change: boolean;
begin
  SetStretchBltMode(Frame.Canvas.Handle, Stretch_Halftone);
  StretchBlt(Frame.Canvas.Handle, 0, 0, Width, Height, Frames.Canvas.Handle, 0, 0, Frames.Width, Frames.Height, SrcCopy);
  for i := 0 to Height - 1 do begin
    Line := Frame.ScanLine[i];
    for j := 0 to Width - 1 do begin
      ImageNew[j, i, 0] := (Line^.rgbtBlue * 11 + Line^.rgbtGreen * 59 + Line^.rgbtRed * 30) div 100;
      Inc(Line);
    end;
  end;
  {初始化金字塔图像第一层ImageNew[x,y,0]}
  MakePyramid(ImageNew, Width, Height, L);
  {生成金字塔图像}
  setlength(Ix, 15, 15); setlength(Iy, 15, 15);
  {申请差分图像临时数组}
  for m := 0 to FeatureCount - 1 do begin
    {算法细节见:
    Jean-Yves Bouguet "Pyramidal Implementation of the Lucas Kanade Feature Tracker Description of the algorithm"}
    gx := 0; gy := 0;
    for ll := L - 1 downto 0 do begin
      px := Features[m].Info.X shr ll;
      py := Features[m].Info.Y shr ll;
      {对应当前金字塔图像的u点：u[L]:=u/2^L}
      nWidth := Width shr ll; nHeight := Height shr ll;
      A := 0; B := 0; C := 0;
      for i := max(px - 7, 1) to min(px + 7, nWidth - 2) do
        for j := max(py - 7, 1) to min(py + 7, nHeight - 2) do begin
          fi := i - px + 7; fj := j - py + 7;
          Ix[fi, fj] := (ImageOld[i + 1, j, ll] - ImageOld[i - 1, j, ll]) / 2;
          Iy[fi, fj] := (ImageOld[i, j + 1, ll] - ImageOld[i, j - 1, ll]) / 2;
          A := A + Ix[fi, fj] * Ix[fi, fj]; B := B + Ix[fi, fj] * Iy[fi, fj];
          C := C + Iy[fi, fj] * Iy[fi, fj];
        end;
      {计算2阶矩阵G:
        |Ix(x,y)*Ix(x,y)  Ix(x,y)*Iy(x,y)|
      ∑|Ix(x,y)*Iy(x,y)  Iy(x,y)*Iy(x,y)|}
      D := A * C - B * B;
      vx := 0; vy := 0; dx := 0; dy := 0;
      if abs(D) > 1E-8 then begin
        for k := 1 to 10 do begin
          E := 0; F := 0;
          for i := max(px - 7, 1) to min(px + 7, nWidth - 2) do
            for j := max(py - 7, 1) to min(py + 7, nHeight - 2) do begin
              fi := i - px + 7; fj := j - py + 7;
              kx := i + gx + dx; ky := j + gy + dy;
              if kx < 0 then kx := 0; if kx > nWidth - 1 then kx := nWidth - 1;
              if ky < 0 then ky := 0; if ky > nHeight - 1 then ky := nHeight - 1;
              Ik := ImageOld[i, j, ll] - ImageNew[kx, ky, ll];
              E := E + Ik * Ix[fi, fj];
              F := F + Ik * Iy[fi, fj];
            end;
          {计算2x1阶矩阵b:
            |Ik(x,y)*Ix(x,y)|
          ∑|Ik(x,y)*Iy(x,y)|}
          nx := (C * E - B * F) / D;
          ny := (B * E - A * F) / (-D);
          {计算η=G^-1*b}
          vx := vx + nx; vy := vy + ny;
          dx := trunc(vx); dy := trunc(vy);
          {得到相对运动向量d}
        end;
      end;
      gx := (gx + dx) shl 1; gy := (gy + dy) shl 1;
      {得到下一层的预测运动向量g}
    end;
    gx := gx div 2; gy := gy div 2;
    px := px + gx; py := py + gy;
    Features[m].Info.X := px;
    Features[m].Info.Y := py;
    Features[m].Vector.X := gx;
    Features[m].Vector.Y := gy;
    if (px > Width - 1) or (px < 0) or (py > Height - 1) or (py < 0) then Features[m].Index := 1;
    {失去特征点处理}
  end;

  for k := 0 to L - 1 do begin
    nWidth := Width shr k; nHeight := Height shr k;
    for i := 0 to nWidth - 1 do
      for j := 0 to nHeight - 1 do
        ImageOld[i, j, k] := ImageNew[i, j, k];
  end;
  {复制J到I}
  repeat
 &nbs

上一页  [1] [2] [3]  下一页