{
作者 ： 刘留
参考文献为 ：
Jean-Yves Bouguet "Implementación piramidal de Lucas Kanade Tracker Descripción del algoritmo"
http://www.aivisoft.net/
Geo.cra [at] gmail [dot] com
}
Unidad OpticalFlowlk;
interfaz
usos
Matemáticas, Windows, Sysutils, Variantes, Clases, Gráficos, Unitypes, ColorConv;
tipo
Topticalflowlk = clase
Privado
Imageold, ImageNew: ttriplelongintaRray;
ImageGray, DX, DY, DXY: TDOBLELGINTINTArray;
Valores propios: TDoulextedenseArray;
WBPoint: TDoubleBooleanArray;
Altura, ancho, L, Radiox, Radioy: Longint;
Procedimiento CornerDetect (Swidth, Sheight: Longint; Calidad: Extendida);
Procedimiento MakePyramid (VAR Imágenes: TtriplelongInTintaRray; Swidth, Sheight, SL: Longint);
público
Marco: TBITMAP;
Características: TsinglePointInfoArray;
FeatureCount, SupportCount: Longint;
Velocidad, radio: extendido;
procedimiento init (swidth, sheight, sl: longint);
Procedimiento InitFeatures (marcos: TBITMAP);
procedimiento Calopticalflowlk (marcos: TBITMAP);
destructor destruir; anular;
fin;
implementación
procedimiento topticalflowlk.cornerDetect (Swidth, Sheight: Longint; Calidad: extendido);
varilla
I, J, FI, FJ: Longint;
a, b, c, suma, minacept, valiente máximo: extendido;
comenzar
FeatRecount: = 0;
{
下面采用 Buena característica para rastrear 介绍的方法
J. Shi y C. Tomasi "Buenas características para la pista", CVPR 94
}
para i: = 1 a swidth - 2 do
para j: = 1 a sheight - 2 comienza
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;
para i: = 2 a swidth - 3 do
para j: = 2 a sheight - 3 empiezan
a: = 0; b: = 0; c: = 0;
para fi: = i - 1 a i + 1 hacer
para fj: = j - 1 a j + 1 comenzar
a: = a + sqr (dx [fi, fj]);
b: = b + dxy [fi, fj];
c: = c + sqr (dy [fi, fj]);
fin;
a: = a / 2; c: = c / 2;
EigenValues [i, j]: = (a + c - sqrt ((a - c) * (a - c) + b * b));
Si los valores propios [i, j]> maxeigenvalue entonces 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 * calidad;
{设置最小允许阀值}
para i: = 8 a swidth - 9 do
para j: = 8 a sheight - 9 do
Si los valores propios [i, j]> minaccept, entonces comience
Wbpoint [i, j]: = true;
Inc (FeatRecount);
fin
Wbpoint [i, j]: = false;
para i: = 8 a swidth - 9 do
para j: = 8 a sheight - 9 do
Si wbpoint [i, j] entonces comience
suma: = valores propios [i, j];
para fi: = i - 8 a i + 8 comience
para fj: = j - 8 a j + 8 hacer
Si sqr (fi - i) + sqr (fj - j) <= 64 entonces
if (eigenValues [fi, fj]> = sum) y ((fi <> i) o (fj <> j)) y (wbpoint [fi, fj]) luego comience
Wbpoint [i, j]: = false;
DEC (FeatRecount);
romper;
fin;
Si no es wbpoint [i, j], entonces rompa;
fin;
fin;
{用非最大化抑制来抑制假角点}
SetLength (características, FeatRecount); fi: = 0;
para i: = 8 a swidth - 9 do
para j: = 8 a sheight - 9 do
Si wbpoint [i, j] entonces comience
Características [fi] .info.x: = i;
Características [fi] .info.y: = j;
Características [fi] .dex: = 0;
Inc (fi);
fin;
{输出最终的点序列}
fin;
procedimiento topticalflowlk.init (swidth, sheight, sl: longint);
comenzar
Ancho: = swidth; Altura: = Sheight; L: = sl;
setLength (imageold, ancho, altura, l);
setLength (imageNew, ancho, altura, l);
Marco: = tbitmap.create;
Frame.width: = ancho; Frame.Height: = altura;
Frame.pixEftormat: = Pf24bit;
setLength (ImageGray, ancho, altura);
setLength (valores propios, ancho, altura);
setLength (dx, ancho, altura);
setLength (dy, ancho, altura);
setLength (dxy, ancho, altura);
setLength (wbpoint, ancho, altura);
FeatRecount: = 0;
fin;
Procedimiento Topticalflowlk.makePyramid (Var Images: TtriplelongInTarray; Swidth, Sheight, SL: Longint);
varilla
I, J, K, II, JJ, Nwidth, Nheight, Owidth, Oheight: Longint;
comenzar
{生成金字塔图像}
Owidth: = swidth; OHEIGHT: = Sheight;
para k: = 1 a sl - 1 comience
nwidth: = (owidth + 1) shr 1; nheight: = (oheight + 1) shR 1;
para i: = 1 a nwidth - 2 comience
II: = I shl 1;
para j: = 1 a nheight - 2 comience
JJ: = J SHL 1;
Imágenes [I, J, K]: = (Imágenes [II, JJ, K - 1] SHL 2 +
Imágenes [II - 1, JJ, K - 1] SHL 1 + Imágenes [II + 1, JJ, K - 1] SHL 1 + Imágenes [II, JJ - 1, K - 1] SHL 1 + Imágenes [II, JJ + 1, k - 1] shl 1 +
Imágenes [II - 1, JJ - 1, K - 1] + Imágenes [II + 1, JJ - 1, K - 1] + Imágenes [II - 1, JJ + 1, K - 1] + Imágenes [II + 1 , jj + 1, k - 1]) shr 4;
{高斯原则 ， shl 右移位 ， shr 左移位}
fin;
fin;
para i: = 1 a nwidth - 2 comience
II: = I shl 1;
Imágenes [i, 0, k]: = (imágenes [ii, 0, k - 1] shl 2 +
Imágenes [ii - 1, 0, k - 1] shl 1 + imágenes [ii + 1, 0, k - 1] shl 1 + imágenes [ii, 0, k - 1] shl 1 + imágenes [ii, 1, k - 1] SHL 1 +
Imágenes [II - 1, 0, K - 1] + Imágenes [II + 1, 0, K - 1] + Imágenes [II - 1, 1, K - 1] + Imágenes [II + 1, 1, K - 1 ]) shr 4;
Imágenes [i, nheight - 1, k]: = (imágenes [ii, oheight - 1, k - 1] shl 2 +
Imágenes [ii - 1, oheight - 1, k - 1] shl 1 + imágenes [ii + 1, oheight - 1, k - 1] shl 1 + imágenes [ii, oheight - 2, k - 1] shl 1 + imágenes [ii, oheight - 1, k - 1] shl 1 +
Imágenes [ii - 1, oheight - 2, k - 1] + imágenes [ii + 1, oheight - 2, k - 1] + imágenes [ii - 1, oheight - 1, k - 1] + imágenes [ii + 1 , oheight - 1, k - 1]) shr 4;
fin;
{处理上下边}
para j: = 1 a nheight - 2 comience
JJ: = J SHL 1;
Imágenes [0, J, K]: = (Imágenes [0, JJ, K - 1] SHL 2 +
Imágenes [0, JJ, K - 1] SHL 1 + Imágenes [1, JJ, K - 1] SHL 1 + Imágenes [0, JJ - 1, K - 1] SHL 1 + Imágenes [0, JJ + 1, K - 1] SHL 1 +
Imágenes [0, JJ - 1, K - 1] + Imágenes [1, JJ - 1, K - 1] + Imágenes [0, JJ + 1, K - 1] + Imágenes [1, JJ + 1, K - 1 ]) shr 4;
Imágenes [NWIDTH - 1, J, K]: = (Imágenes [Owidth - 1, JJ, K - 1] SHL 2 +
Imágenes [Owidth - 2, JJ, K - 1] SHL 1 + Imágenes [Owidth - 1, JJ, K - 1] SHL 1 + Imágenes [Owidth - 1, JJ - 1, K - 1] SHL 1 + Imágenes [OWIDTH - 1, jj + 1, k - 1] shl 1 +
Imágenes [Owidth - 2, JJ - 1, K - 1] + Imágenes [Owidth - 1, JJ - 1, K - 1] + Imágenes [Owidth - 2, JJ + 1, K - 1] + Imágenes [Owidth - 1 , jj + 1, k - 1]) shr 4;
fin;
{处理左右边}
Imágenes [0, 0, k]: = (imágenes [0, 0, k - 1] shl 2 +
Imágenes [0, 0, K - 1] SHL 1 + Imágenes [1, 0, K - 1] SHL 1 + Imágenes [0, 0, K - 1] SHL 1 + Imágenes [0, 1, K - 1] SHL 1 +
Imágenes [0, 0, K - 1] + Imágenes [1, 0, K - 1] + Imágenes [0, 1, K - 1] + Imágenes [1, 1, K - 1]) SHR 4;
{处理左上点}
Imágenes [0, nheight - 1, k]: = (imágenes [0, oheight - 1, k - 1] shl 2 +
Imágenes [0, oheight - 1, k - 1] shl 1 + imágenes [1, oheight - 1, k - 1] shl 1 + imágenes [0, oheight - 2, k - 1] shl 1 + imágenes [0, oheight - 1, k - 1] shl 1 +
Imágenes [0, oheight - 2, k - 1] + imágenes [1, oheight - 2, k - 1] + imágenes [0, oheight - 1, k - 1] + imágenes [1, oheight - 1, k - 1 ]) shr 4;
{处理左下点}
Imágenes [NWIDTH - 1, 0, K]: = (Imágenes [Owidth - 1, 0, K - 1] SHL 2 +
Imágenes [Owidth - 2, oheight - 1, k - 1] shl 1 + imágenes [owidth - 1, oheight - 1, k - 1] shl 1 + imágenes [owidth - 1, oheight - 1, k - 1] shl 1 + Imágenes [Owidth - 1, oheight - 1, k - 1] shl 1 +
Imágenes [Owidth - 2, oheight - 1, k - 1] + imágenes [Owidth - 1, oheight - 1, k - 1] + imágenes [owidth - 2, oheight - 1, k - 1] + imágenes [owidth - 1 , oheight - 1, k - 1]) shr 4;
{处理右上点}
Imágenes [nwidth - 1, nheight - 1, k]: = (imágenes [owidth - 1, oheight - 1, k - 1] shl 2 +
Imágenes [owidth - 2, oheight - 1, k - 1] shl 1 + imágenes [owidth - 1, oheight - 1, k - 1] shl 1 + imágenes [owidth - 1, oheight - 2, k - 1] shl 1 + Imágenes [Owidth - 1, oheight - 1, k - 1] shl 1 +
Imágenes [Owidth - 2, oheight - 2, k - 1] + imágenes [Owidth - 1, oheight - 2, k - 1] + imágenes [owidth - 2, oheight - 1, k - 1] + imágenes [owidth - 1 , oheight - 1, k - 1]) shr 4;
{处理右下点}
fin;
fin;
procedimiento topticalflowlk.initfeatures (marcos: tbitmap);
varilla
I, J: Longint;
Línea: prgbtriple;
comenzar
SetStretchBlTMode (frame.canvas.handle, estiramiento_halftone);
STRINGBLT (frame.canvas.handle, 0, 0, ancho, altura, frames.canvas.handle, 0, 0, frames.width, frames.Height, srcCopy);
para i: = 0 a la altura - 1 Comience
Línea: = frame.scanline [i];
para j: = 0 al ancho - 1 Comenzar
Imageold [j, i, 0]: = (línea^.rgbtblue * 11 + línea^.rgbtgreen * 59 + línea^.rgbtred * 30) div 100;
ImageGray [j, i]: = imageold [j, i, 0];
Inclinación);
fin;
fin;
{初始化金字塔图像第一层 Imageold [x, y, 0]}
MakePyramid (imageold, ancho, altura, l);
{生成金字塔图像}
Cornerdetect (ancho, altura, 0.01);
{进行强角点检测}
fin;
procedimiento Topticalflowlk.calopticalflowlk (marcos: tbitmap);
varilla
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: extendido;
Ix, iy: tDoulExtedDedenArray;
Línea: prgbtriple;
Cambio: booleano;
comenzar
SetStretchBlTMode (frame.canvas.handle, estiramiento_halftone);
STRINGBLT (frame.canvas.handle, 0, 0, ancho, altura, frames.canvas.handle, 0, 0, frames.width, frames.Height, srcCopy);
para i: = 0 a la altura - 1 Comience
Línea: = frame.scanline [i];
para j: = 0 al ancho - 1 Comenzar
ImageNew [J, i, 0]: = (línea^.rgbtblue * 11 + línea^.rgbtgreen * 59 + línea^.rgbtred * 30) div 100;
Inclinación);
fin;
fin;
{初始化金字塔图像第一层 ImageNew [x, y, 0]}
MakePyramid (imageNew, ancho, altura, l);
{生成金字塔图像}
setLength (ix, 15, 15); setLength (iy, 15, 15);
{申请差分图像临时数组}
para m: = 0 a treeRecount - 1 Do Comen
{算法细节见:
Jean-Yves Bouguet "Implementación piramidal de Lucas Kanade Tracker Descripción del algoritmo"}
gx: = 0; Gy: = 0;
para ll: = l - 1 downto a 0 comenzar
px: = características [m] .info.x shr ll;
py: = características [m] .info.y shr ll;
{对应当前金字塔图像的 U 点 ： U [l]: = u/2^l}
nwidth: = ancho shr ll; nheight: = altura shr ll;
A: = 0; B: = 0; C: = 0;
para i: = max (px - 7, 1) a min (px + 7, nwidth - 2) hacer
para j: = max (py - 7, 1) a min (py + 7, nheight - 2) comience
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 entonces comience
para k: = 1 a 10 comience
E: = 0; F: = 0;
para i: = max (px - 7, 1) a min (px + 7, nwidth - 2) hacer
para j: = max (py - 7, 1) a min (py + 7, nheight - 2) comience
fi: = i - px + 7; FJ: = j - py + 7;
kx: = i + gx + dx; Ky: = j + gy + dy;
si kx <0 entonces kx: = 0; si kx> nwidth - 1 entonces kx: = nwidth - 1;
Si Ky <0 entonces Ky: = 0; si ky> nheight - 1 entonces 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;
Características [m] .info.x: = px;
Características [m] .info.y: = py;
Características [m] .vector.x: = gx;
Características [m] .vector.y: = gy;
if (px> ancho - 1) o (px <0) o (py> altura - 1) o (py <0) luego presenta [m] .index: = 1;
{失去特征点处理}
fin;
para k: = 0 a l - 1 comience
nwidth: = ancho shr k; nheight: = altura shr k;
para i: = 0 a nwidth - 1 do
para j: = 0 a nheight - 1
Imageold [i, j, k]: = imageNew [i, j, k];
fin;
{复制 J 到 I}
repetir
Cambio: = falso;
para i: = 0 a treeRecount - 1 Do Comen
if características [i] .index = 1 entonces
para j: = i + 1 a treeRecount - 1 Do Comen
Características [j - 1]: = características [j];
Cambio: = verdadero;
fin;
Si cambia, entonces rompa;
fin;
si cambia entonces dec (FeatRecount);
hasta que no cambie;
SetLength (características, FeatRecount);
{删除失去的特征点}
Velocidad: = 0; Radio: = 0; Radiox: = 0; Radioy: = 0;
Si FeatureCount> 0 entonces comienza
SupportCount: = 0;
para i: = 0 a treeRecount - 1 do
if (características [i] .vector.x <> 0) y (características [i] .vector.y <> 0) Entonces comience
Radiox: = Radiox + características [i] .vector.x;
Radioy: = Características radiotinas + [i] .vector.y;
Velocidad: = velocidad + sqrt (sqr (características [i] .vector.x) + sqr (características [i] .vector.y));
Inc (SupportCount);
fin;
Si SupportCount> 0 entonces comienza
Velocidad: = velocidad / soporteCount; //*0.0352778;
Radio: = arctan2 (radiotía, radiox);
fin;
fin;
{计算平均速度和整体方向}
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 1;
Frame.canvas.pen.color: = cllime;
Frame.canvas.brush.style: = bsclear;
para i: = 0 a treeRecount - 1 do
Frame.canvas.ellipse (características [i] .info.x - 4, características [i] .info.y - 4, características [i] .info.x + 4, características [i] .info.y + 4) ;
{用绿圈标识特征点}
Frame.canvas.pen.color: = clyellow;
para i: = 0 a treeRecount - 1 Do Comen
Frame.canvas.moveto (características [i] .info.x - características [i] .vector.x, características [i] .info.y - características [i] .vector.y);
Frame.canvas.lineto (características [i] .info.x, características [i] .info.y);
fin;
{用黄色线条表示运动向量}
if (SupportCount> 0) y (Speed> 3) Entonces comience
Frame.canvas.pen.style: = pssolid;
Frame.canvas.pen.width: = 6;
Frame.canvas.pen.color: = clwhite;
Frame.canvas.ellipse (ancho div 2 - 100, altura div 2 - 100, ancho div 2 + 100, altura div 2 + 100);
Frame.canvas.moveto (ancho div 2, altura div 2);
Frame.canvas.pen.color: = clblue;
Frame.canvas.lineto (ancho div 2 + trunc (90 * cos (radio)), altura div 2 + trunc (90 * sin (radio)));
Frame.canvas.pen.style: = psclear;
Frame.canvas.brush.style: = bssolid;
Frame.canvas.brush.color: = clred;
Frame.canvas.ellipse (ancho div 2 + trunc (90 * cos (radio)) - 8, altura div 2 + trunc (90 * sin (radio)) - 8, ancho div 2 + trunc (90 * cos (radio) ) + 8, altura div 2 + trunc (90 * sin (radio)) + 8);
fin;
fin;
Destructor Topticalflowlk.DestrOY;
comenzar
setLength (imageold, 0);
setLength (imageNew, 0);
setLength (ImageGray, 0);
setLength (EigenValues, 0);
setLength (dx, 0);
setLength (dy, 0);
setLength (dxy, 0);
setLength (wbpoint, 0);
heredado;
fin;
fin.