{
作者：劉留
參考文獻為：
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）;
遺傳;
結尾;
結尾。