{
作者 ： 刘留
：
Jean-Yves Bouguet "Implémentation pyramidale de la description du tracker de fonctionnalités Lucas Kanade de l'algorithme"
http://www.aivisoft.net/
Geo.cra [at] gmail [dot] com
}
unité OpticalFlowlk;
interface
usages
Mathématiques, fenêtres, synutils, variantes, classes, graphiques, unitypes, colorconv;
taper
Topticalflowlk = classe
Privé
Imageold, ImageNew: TtriplelongIntArray;
ImageGray, DX, Dy, DXY: TdoublelongIntArray;
Values Eigen: TdoubleExtendedArray;
WBPoint: TdoubleboolEanArray;
Hauteur, largeur, l, radiox, radioy: longint;
Procédure Cornerdetect (Swidth, Sheight: Longint; Qualité: Extended);
Procédure MakePyramid (Var Images: TtriplelongIntArray; Swidth, Sheight, Sl: Longint);
publique
Cadre: Tbitmap;
Caractéristiques: TsinglePointInfoArray;
FeatureCount, SupportCount: Longint;
Vitesse, radio: étendu;
Procédure init (Swidth, Sheight, SL: Longint);
Procédure InitFeures (cadres: tbitmap);
Procédure CalopticalFlowlk (cadres: tbitmap);
destructeur détruire; outrepasser;
fin;
mise en œuvre
procédure topticalflowlk.cornerdetect (serpente, sheight: longint; qualité: étendue);
var
I, J, Fi, FJ: Longint;
a, b, c, sum, minaccept, maxeigenvalue: étendu;
commencer
FeatureCount: = 0;
{
下面采用 Bonne fonctionnalité pour suivre 介绍的方法
J. Shi et C. Tomasi "Bonnes fonctionnalités à suivre", CVPR 94
}
pour i: = 1 à la peau - 2 do
pour j: = 1 à sheight - 2 commencez
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]);
fin;
{求取 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;
pour i: = 2 à la peau - 3 do
pour j: = 2 à sheight - 3 commencez
a: = 0; b: = 0; C: = 0;
pour fi: = i - 1 à i + 1 do
pour fj: = j - 1 à j + 1 commencez
a: = a + sqr (dx [fi, fj]);
b: = b + dxy [fi, fj];
c: = c + sqr (dy [fi, fj]);
fin;
a: = a / 2; c: = c / 2;
Les valeurs eigen [i, j]: = (a + c - sqrt ((a - c) * (a - c) + b * b));
Si les valeurs propres [i, j]> maxeigenvalue alors maxeigenvalue: = eigenvalues [i, j];
fin;
{求取矩阵
| ∑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 * Qualité;
{设置最小允许阀值}
pour i: = 8 à la peau - 9 faire
pour j: = 8 à sheight - 9 faire
Si les valeurs propres [i, j]> minaccept alors commencez
Wbpoint [i, j]: = true;
Inc (featCount);
fin d'autre
Wbpoint [i, j]: = false;
pour i: = 8 à la peau - 9 faire
pour j: = 8 à sheight - 9 faire
Si wbpoint [i, j] alors commencez
sum: = eigenvalues [i, j];
pour fi: = i - 8 à i + 8, commencez
pour fj: = j - 8 à j + 8 do
Si sqr (fi - i) + sqr (fj - j) <= 64 alors
si (valeurs propres [fi, fj]> = sum) et ((fi <> i) ou (fj <> j)) et (wbpoint [fi, fj]) alors commencez
Wbpoint [i, j]: = false;
Dec (FeatureCount);
casser;
fin;
Si ce n'est pas wbpoint [i, j] alors cassez;
fin;
fin;
{用非最大化抑制来抑制假角点}
SetLength (fonctionnalités, FeatureCount); fi: = 0;
pour i: = 8 à la peau - 9 faire
pour j: = 8 à sheight - 9 faire
Si wbpoint [i, j] alors commencez
Fonctionnalités [fi] .info.x: = i;
Caractéristiques [fi] .info.y: = j;
Fonctionnalités [fi] .Index: = 0;
Inc (fi);
fin;
{输出最终的点序列}
fin;
procédure topticalflowlk.init (Swidth, Sheight, Sl: Longint);
commencer
Largeur: = Swidth; Hauteur: = sheight; L: = SL;
setLength (imageold, largeur, hauteur, l);
SetLength (ImageNew, largeur, hauteur, l);
Cadre: = tbitmap.create;
Frame.width: = largeur; Frame.Height: = hauteur;
Frame.Pixelformat: = pf24bit;
setLength (imagegray, largeur, hauteur);
SetLength (valeurs propres, largeur, hauteur);
setLength (dx, largeur, hauteur);
setLength (dy, largeur, hauteur);
setLength (dxy, largeur, hauteur);
setLength (wbpoint, largeur, hauteur);
FeatureCount: = 0;
fin;
procédure topticalflowlk.makepyramid (var images: ttriplelongInTarray; Swidth, sheight, sl: longint);
var
I, J, K, II, JJ, nwidth, nheight, owidth, oheight: longint;
commencer
{生成金字塔图像}
Owidth: = Swidth; oheight: = sheight;
pour k: = 1 à sl - 1 commence
nwidth: = (owidth + 1) shR 1; nheight: = (oheight + 1) shR 1;
pour i: = 1 à nwidth - 2 commencez
ii: = i shl 1;
pour j: = 1 à nheight - 2 commencez
JJ: = J Shl 1;
Images [i, j, k]: = (images [ii, jj, k - 1] shl 2 +
Images [ii - 1, jj, k - 1] shl 1 + images [ii + 1, jj, k - 1] shl 1 + images [ii, jj - 1, k - 1] shl 1 + images [ii, jj + 1, k - 1] shl 1 +
Images [ii - 1, jj - 1, k - 1] + images [ii + 1, jj - 1, k - 1] + images [ii - 1, jj + 1, k - 1] + images [ii + 1 , JJ + 1, K - 1]) SHR 4;
{高斯原则 ， Shl 右移位 ， SHR 左移位}
fin;
fin;
pour i: = 1 à nwidth - 2 commencez
ii: = i shl 1;
Images [i, 0, k]: = (images [ii, 0, k - 1] shl 2 +
Images [ii - 1, 0, k - 1] shl 1 + images [ii + 1, 0, k - 1] shl 1 + images [ii, 0, k - 1] shl 1 + images [ii, 1, k - 1] shl 1 +
Images [ii - 1, 0, k - 1] + images [ii + 1, 0, k - 1] + images [ii - 1, 1, k - 1] + images [ii + 1, 1, k - 1 ]) SHR 4;
Images [i, nheight - 1, k]: = (images [ii, oheight - 1, k - 1] shl 2 +
Images [II - 1, Oheight - 1, K - 1] SHL 1 + Images [II + 1, Oheight - 1, K - 1] SHL 1 + Images [II, Oheight - 2, K - 1] SHL 1 + [ii, oheight - 1, k - 1] shl 1 +
Images [ii - 1, oheight - 2, k - 1] + images [ii + 1, oheight - 2, k - 1] + images [ii - 1, oheight - 1, k - 1] + images [ii + 1 , Oheight - 1, K - 1]) SHR 4;
fin;
{处理上下边}
pour j: = 1 à nheight - 2 commencez
JJ: = J Shl 1;
Images [0, j, k]: = (images [0, jj, k - 1] shl 2 +
Images [0, jj, k - 1] shl 1 + images [1, jj, k - 1] shl 1 + images [0, jj - 1, k - 1] shl 1 + images [0, jj + 1, k - 1] shl 1 +
Images [0, jj - 1, k - 1] + images [1, jj - 1, k - 1] + images [0, jj + 1, k - 1] + images [1, jj + 1, k - 1 ]) SHR 4;
Images [nwidth - 1, j, k]: = (images [owidth - 1, jj, k - 1] shl 2 +
Images [owidth - 2, jj, k - 1] shl 1 + images [owidth - 1, jj, k - 1] shl 1 + images [owidth - 1, jj - 1, k - 1] shl 1 + images [owidth - 1, jj + 1, k - 1] shl 1 +
Images [owidth - 2, jj - 1, k - 1] + images [owidth - 1, jj - 1, k - 1] + images [owidth - 2, jj + 1, k - 1] + images [owidth - 1 , JJ + 1, K - 1]) SHR 4;
fin;
{处理左右边}
Images [0, 0, k]: = (images [0, 0, k - 1] shl 2 +
Images [0, 0, k - 1] shl 1 + images [1, 0, k - 1] shl 1 + images [0, 0, k - 1] shl 1 + images [0, 1, k - 1] shl 1 +
Images [0, 0, k - 1] + images [1, 0, k - 1] + images [0, 1, k - 1] + images [1, 1, k - 1]) shR 4;
{处理左上点}
Images [0, nheight - 1, k]: = (images [0, oheight - 1, k - 1] shl 2 +
Images [0, oheight - 1, k - 1] shl 1 + images [1, oheight - 1, k - 1] shl 1 + images [0, oheight - 2, k - 1] shl 1 + images [0, oheight - 1, k - 1] shl 1 +
Images [0, oheight - 2, k - 1] + images [1, oheight - 2, k - 1] + images [0, oheight - 1, k - 1] + images [1, oheight - 1, k - 1 ]) SHR 4;
{处理左下点}
Images [nwidth - 1, 0, k]: = (images [owidth - 1, 0, k - 1] shl 2 +
Images [Owidth - 2, Oheight - 1, K - 1] Shl 1 + Images [Owidth - 1, Oheight - 1, K - 1] Shl 1 + Images [Owidth - 1, Oheight - 1, K - 1] Shl 1 + Images [owidth - 1, oheight - 1, k - 1] shl 1 +
Images [Owidth - 2, Oheight - 1, K - 1] + Images [Owidth - 1, Oheight - 1, K - 1] + Images [Owidth - 2, Oheight - 1, K - 1] + Images [Owidth - 1 , Oheight - 1, K - 1]) SHR 4;
{处理右上点}
Images [nwidth - 1, nheight - 1, k]: = (images [owidth - 1, oheight - 1, k - 1] shl 2 +
Images [Owidth - 2, Oheight - 1, K - 1] Shl 1 + Images [Owidth - 1, Oheight - 1, K - 1] Shl 1 + Images [Owidth - 1, Oheight - 2, K - 1] Shl 1 + Images [owidth - 1, oheight - 1, k - 1] shl 1 +
Images [Owidth - 2, Oheight - 2, K - 1] + Images [Owidth - 1, Oheight - 2, K - 1] + Images [Owidth - 2, Oheight - 1, K - 1] + Images [Owidth - 1 , Oheight - 1, K - 1]) SHR 4;
{处理右下点}
fin;
fin;
procédure topticalflowlk.initfeatures (cadres: tbitmap);
var
I, J: Longint;
Ligne: prgbtriple;
commencer
SetStRetchBlTMode (trame.canvas.handle, stretch_halftone);
Stretchblt (trame.canvas.handle, 0, 0, largeur, hauteur, cadres.canvas.handle, 0, 0, frames.width, frames.height, srccopy);
pour i: = 0 à la hauteur - 1 commencez
Ligne: = frame.scanline [i];
Pour j: = 0 à la largeur - 1 commence
Imageold [J, i, 0]: = (Line ^ .rgBtBlue * 11 + Line ^ .rgbtGreen * 59 + Line ^ .rgbTred * 30) div 100;
ImageGray [j, i]: = imageold [j, i, 0];
Inclinaison);
fin;
fin;
{初始化金字塔图像第一层 Imageold [x, y, 0]}
MakePyramid (image, largeur, hauteur, l);
{生成金字塔图像}
Cornerdetect (largeur, hauteur, 0,01);
{进行强角点检测}
fin;
procédure topticalflowlk.calopticalflowlk (cadres: 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: Extension;
Ix, iy: tdoubleextendArray;
Ligne: prgbtriple;
Changement: booléen;
commencer
SetStRetchBlTMode (trame.canvas.handle, stretch_halftone);
Stretchblt (trame.canvas.handle, 0, 0, largeur, hauteur, cadres.canvas.handle, 0, 0, frames.width, frames.height, srccopy);
pour i: = 0 à la hauteur - 1 commencez
Ligne: = frame.scanline [i];
Pour j: = 0 à la largeur - 1 commence
ImageNew [J, i, 0]: = (Line ^ .rgBtBlue * 11 + Line ^ .rgbtGreen * 59 + Line ^ .rgbTred * 30) div 100;
Inclinaison);
fin;
fin;
{初始化金字塔图像第一层 ImageNew [x, y, 0]}
Makepyramide (Imagenew, largeur, hauteur, l);
{生成金字塔图像}
setLength (ix, 15, 15); setLength (iy, 15, 15);
{申请差分图像临时数组}
Pour m: = 0 à FeatureCount - 1 commencez
{算法细节见:
Jean-Yves Bouguet "Implémentation pyramidale du tracker de fonctionnalité Lucas Kanade Description de l'algorithme"}
gx: = 0; gy: = 0;
pour ll: = l - 1 down à 0
px: = fonctionnalités [m] .info.x shr ll;
py: = caractéristiques [m] .info.y shr ll;
{对应当前金字塔图像的 u 点 ： u [l]: = u / 2 ^ l}
nwidth: = largeur shR ll; nheight: = hauteur shr ll;
A: = 0; B: = 0; C: = 0;
pour i: = max (px - 7, 1) à min (px + 7, nwidth - 2) faire
pour j: = max (py - 7, 1) à min (py + 7, nheight - 2) commencez
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];
fin;
{计算 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;
Si ABS (D)> 1E-8 alors commencez
pour k: = 1 à 10 commence
E: = 0; F: = 0;
pour i: = max (px - 7, 1) à min (px + 7, nwidth - 2) faire
pour j: = max (py - 7, 1) à min (py + 7, nheight - 2) commencez
fi: = i - px + 7; fj: = J - py + 7;
kx: = i + gx + dx; ky: = j + gy + dy;
Si kx <0 alors kx: = 0; Si kx> nwidth - 1 alors kx: = nwidth - 1;
Si ky <0 alors ky: = 0; Si ky> nheight - 1 alors ky: = nheight - 1;
Ik: = imageold [i, j, ll] - ImageNew [kx, ky, ll];
E: = e + ik * ix [fi, fj];
F: = f + ik * iy [fi, fj];
fin;
{计算 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}
fin;
fin;
gx: = (gx + dx) shl 1; gy: = (gy + dy) shl 1;
{得到下一层的预测运动向量 g}
fin;
gx: = gx div 2; gy: = gy div 2;
px: = px + gx; py: = py + gy;
Fonctionnalités [m] .info.x: = px;
Caractéristiques [m] .info.y: = py;
Fonctionnalités [m] .vector.x: = gx;
Caractéristiques [m] .vector.y: = gy;
if (px> largeur - 1) ou (px <0) ou (py> hauteur - 1) ou (py <0) alors dispose [m] .index: = 1;
{失去特征点处理}
fin;
pour k: = 0 à l - 1 commencez
nwidth: = largeur shR K; nheight: = hauteur shr k;
pour i: = 0 à nwidth - 1 do
pour j: = 0 à nheight - 1 do
Imageold [i, j, k]: = ImageNew [i, j, k];
fin;
{复制 j 到 i}
répéter
Changement: = false;
pour i: = 0 à FeatureCount - 1 commencez
Si les fonctionnalités [i] .index = 1 alors
pour j: = i + 1 à FeatureCount - 1 commencez
Fonctionnalités [J - 1]: = fonctionnalités [J];
Changement: = true;
fin;
Si changez, cassez;
fin;
Si changez alors dec (FeatureCount);
jusqu'à ce que je ne change pas;
SetLength (fonctionnalités, FeatureCount);
{删除失去的特征点}
Vitesse: = 0; Radio: = 0; Radiox: = 0; Radioy: = 0;
Si featureCount> 0 alors commencez
SupportCount: = 0;
pour i: = 0 à FeatureCount - 1 do
if (fonctionnalités [i] .vector.x <> 0) et (fonctionnalités [i] .vector.y <> 0) alors commencez
RadioX: = RadioX + Fonctions [i] .vector.x;
Radioy: = radioy + caractéristiques [i] .vector.y;
Vitesse: = vitesse + sqrt (sqr (fonctionnalités [i] .vector.x) + sqr (fonctionnalités [i] .vector.y));
Inc (supportCount);
fin;
Si supportCount> 0 alors commencez
Vitesse: = vitesse / supportCount; //*0.0352778;
Radio: = arctan2 (radioy, radiox);
fin;
fin;
{计算平均速度和整体方向}
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 1;
Frame.canvas.pen.color: = cllime;
Frame.canvas.brush.style: = bsclear;
pour i: = 0 à FeatureCount - 1 do
Frame.canvas.ellipse (fonctionnalités [i] .info.x - 4, fonctionnalités [i] .info.y - 4, fonctionnalités [i] .info.x + 4, fonctionnalités [i] .info.y + 4) ;
{用绿圈标识特征点}
Frame.canvas.pen.color: = Clyellow;
pour i: = 0 à FeatureCount - 1 commencez
Frame.canvas.moveto (fonctionnalités [i] .info.x - fonctionnalités [i] .vector.x, fonctionnalités [i] .info.y - fonctionnalités [i] .vector.y);
Frame.canvas.lineto (fonctionnalités [i] .info.x, fonctionnalités [i] .info.y);
fin;
{用黄色线条表示运动向量}
if (supportCount> 0) et (vitesse> 3) alors commencez
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 6;
Frame.canvas.pen.color: = clwhite;
Frame.Canvas.Ellipse (largeur div 2 - 100, hauteur div 2 - 100, largeur div 2 + 100, hauteur div 2 + 100);
Frame.canvas.moveto (largeur div 2, hauteur div 2);
Frame.canvas.pen.color: = clblue;
Frame.canvas.lineto (largeur div 2 + trunc (90 * cos (radio)), hauteur div 2 + truncd (90 * sin (radio)));
Frame.canvas.pen.style: = psClear;
Frame.canvas.brush.style: = bssolid;
Frame.Canvas.brush.color: = Clire;
Frame.canvas.ellipse (largeur div 2 + trunc (90 * cos (radio)) - 8, hauteur div 2 + trunc (90 * sin (radio)) - 8, largeur div 2 + trunc (90 * cos (radio) ) + 8, hauteur div 2 + trunc (90 * sin (radio)) + 8);
fin;
fin;
destructor topticalflowlk.destroy;
commencer
setLength (imageold, 0);
SetLength (ImageNew, 0);
setLength (imagegray, 0);
SetLength (Values Eigen, 0);
setLength (dx, 0);
setLength (dy, 0);
setLength (dxy, 0);
setLength (wbpoint, 0);
hérité;
fin;
fin.