{
作者 : 刘留
参考文献为 :
Jean-Yves Bouguet "Implementasi Piramidal dari Deskripsi Pelacak Fitur Lucas Kanade dari Algoritma"
http://www.aivisoft.net/
Geo.cra [at] gmail [dot] com
}
Unit OpticalFlowlk;
antarmuka
penggunaan
Matematika, Windows, Sysutils, Varian, Kelas, Grafik, Unitypes, ColorConV;
jenis
Topticalflowlk = kelas
Pribadi
ImageOld, Imagenew: ttriplelongintarray;
ImageGray, DX, DY, DXY: tdoublelongintarray;
Nilai eigen: tdoubleextEdeRray;
Wbpoint: tdoublebooleanArray;
Tinggi, Lebar, L, Radiox, Radioy: Longint;
Prosedur CornerDetect (Swidth, Sheight: Longint; Kualitas: Diperpanjang);
prosedur makePyramid (gambar var: ttriplelongintarray; swidth, sheight, sl: longint);
publik
Bingkai: tbitmap;
Fitur: TsinglePointInfoarray;
FeatureCount, SupportCount: Longint;
Kecepatan, Radio: Diperpanjang;
Prosedur init (Swidth, Sheight, SL: Longint);
Prosedur initfeatures (frame: tbitmap);
prosedur calopticalflowlk (bingkai: tbitmap);
Destructor menghancurkan; mengesampingkan;
akhir;
pelaksanaan
Prosedur Topticalflowlk.cornerdetect (Swidth, Sheight: Longint; kualitas: diperpanjang);
var
i, j, fi, fj: longint;
a, b, c, jumlah, minak, maxeigenvalue: diperpanjang;
mulai
FeatureCount: = 0;
{
下面采用 Fitur bagus untuk melacak 介绍的方法
J. Shi dan C. Tomasi "Fitur Baik Untuk Dilacak", CVPR 94
}
untuk i: = 1 ke swidth - 2 lakukan
Untuk J: = 1 ke Sheight - 2 Mulai
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]);
akhir;
{求取 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;
untuk i: = 2 to swidth - 3 do
untuk j: = 2 ke sheight - 3 mulai
A: = 0; B: = 0; C: = 0;
untuk fi: = i - 1 hingga i + 1 lakukan
untuk fj: = j - 1 ke j + 1 mulai
A: = a + sqr (dx [fi, fj]);
B: = B + DXY [FI, FJ];
c: = c + sqr (dy [fi, fj]);
akhir;
A: = A / 2; C: = C / 2;
EigenValues [i, j]: = (a + c - sqrt ((a - c) * (a - c) + b * b));
jika nilai eigen [i, j]> maxeigenvalue maka maxeigenvalue: = nilai eigen [i, j];
akhir;
{求取矩阵
| ∑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 }
Minyak: = MaxeigenValue * kualitas;
{设置最小允许阀值}
untuk i: = 8 to swidth - 9 do
untuk j: = 8 ke sheight - 9 do
Jika nilai eigen [i, j]> minak, mulailah
Wbpoint [i, j]: = true;
Inc (featureCount);
akhir yang lain
Wbpoint [i, j]: = false;
untuk i: = 8 to swidth - 9 do
untuk j: = 8 ke sheight - 9 do
Jika wbpoint [i, j] maka mulailah
jumlah: = nilai eigen [i, j];
untuk fi: = i - 8 hingga i + 8 Mulai
untuk fj: = j - 8 ke j + 8 do
jika sqr (fi - i) + sqr (fj - j) <= 64 Kemudian
if (nilai eigen [fi, fj]> = sum) dan ((fi <> i) atau (fj <> j)) dan (wbpoint [fi, fj]) lalu mulai
Wbpoint [i, j]: = false;
Dec (featureCount);
merusak;
akhir;
Jika tidak wbpoint [i, j] maka hancurkan;
akhir;
akhir;
{用非最大化抑制来抑制假角点}
setLength (fitur, featureCount); fi: = 0;
untuk i: = 8 to swidth - 9 do
untuk j: = 8 ke sheight - 9 do
Jika wbpoint [i, j] maka mulailah
Fitur [fi] .info.x: = i;
Fitur [fi] .info.y: = j;
Fitur [fi] .index: = 0;
Inc (FI);
akhir;
{输出最终的点序列}
akhir;
Prosedur Topticalflowlk.init (Swidth, Sheight, SL: Longint);
mulai
Lebar: = Swidth; Tinggi: = Sheight; L: = SL;
setLength (ImageOld, Width, Height, L);
setLength (imagenew, lebar, tinggi, l);
Bingkai: = tbitmap.create;
Frame.width: = lebar; Frame.height: = tinggi;
Frame.pixelformat: = pf24bit;
setLength (gambargray, lebar, tinggi);
setLength (nilai eigen, lebar, tinggi);
setLength (dx, lebar, tinggi);
setLength (dy, lebar, tinggi);
setLength (dxy, lebar, tinggi);
setLength (wbpoint, lebar, tinggi);
FeatureCount: = 0;
akhir;
Prosedur Topticalflowlk.makepyramid (gambar var: ttriplelongintarray; swidth, sheight, sl: longint);
var
i, j, k, ii, jj, nwidth, nheight, owidth, oheight: longint;
mulai
{生成金字塔图像}
Owidth: = Swidth; OHEIGHT: = Sheight;
untuk k: = 1 ke sl - 1 mulai
nwidth: = (owidth + 1) shr 1; nheight: = (oheight + 1) shr 1;
untuk i: = 1 ke nwidth - 2 mulai
II: = I SHL 1;
untuk j: = 1 hingga nheight - 2 Mulai
JJ: = J shl 1;
Gambar [i, j, k]: = (gambar [ii, jj, k - 1] shl 2 +
Gambar [II - 1, JJ, K - 1] SHL 1 + Gambar [II + 1, JJ, K - 1] SHL 1 + Gambar [II, JJ - 1, K - 1] SHL 1 + Gambar [II, JJ + 1, k - 1] shl 1 +
Gambar [II - 1, JJ - 1, K - 1] + Gambar [II + 1, JJ - 1, K - 1] + Gambar [ii - 1, jj + 1, k - 1] + gambar [ii + 1 , jj + 1, k - 1]) shr 4;
{高斯原则 , shl 右移位 , shr 左移位}
akhir;
akhir;
untuk i: = 1 ke nwidth - 2 mulai
II: = I SHL 1;
Gambar [i, 0, k]: = (gambar [ii, 0, k - 1] shl 2 +
Gambar [II - 1, 0, K - 1] SHL 1 + Gambar [II + 1, 0, K - 1] SHL 1 + Gambar [II, 0, K - 1] SHL 1 + Gambar [ii, 1, k - 1] SHL 1 +
Gambar [II - 1, 0, K - 1] + Gambar [ii + 1, 0, k - 1] + Gambar [ii - 1, 1, k - 1] + gambar [ii + 1, 1, k - 1 ]) shr 4;
Gambar [i, nheight - 1, k]: = (gambar [ii, oheight - 1, k - 1] shl 2 +
Gambar [II - 1, OHEIGHT - 1, K - 1] SHL 1 + Gambar [II + 1, OHEIGHT - 1, K - 1] SHL 1 + Gambar [II, OHEIGHT - 2, K - 1] SHL 1 + Gambar [II, OHEIGHT - 1, K - 1] SHL 1 +
Gambar [II - 1, OHEIGHT - 2, K - 1] + Gambar [II + 1, OHEIGHT - 2, K - 1] + Gambar [II - 1, OHEIGHT - 1, K - 1] + Gambar [II + 1 , oheight - 1, k - 1]) shr 4;
akhir;
{处理上下边}
untuk j: = 1 hingga nheight - 2 Mulai
JJ: = J shl 1;
Gambar [0, j, k]: = (gambar [0, jj, k - 1] shl 2 +
Gambar [0, JJ, K - 1] SHL 1 + Gambar [1, JJ, K - 1] SHL 1 + Gambar [0, JJ - 1, K - 1] SHL 1 + Gambar [0, jj + 1, k - 1] SHL 1 +
Gambar [0, JJ - 1, K - 1] + Gambar [1, jj - 1, k - 1] + gambar [0, jj + 1, k - 1] + gambar [1, jj + 1, k - 1 ]) shr 4;
Gambar [nwidth - 1, j, k]: = (gambar [owidth - 1, jj, k - 1] shl 2 +
Gambar [owidth - 2, jj, k - 1] shl 1 + gambar [owidth - 1, jj, k - 1] shl 1 + gambar [owidth - 1, jj - 1, k - 1] shl 1 + gambar [owidth - 1, jj + 1, k - 1] shl 1 +
Gambar [owidth - 2, jj - 1, k - 1] + gambar [owidth - 1, jj - 1, k - 1] + gambar [owidth - 2, jj + 1, k - 1] + gambar [owidth - 1 , jj + 1, k - 1]) shr 4;
akhir;
{处理左右边}
Gambar [0, 0, k]: = (gambar [0, 0, k - 1] shl 2 +
Gambar [0, 0, k - 1] SHL 1 + Gambar [1, 0, K - 1] SHL 1 + Gambar [0, 0, K - 1] SHL 1 + Gambar [0, 1, K - 1] SHL 1 +
Gambar [0, 0, k - 1] + gambar [1, 0, k - 1] + gambar [0, 1, k - 1] + gambar [1, 1, k - 1]) shr 4;
{处理左上点}
Gambar [0, nheight - 1, k]: = (gambar [0, oheight - 1, k - 1] shl 2 +
Gambar [0, OHEIGHT - 1, K - 1] SHL 1 + Gambar [1, OHEIGHT - 1, K - 1] SHL 1 + Gambar [0, OHEIGHT - 2, K - 1] SHL 1 + Gambar [0, OHEIGHT - 1, K - 1] SHL 1 +
Gambar [0, OHEIGHT - 2, K - 1] + Gambar [1, OHEIGHT - 2, K - 1] + Gambar [0, OHEIGHT - 1, K - 1] + Gambar [1, OHEIGHT - 1, K - 1 ]) shr 4;
{处理左下点}
Gambar [nwidth - 1, 0, k]: = (gambar [owidth - 1, 0, k - 1] shl 2 +
Gambar [owidth - 2, oheight - 1, k - 1] shl 1 + gambar [owidth - 1, oheight - 1, k - 1] shl 1 + gambar [owidth - 1, oheight - 1, k - 1] shl 1 + Gambar [OWIDTH - 1, OHEIGHT - 1, K - 1] SHL 1 +
Gambar [owidth - 2, oheight - 1, k - 1] + gambar [owidth - 1, oheight - 1, k - 1] + gambar [owidth - 2, oheight - 1, k - 1] + gambar [owidth - 1 , oheight - 1, k - 1]) shr 4;
{处理右上点}
Gambar [nwidth - 1, nheight - 1, k]: = (gambar [owidth - 1, oheight - 1, k - 1] shl 2 +
Gambar [owidth - 2, oheight - 1, k - 1] shl 1 + gambar [owidth - 1, oheight - 1, k - 1] shl 1 + gambar [owidth - 1, oheight - 2, k - 1] shl 1 + Gambar [OWIDTH - 1, OHEIGHT - 1, K - 1] SHL 1 +
Gambar [owidth - 2, oheight - 2, k - 1] + gambar [owidth - 1, oheight - 2, k - 1] + gambar [owidth - 2, oheight - 1, k - 1] + gambar [owidth - 1 , oheight - 1, k - 1]) shr 4;
{处理右下点}
akhir;
akhir;
Prosedur Topticalflowlk.Initfeatures (frame: tbitmap);
var
I, J: Longint;
Baris: prgbtriple;
mulai
SetStretchBlTmode (frame.canvas.handle, stretch_halftone);
Stretchblt (frame.canvas.handle, 0, 0, lebar, tinggi, frames.canvas.handle, 0, 0, frames.width, frames.height, srccopy);
untuk i: = 0 hingga tinggi - 1 Mulai
Baris: = frame.scanline [i];
untuk j: = 0 ke lebar - 1 mulai
ImageOld [j, i, 0]: = (baris^.rgbtBlue * 11 + baris^.rgbtgreen * 59 + baris^.rgbtred * 30) div 100;
ImageGray [j, i]: = ImageOld [j, i, 0];
Lereng);
akhir;
akhir;
{初始化金字塔图像第一层 ImageOld [x, y, 0]}
MakePyramid (ImageOld, Width, Height, L);
{生成金字塔图像}
CornerDetect (lebar, tinggi, 0,01);
{进行强角点检测}
akhir;
Prosedur Topticalflowlk.CalopticalFlowlk (frame: 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: tdoubleextEdeRarray;
Baris: prgbtriple;
Ubah: Boolean;
mulai
SetStretchBlTmode (frame.canvas.handle, stretch_halftone);
Stretchblt (frame.canvas.handle, 0, 0, lebar, tinggi, frames.canvas.handle, 0, 0, frames.width, frames.height, srccopy);
untuk i: = 0 hingga tinggi - 1 Mulai
Baris: = frame.scanline [i];
untuk j: = 0 ke lebar - 1 mulai
Imagenew [j, i, 0]: = (baris^.rgbtBlue * 11 + baris^.rgbtgreen * 59 + baris^.rgbtred * 30) div 100;
Lereng);
akhir;
akhir;
{初始化金字塔图像第一层 imagenew [x, y, 0]}
MakePyramid (imagenew, lebar, tinggi, l);
{生成金字塔图像}
setLength (ix, 15, 15); setLength (iy, 15, 15);
{申请差分图像临时数组}
untuk m: = 0 ke featureCount - 1 do begin
{算法细节见:
Jean-Yves Bouguet "Implementasi Piramidal dari Deskripsi Pelacak Fitur Lucas Kanade dari Algoritma"}
GX: = 0; Gy: = 0;
untuk ll: = l - 1 downto 0 Mulai
px: = fitur [m] .info.x shr ll;
py: = fitur [m] .info.y shr ll;
{对应当前金字塔图像的 u 点 : u [l]: = u/2^l}
nwidth: = lebar shr ll; nheight: = tinggi shr ll;
A: = 0; B: = 0; C: = 0;
untuk i: = maks (px - 7, 1) hingga min (px + 7, nwidth - 2) lakukan
untuk j: = maks (py - 7, 1) hingga min (py + 7, nheight - 2) Mulai
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];
akhir;
{计算 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;
Jika abs (d)> 1e-8 maka mulailah
untuk k: = 1 hingga 10 mulai
E: = 0; F: = 0;
untuk i: = maks (px - 7, 1) hingga min (px + 7, nwidth - 2) lakukan
untuk j: = maks (py - 7, 1) hingga min (py + 7, nheight - 2) Mulai
fi: = i - px + 7; fj: = j - py + 7;
kx: = i + gx + dx; KY: = j + gy + dy;
jika kx <0 maka kx: = 0; jika kx> nwidth - 1 lalu kx: = nwidth - 1;
Jika Ky <0 maka Ky: = 0; jika ky> nheight - 1 lalu ky: = nheight - 1;
IK: = ImageOld [i, j, ll] - imagenew [kx, ky, ll];
E: = e + ik * ix [fi, fj];
F: = f + ik * iy [fi, fj];
akhir;
{计算 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}
akhir;
akhir;
gx: = (gx + dx) shl 1; gy: = (gy + dy) shl 1;
{得到下一层的预测运动向量 g}
akhir;
GX: = GX Div 2; gy: = gy div 2;
px: = px + gx; py: = py + gy;
Fitur [m] .info.x: = px;
Fitur [m] .info.y: = py;
Fitur [m] .vector.x: = gx;
Fitur [m] .vector.y: = gy;
if (px> lebar - 1) atau (px <0) atau (py> tinggi - 1) atau (py <0) kemudian fitur [m] .index: = 1;
{失去特征点处理}
akhir;
untuk k: = 0 ke l - 1 mulai
nwidth: = lebar shr k; nheight: = tinggi shr k;
untuk i: = 0 ke nwidth - 1 lakukan
untuk j: = 0 hingga nheight - 1 lakukan
ImageOld [i, j, k]: = imagenew [i, j, k];
akhir;
{复制 j 到 i}
mengulang
Ubah: = false;
untuk i: = 0 ke featureCount - 1 do begin
Jika fitur [i] .index = 1 maka
untuk j: = i + 1 ke featureCount - 1 do begin
Fitur [j - 1]: = fitur [j];
Ubah: = true;
akhir;
Jika perubahan maka hancurkan;
akhir;
jika perubahan maka dec (featureCount);
sampai tidak berubah;
setLength (fitur, featureCount);
{删除失去的特征点}
Kecepatan: = 0; Radio: = 0; Radiox: = 0; RADIOY: = 0;
Jika featureCount> 0 lalu mulai
SupportCount: = 0;
untuk i: = 0 untuk featureCount - 1 do
if (fitur [i] .vector.x <> 0) dan (fitur [i] .vector.y <> 0) Kemudian mulailah
Radiox: = Radiox + Fitur [i] .vector.x;
Radioy: = Radioy + Fitur [i] .vector.y;
Kecepatan: = kecepatan + sqrt (sqr (fitur [i] .vector.x) + sqr (fitur [i] .vector.y));
Inc (SupportCount);
akhir;
Jika SupportCount> 0 maka mulailah
Kecepatan: = kecepatan / supportCount; //*0.0352778;
Radio: = arctan2 (Radioy, Radiox);
akhir;
akhir;
{计算平均速度和整体方向}
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 1;
Frame.canvas.pen.color: = cllime;
Frame.canvas.brush.style: = bsclear;
untuk i: = 0 untuk featureCount - 1 do
Frame.canvas.ellipse (fitur [i] .info.x - 4, fitur [i] .info.y - 4, fitur [i] .info.x + 4, fitur [i] .info.y + 4) ;
{用绿圈标识特征点}
Frame.canvas.pen.color: = clyellow;
untuk i: = 0 ke featureCount - 1 do begin
Frame.canvas.moveto (fitur [i] .info.x - fitur [i] .vector.x, fitur [i] .info.y - fitur [i] .vector.y);
Frame.canvas.lineto (fitur [i] .info.x, fitur [i] .info.y);
akhir;
{用黄色线条表示运动向量}
if (supportCount> 0) dan (kecepatan> 3) lalu mulailah
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 6;
Frame.canvas.pen.color: = clwhite;
Frame.canvas.ellipse (lebar div 2 - 100, tinggi div 2 - 100, lebar div 2 + 100, tinggi div 2 + 100);
Frame.canvas.moveto (lebar div 2, tinggi div 2);
Frame.canvas.pen.color: = clblue;
Frame.canvas.lineto (lebar div 2 + trunc (90 * cos (radio)), tinggi div 2 + trunc (90 * sin (radio)));
Frame.canvas.pen.style: = psclear;
Frame.canvas.brush.style: = bssolid;
Frame.canvas.brush.color: = clred;
Frame.canvas.ellipse (lebar div 2 + trunc (90 * cos (radio)) - 8, tinggi div 2 + trunc (90 * sin (radio)) - 8, lebar div 2 + trunc (90 * cos (radio) ) + 8, tinggi div 2 + trunc (90 * sin (radio)) + 8);
akhir;
akhir;
Destructor Topticalflowlk.Destroy;
mulai
setLength (ImageOld, 0);
setLength (imagenew, 0);
setLength (ImageGray, 0);
setLength (nilai eigen, 0);
setLength (dx, 0);
setLength (dy, 0);
setLength (dxy, 0);
setLength (wbpoint, 0);
diwariskan;
akhir;
akhir.