{
作者:刘留
参考文献为:
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);
継承;
終わり;
終わり。