{
作者：刘留
参考文献为：
Jean-Yves Bouguet“卢卡斯·卡纳德（Lucas Kanade）特征跟踪器算法描述的锥体实现”
http://www.aivisoft.net/
geo.cra [at] gmail [dot] com
}
单元光量；
界面
用途
数学，Windows，Sysutils，变体，类，图形，Unitypes，ColorConv;
类型
UpticalFlowlk =类
私人的
Imageold，ImageNw：ttriplelongintarray;
ImageGray，dx，dy，dxy：tdoublelongintarray;
特征值：tdoubleendendedArray;
wbpoint：tdoublebooleanarray;
高度，宽度，L，Roadix，Radioy：longint;
程序CornerDetect（Swidth，Sheight：Longint；质量：扩展）；
过程制作pypyramid（var图像：ttriplelongintarray; swidth，sheight，sl：longint）;
民众
框架：tbitmap;
功能：tsinglepointInfoArray；
featurecount，supportcount：longint;
速度，收音机：扩展；
程序init（swidth，sheight，sl：longint）;
程序initfeatures（帧：tbitMap）;
过程calopticalflowlk（帧：tbitMap）;
毁灭者销毁；覆盖；
结尾;
执行
过程topticalflowlk.cornerdetect（swidth，sheight：longint；质量：扩展）;
var
I，J，FI，FJ：Longint;
a，b，c，sum，minaccept，maxeigenvalue：扩展；
开始
featurecount：= 0;
{
下面采用跟踪的好功能介绍的方法
J. Shi和C. Tomasi“追踪的好功能”，CVPR 94
}
对于我：= 1到旋转-2做
j：= 1 to sheight -2开始
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]）;
结尾;
{求取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;
对于我：= 2到旋转-3做
j：= 2 to sheight -3确实开始
答：= 0; B：= 0; C：= 0;
对于fi：= i -1 to i + 1做
对于fj：= j -1到j + 1确实开始
a：= a + sqr（dx [fi，fj]）;
B：= B + DXY [FI，FJ];
c：= c + sqr（dy [fi，fj]）;
结尾;
答：= a / 2; C：= C / 2;
eigenvalues [i，j]：=（a + c -c -sqrt（（a -c） *（a -c） + b * b））;
如果eigenvalues [i，j]> maxeigenvalue，则maxeigenvalue：= eigenvalues [i，j];
结尾;
{求取矩阵
| ∑DX*dx ∑dxy |
M = | |
| ∑Dxy ∑dy*dy |
的特征值
λ= ∑Dx*dx + ∑dy*dy-（（∑dx*dx + ∑dy*dy）^2-4*（∑dx*dx*dx*∑dy*∑dy*dy -t dy -∑dxy*∑dxy）^1/2）^1/2 }
minaccept：= maxeigenvalue *质量；
{设置最小允许阀值}
对于我：= 8到旋转-9做
对于J：= 8到Sheight -9做
如果特征值[i，j]> minaccept，然后开始
wbpoint [i，j]：= true;
Inc（featurecount）；
结束
wbpoint [i，j]：= false;
对于我：= 8到旋转-9做
对于J：= 8到Sheight -9做
如果wbpoint [i，j]，则开始
sum：=特征值[i，j];
对于fi：= i -8到i + 8确实开始
对于FJ：= J -8到J + 8
如果sqr（fi -i） + sqr（fj -j）<= 64，则
if（eigenvalues [fi，fj]> = sum）和（（fi <> i）或（fj <> j））和（wbpoint [fi，fj]）然后开始
wbpoint [i，j]：= false;
DEC（featurecount）；
休息;
结尾;
如果不是wbpoint [i，j]，则断开；
结尾;
结尾;
{用非最大化抑制来抑制假角点}
setLength（功能，featurecount）; fi：= 0;
对于我：= 8到旋转-9做
对于J：= 8到Sheight -9做
如果wbpoint [i，j]，则开始
功能[fi] .info.x：= i;
功能[fi] .info.y：= j;
功能[fi] .index：= 0;
Inc（FI）;
结尾;
{输出最终的点序列}
结尾;
过程topticalflowlk.init（swidth，sheight，sl：longint）;
开始
宽度：= swidth;身高：= Sheight; l：= sl;
setLength（imageold，宽度，高度，l）;
setLength（Imagenew，宽度，高度，L）;
帧：= tbitmap.create;
frame.width：= width; frame.height：=高度;
frame.pixelformat：= pf24bit;
setLength（imageGray，宽度，高度）;
setLength（特征值，宽度，高度）;
setLength（DX，宽度，高度）;
setLength（dy，width，高度）;
setLength（dxy，宽度，高度）;
setLength（wbpoint，宽度，高度）;
featurecount：= 0;
结尾;
过程topticalflowlk.makepyramid（var图像：ttriplelongintarray; swidth，sheight，sl：longint）;
var
i，j，k，ii，jj，nwidth，nheight，owidth，oheight：longint;
开始
{生成金字塔图像}
owitth：= swidth; Oheight：= Sheight;
对于K：= 1到SL -1确实开始
nwidth：=（owidth + 1）shr 1; nheight：=（oheight + 1）shr 1;
对于我：= 1至nwidth -2确实开始
II：= I SHL 1;
对于J：= 1到Nheight -2确实开始
JJ：= J SHL 1;
图像[I，J，K]：=（图像[II，JJ，K -1] SHL 2 +
图像[II -1，JJ，K -1] SHL 1 +图像[II + 1，JJ，K -1] SHL 1 +图像[II，JJ -1，K -1] SHL 1 +图像[II，JJ + 1，k -1] SHL 1 +
图像[II -1，JJ -1，K -1] +图像[II + 1，JJ -1，K -1] +图像[II -1，JJ + 1，K -1] +图像[II + 1 ，JJ + 1，K -1]）SHR 4;
{高斯原则，shl右移位，shr左移位}
结尾;
结尾;
对于我：= 1至nwidth -2确实开始
II：= I SHL 1;
图像[I，0，K]：=（图像[II，0，K -1] SHL 2 +
图像[II -1，0，K -1] SHL 1 +图像[II + 1，0，K -1] SHL 1 +图像[II，0，K -1] SHL 1 +图像[II，1，K -1] SHL 1 +
图像[II -1，0，K -1] +图像[II + 1，0，K -1] +图像[II -1，1，1，K -1] +图像[II + 1，1，K -1，K -1 ]）SHR 4;
图像[i，nheight -1，k]：=（图像[ii，oheight -1，k -1] shl 2 +
图像[ii -1，oheight -1，k -1] shl 1 +图像[ii + 1，oheight -1，k -1] shl 1 +图像[ii，oheight -2，k -1] shl 1 +图像[ii，Oheight -1，K -1] SHL 1 +
图像[ii -1，Oheight -2，K -1] +图像[II + 1，OHEIGHT -2，K -1] +图像[II -1，OHEIGHT -1，K -1] +图像[II + 1 ，Oheight -1，K -1]）SHR 4;
结尾;
{处理上下边}
对于j：= 1到nheight -2次开始
JJ：= J SHL 1;
图像[0，j，k]：=（图像[0，JJ，K -1] SHL 2 +
图像[0，JJ，K -1] SHL 1 +图像[1，JJ，K -1] SHL 1 +图像[0，JJ -1，K -1] SHL 1 +图像[0，JJ + 1，K -1] SHL 1 +
图像[0，JJ -1，K -1] +图像[1，JJ -1，K -1] +图像[0，JJ + 1，K -1] +图像[1，JJ + 1，K -1，K -1 ]）SHR 4;
图像[nwidth -1，j，k]：=（图像[OWIDTH -1，JJ，K -1] SHL 2 +
图像[OWIDTH -2，JJ，K -1] SHL 1 +图像[OWIDTH -1，JJ，K -1] SHL 1 +图像[OWIDTH -1，JJ -1，K -1，K -1] SHL 1 +图像-1，JJ + 1，K -1] SHL 1 +
图像[OWIDTH -2，JJ -1，K -1] +图像[OWIDTH -1，JJ -1，K -1] +图像[OWIDTH -2，JJ + 1，K -1] +图像[OWIDTH -1 ，JJ + 1，K -1]）SHR 4;
结尾;
{处理左右边}
图像[0，0，K]：=（图像[0，0，K -1] SHL 2 +
图像[0，0，K -1] SHL 1 +图像[1，0，K -1] SHL 1 +图像[0，0，K -1] SHL 1 +图像[0，1，K -1] SHL 1 +
图像[0，0，K -1] +图像[1，0，K -1] +图像[0，1，K -1] +图像[1，1，K -1]）SHR 4;
{处理左上点}
图像[0，Nheight -1，K]：=（图像[0，OHEIGHT -1，K -1] SHL 2 +
图像[0，OHEIGHT -1，K -1] SHL 1 +图像[1，OHEIGHT -1，K -1] SHL 1 +图像[0，OHEIGHT -2，K -1] SHL 1 +图像[0，OHEIGHT -1，k -1] SHL 1 +
图像[0，OHEIGHT -2，K -1] +图像[1，OHEIGHT -2，K -1] +图像[0，OHEIGHT -1，K -1] +图像[1，OHEIGHT -1，K -1，K -1 ]）SHR 4;
{处理左下点}
图像[NWIDTH -1，0，K]：=（图像[OwIdth -1，0，K -1] SHL 2 +
图像[OwIdth -2，Oheight -1，K -1] SHL 1 +图像[OwIdth -1，Oheight -1，K -1] Shl 1 +图像[OwIdth -1，Oheight -1，K -1，K -1] SHL 1] SHL 1 +图像[OWIDTH -1，OHEIGHT -1，K -1] SHL 1 +
图像[OWIDTH -2，OHEIGHT -1，K -1] +图像[OWIDTH -1，OHEIGHT -1，K -1] +图像[OWIDTH -2，OHEIGHT -1，K -1，K -1] +图像[OwIdth -1 ，Oheight -1，K -1]）SHR 4;
{处理右上点}
图像[nwidth -1，nheight -1，k]：=（图像[owidth -1，oheight -1，k -1] shl 2 +
图像[OwIdth -2，Oheight -1，K -1] SHL 1 +图像[OwIdth -1，Oheight -1，K -1] SHL 1 +图像[OwIdth -1，Oheight -2，K -1] Shl 1] Shl 1 +图像[OWIDTH -1，OHEIGHT -1，K -1] SHL 1 +
图像[OWIDTH -2，OHEIGHT -2，K -1] +图像[OWIDTH -1，OHEIGHT -2，K -1] +图像[OWIDTH -2，OHEIGHT -1，K -1，K -1] +图像[OwIdth -1 ，Oheight -1，K -1]）SHR 4;
{处理右下点}
结尾;
结尾;
过程UpticalFlowlk.InitFeatures（帧：TBITMAP）;
var
我，J：longint;
线：prgbtriple;
开始
setStretchbltmode（frame.canvas.handle，strect_halftone）;
StretchBlt（Frame.Canvas.Handle，0，0，0，宽度，高度，框架。
对于I：= 0到高度 - 1个开始
线：= frame.scanline [i];
对于J：= 0到宽度-1个开始
imageold [J，i，0]：=（line^.rgbtblue * 11 + line^.rgbtgreen * 59 + line^.rgbtred * 30）div 100;
ImageGray [J，I]：= ImageOLD [J，i，0];
倾斜）;
结尾;
结尾;
{初始化金字塔图像第一层imageold [x，y，0]}
makepyramid（imageold，宽度，高度，l）;
{生成金字塔图像}
CornerDetect（宽度，高度，0.01）;
{进行强角点检测}
结尾;
过程UpticalFlowlk.caloptictallk（帧：TBITMAP）;
var
i，j，fi，fj，k，ll，m，dx，dy，gx，gy，px，py，kx，ky，ed，ed，edc，nwidth，nheight，nheight：longint;
NX，NY，VX，VY，A，B，C，C，D，E，F，IK：扩展；
ix，iy：tdoublestendendendaray;
线：prgbtriple;
更改：布尔人；
开始
setStretchbltmode（frame.canvas.handle，strect_halftone）;
StretchBlt（Frame.Canvas.Handle，0，0，0，宽度，高度，框架。
对于I：= 0到高度 - 1个开始
线：= frame.scanline [i];
对于J：= 0到宽度-1个开始
Imagenew [j，i，0]：=（line^.rgbtblue * 11 + line^.rgbtgreen * 59 + line^.rgbtred * 30）div 100;
倾斜）;
结尾;
结尾;
{初始化金字塔图像第一层imagenew [x，y，0]}
makepyramid（成像，宽度，高度，L）;
{生成金字塔图像}
setLength（ix，15，15）; setLength（IY，15，15）;
{申请差分图像临时数组}
对于m：= 0到特征 - 1个开始
{算法细节见：
Jean-Yves Bouguet“ Lucas Kanade特征跟踪器算法的Tracker描述”}
GX：= 0; gy：= 0;
对于ll：= l -1降至0的开始
px：=特征[m] .info.x shr ll;
py：=特征[m] .info.y shr ll;
{对应当前金字塔图像的u点：u [l]：= u/2^l}
nwidth：= width shr ll; nheight：=高度shr ll;
答：= 0; B：= 0; C：= 0;
对于i：= max（px -7，1）至min（px + 7，nwidth -2）
对于J：= Max（Py -7，1）至最小（Py + 7，Nheight -2）确实开始
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];
结尾;
{计算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;
如果ABS（D）> 1E-8，则开始
对于K：= 1至10的开始
E：= 0; F：= 0;
对于i：= max（px -7，1）至min（px + 7，nwidth -2）
对于J：= Max（Py -7，1）至最小（Py + 7，Nheight -2）确实开始
fi：= i -px + 7; FJ：= J -PY + 7;
kx：= i + gx + dx; ky：= j + gy + dy;
如果kx <0，则kx：= 0;如果kx> nwidth -1，则kx：= nwidth -1;
如果ky <0，则ky：= 0;如果ky> nheight -1，则ky：= nheight -1;
ik：= imageold [i，j，ll] - imagew [kx，ky，ll];
E：= E + ik * ix [fi，fj];
f：= f + ik * iy [fi，fj];
结尾;
{计算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}
结尾;
结尾;
GX：=（GX + DX）SHL 1; gy：=（gy + dy）shl 1;
{得到下一层的预测运动向量g}
结尾;
GX：= GX Div 2; gy：= gy div 2;
px：= px + gx; py：= py + gy;
功能[m] .info.x：= px;
功能[m] .info.y：= py;
功能[m] .vector.x：= gx;
功能[m] .vector.y：= gy;
if（px> width -1）或（px <0）或（py> height -1）或（py <0），然后功能[m] .index：= 1;
{失去特征点处理}
结尾;
对于k：= 0到l -1确实开始
nwidth：=宽度shr k; nheight：=高度shr k;
对于I：= 0到NWIDTH -1做
对于J：= 0到Nheight -1做
imageold [i，j，k]：= ImageNew [i，j，k];
结尾;
{复制j i}
重复
更改：= false;
对于i：= 0 to temuturecount -1次开始
如果功能[i] .index = 1，则
对于J：= i + 1 to farmaturecount -1确实开始
功能[J -1]：=功能[J];
更改：= true;
结尾;
如果改变，然后休息；
结尾;
如果更改，则DEC（featurecount）；
直到不变；
setLength（功能，featurecount）;
{删除失去的特征点}
速度：= 0;广播：= 0; roadyx：= 0; radioy：= 0;
如果特征> 0，然后开始
supportCount：= 0;
对于I：= 0到特征 - 1做
if（功能[i] .vector.x <> 0）和（功能[i] .vector.y <> 0）然后开始
roadax：= radiox +功能[i] .vector.x;
radioy：= radioy +特征[i] .vector.y;
速度：= speed + sqrt（sqr（sqr（功能[i] .vector.x） + sqr（功能[i] .vector.y））;
Inc（支持量）；
结尾;
如果支持> 0，则开始
速度：=速度 /支持量； //* 0.0352778;
无线电：= Arctan2（Radioy，radax）;
结尾;
结尾;
{计算平均速度和整体方向}
frame.canvas.pen.Style：= pssolid;
frame.canvas.pen.width：= 1;
frame.canvas.pen.Color：= cllime;
frame.canvas.brush.Style：= bsclear;
对于I：= 0到特征 - 1做
frame.canvas.ellipse（功能[i] .info.x -4，功能[i] .info.y -4，功能[i] .info.x + 4，功能[i] .info.y + 4） ;
{用绿圈标识特征点}
frame.canvas.pen.Color：= clyellow;
对于i：= 0 to temuturecount -1次开始
frame.canvas.moveto（功能[i] .info.x-功能[i] .vector.x，功能[i] .info.y-功能[i] .vector.y）;
frame.canvas.lineto（功能[i] .info.x，功能[i] .info.y）;
结尾;
{用黄色线条表示运动向量}
if（supportCount> 0）和（速度> 3），然后开始
frame.canvas.pen.Style：= pssolid;
frame.canvas.pen.width：= 6;
frame.canvas.pen.Color：= clwhite;
Frame.Canvas.Ellipse（宽度Div 2-100，高度Div 2-100，宽度DIV 2 + 100，高度DIV 2 + 100）;
Frame.Canvas.Moveto（宽度Div 2，高度Div 2）;
frame.canvas.pen.Color：= clblue;
frame.canvas.lineto（宽度Div 2 + trunc（90 * cos（无线电）），高度DIV 2 + TRUNC（90 * sin（无线电）））;
frame.canvas.pen.Style：= psclear;
frame.canvas.brush.Style：= bssolid;
frame.canvas.brush.color：= clred;
frame.canvas.ellipse（宽度Div 2 + trunc（90 * cos（Radio））-8，高度DIV 2 + trunc（90 * sin（无线电）） - 8，宽度Div 2 + Trunc（90 * COS（无线电）） ） + 8，高度div 2 + trunc（90 * sin（无线电）） + 8）;
结尾;
结尾;
Destructor upticalflowlk.destroy;
开始
setLength（imageold，0）;
setLength（Imagenew，0）;
setLength（imageGray，0）;
setLength（特征值，0）;
setLength（dx，0）;
setLength（dy，0）;
setLength（dxy，0）;
setLength（wbpoint，0）;
遗传;
结尾;
结尾。