{
作者：刘留
参考文献为：
Jean-Yves Bouguet「Lucas Kanadeのピラミッドの実装機能トラッカーの説明の説明」
http://www.aivisoft.net/
geo.cra [at] gmail [dot] com
}
ユニットOpticalFlowlk;
インタフェース
用途
数学、窓、sysutils、バリアント、クラス、グラフィックス、unitypes、colorconv。
タイプ
topticalflowlk = class
プライベート
ImageOld、ImageNew：Ttriplelongintarray;
ImageGray、dx、dy、dxy：tdoublelongintarray;
固有値：tdoubleextedendedArray;
wbpoint：tdoublebooleanarray;
高さ、幅、L、無線、Radioy：Longint;
手順CornerDetect（Swidth、Sheight：Longint; Quality：Extended）;
手順makepyramid（var images：ttriplelongintarray; swidth、sheight、sl：longint）;
公共
フレーム：tbitmap;
機能：Tsinglepointinfoarray;
featurecount、supportcount：longint;
速度、ラジオ：拡張;
手順init（swidth、sheight、sl：longint）;
手順initfeatures（フレーム：tbitmap）;
手順CalopticalFlowlk（frames：tbitmap）;
デストラクタは破壊します。オーバーライド;
終わり;
実装
手順topticalflowlk.cornerdetect（swidth、sheight：longint; quality：extended）;
var
I、J、FI、FJ：Longint;
a、b、c、sum、minaccept、maxeigenvalue：拡張;
始める
featurecount：= 0;
{
cract介绍的方法を追跡するのに適した機能
J. shiとC. tomasi「追跡するのに適した機能」、CVPR 94
}
i：= 1からswidth -2 do
j：= 1から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;
i：= 2からswidth -3 do
j：= 2からsheight -3は始まります
A：= 0; B：= 0; C：= 0;
fi：= i -1 to i + 1 do
FJの場合：= j -1からj + 1から始まります
a：= a + sqr（dx [fi、fj]）;
b：= b + dxy [fi、fj];
c：= c + sqr（dy [fi、fj]）;
終わり;
a：= a / 2; C：= C / 2;
固有値[i、j]：=（a + c -sqrt（（a -c） *（a -c） + b * b））;
固有値[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*∑dy*dy -∑dxy*∑dxy））^1/2 }
Minaccept：= maxeigenvalue * quality;
{设置最小允许阀值}
i：= 8からswidth -9 do
j：= 8 to sheight -9 do
固有値[i、j]> minacceptの場合、開始します
wbpoint [i、j]：= true;
Inc（featurecount）;
ENSEを終了します
wbpoint [i、j]：= false;
i：= 8からswidth -9 do
j：= 8 to sheight -9 do
wbpoint [i、j]の場合、開始します
合計：=固有値[i、j];
fi：= i -8 to i + 8は始まります
FJの場合：= j -8からj + 8
sqr（fi -i） + sqr（fj -j）<= 64の場合then
if（eigenvalues [fi、fj]> = sum）and（（fi <> i）または（fj <> j））および（wbpoint [fi、fj]）を開始します
wbpoint [i、j]：= false;
12月（featurecount）;
壊す;
終わり;
wbpoint [i、j]ではない場合は壊れます。
終わり;
終わり;
{用非最大化抑制来抑制假角点}
setLength（機能、featurecount）; fi：= 0;
i：= 8からswidth -9 do
j：= 8 to sheight -9 do
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、width、height、l）;
setLength（imageNew、width、height、l）;
フレーム：= tbitmap.create;
frame.width：= width; frame.height：= height;
frame.pixelformat：= pf24bit;
setLength（ImageGray、幅、高さ）;
setLength（固有値、幅、高さ）;
setLength（dx、width、height）;
setLength（dy、width、height）;
setLength（dxy、幅、高さ）;
setLength（wbpoint、width、height）;
featurecount：= 0;
終わり;
手順topticalflowlk.makepyramid（var images：ttriplelongintarray; swidth、sheight、sl：longint）;
var
I、J、K、II、JJ、nwidth、nheight、owidth、oheight：longint;
始める
{生成金字塔图像}
owidth：= swidth; oheight：= sheight;
k：= 1からsl -1から始まります
nwidth：=（owidth + 1）shr 1; nheight：=（oheight + 1）shr 1;
i：= 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左移位}
終わり;
終わり;
i：= 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、k -1] +画像[II + 1、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 ]）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] shl 1 +画像[owidth -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 ]）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] 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] +画像[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 +画像[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] +画像[owidth -1 、oheight -1、k -1]）shr 4;
{处理右下点}
終わり;
終わり;
手順topticalflowlk.initfeatures（frames：tbitmap）;
var
I、J：Longint;
ライン：prgbtriple;
始める
setStretchBltMode（frame.canvas.handle、stretch_halftone）;
Stretchblt（frame.canvas.handle、0、0、width、height、frames.canvas.handle、0、0、frames.width、frames.height、srccopy）;
i：= 0から高さ-1は始まります
行：= frame.scanline [i];
j：= 0からwidth -1は始まります
ImageOld [j、i、0]：=（line^.rgbtblue * 11 + line^.rgbtgreen * 59 + line^.rgbtred * 30）div 100;
ImageGray [j、i]：= imageold [j、i、0];
inc（line）;
終わり;
終わり;
{初始化金字塔图像第一层イメージョルド[x、y、0]}
makepyramid（imageold、width、height、l）;
{生成金字塔图像}
cornerdetect（幅、高さ、0.01）;
{进行强角点检测}
終わり;
手順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：拡張;
IX、IY：TDoubleExtEddedArray;
ライン：prgbtriple;
変更：ブール;
始める
setStretchBltMode（frame.canvas.handle、stretch_halftone）;
Stretchblt（frame.canvas.handle、0、0、width、height、frames.canvas.handle、0、0、frames.width、frames.height、srccopy）;
i：= 0から高さ-1は始まります
行：= frame.scanline [i];
j：= 0からwidth -1は始まります
imageenew [j、i、0]：=（line^.rgbtblue * 11 + line^.rgbtgreen * 59 + line^.rgbtred * 30）div 100;
inc（line）;
終わり;
終わり;
{初始化金字塔图像第一层イメージェウ[x、y、0]}
makepyramid（imageenew、width、height、l）;
{生成金字塔图像}
setLength（ix、15、15）; setLength（IY、15、15）;
{申请差分图像临时数组}
m：= 0 to featurecount -1は始まります
{算法细节见：
Jean-Yves Bouguet「ルーカスカナデのピラミッドの実装機能は、アルゴリズムのトラッカーの説明」}}
GX：= 0; Gy：= 0;
ll：= l -1ダウンまで0を開始します
px：= feature [m] .info.x sh sh ll;
py：= feature [m] .info.y sh ll;
{对应当前金字塔图像的u点：u [l]：= u/2^l}
nwidth：= width sh ll; nheight：= height shr ll;
A：= 0; B：= 0; C：= 0;
i：= max（px -7、1）からmin（px + 7、nwidth -2）do
j：= max（py -7、1）からmin（py + 7、nheight -2）が始まります
fi：= i -px + 7; FJ：= J -Py + 7;
ix [fi、fj]：=（imageold [i + 1、j、ll] - yaugsold [i -1、j、ll]） / 2;
iy [fi、fj]：=（imageold [i、j + 1、ll] - yamesold [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）do
j：= max（py -7、1）からmin（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] -imageNew [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;
feature [m] .vector.y：= gy;
if（px> width -1）または（px <0）または（py> height -1）または（py <0）は[m] .index：= 1;
{失去特征点处理}
終わり;
k：= 0からl -1は始まります
nwidth：= width shr k; nheight：= height shr k;
i：= 0からnwidth -1 do
j：= 0からnheight -1 do
ImageOld [i、j、k]：= imageenew [i、j、k];
終わり;
{复制j到i}
繰り返す
変更：= false;
i：= 0 to featurecount -1は始まります
機能[i] .index = 1の場合
j：= i + 1 to featurecount -1 do begin
機能[j -1]：= feature [j];
変更：= true;
終わり;
変更された場合、壊れます。
終わり;
変更された場合は、12月（featurecount）;
変更されないまで;
setLength（機能、featurecount）;
{删除失去的特征点}
速度：= 0;ラジオ：= 0; Radiox：= 0; Radioy：= 0;
featurecount> 0の場合、開始します
supportcount：= 0;
for i：= 0 to featurecount -1 do
if（feation [i] .vector.x <> 0）および（feature [i] .vector.y <> 0）からBEGIN
Radiox：= Radiox + Feature [i] .vector.x;
radioy：= radioy + feature [i] .vector.y;
速度：=速度 + sqrt（sqr（feation [i] .vector.x） + sqr（feature [i] .vector.y））;
Inc（SupportCount）;
終わり;
supportcount> 0の場合、開始します
速度：= speed / supportCount; //*0.0352778;
無線：= arctan2（radioy、radiox）;
終わり;
終わり;
{计算平均速度和整体方向}
frame.canvas.pen.style：= pssolid;
frame.canvas.pen.width：= 1;
frame.canvas.pen.color：= cllime;
frame.canvas.brush.style：= bsclear;
for i：= 0 to featurecount -1 do
frame.canvas.ellipse（feature [i] .info.x -4、機能[i] .info.y -4、feature [i] .info.x + 4、feation [i] .info.y + 4） ;
{用绿圈标识特征点}
frame.canvas.pen.color：= Clyellow;
i：= 0 to featurecount -1は始まります
frame.canvas.moveto（feature [i] .info.x -feature [i] .vector.x、feature [i] .info.y-機能[i] .vector.y）;
frame.canvas.lineto（feature [i] .info.x、feature [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、height div 2-100、width div 2 + 100、height div 2 + 100）;
frame.canvas.moveto（width div 2、height div 2）;
frame.canvas.pen.color：= clblue;
frame.canvas.lineto（width div 2 + trunc（90 * cos（radio））、height div 2 + trunc（90 * sin（radio）））;
frame.canvas.pen.style：= psclear;
frame.canvas.brush.style：= bssolid;
frame.canvas.brush.color：= clred;
frame.canvas.ellipse（width div 2 + trunc（90 * cos（radio））-8、height div 2 + trunc（90 * sin（radio））-8、width div 2 + trunc（90 * cos（無線）） ） + 8、height div 2 + trunc（90 * sin（radio）） + 8）;
終わり;
終わり;
Destructor topticalflowlk.destroy;
始める
setLength（ImageOld、0）;
setLength（ImageNew、0）;
setLength（ImageGray、0）;
setLength（eigenvalues、0）;
setLength（dx、0）;
setLength（dy、0）;
setLength（dxy、0）;
setLength（wbpoint、0）;
継承;
終わり;
終わり。