{
作者 : 刘留
参考文献为 :
Jean-yves Bouguet "Pyramiden-Implementierung der Lucas Kanade-Feature-Tracker-Beschreibung des Algorithmus"
http://www.aivisoft.net/
Geo.cra [at] gmail [dot] com
}
Einheit OpticalFlowlk;
Schnittstelle
Verwendung
Mathematik, Fenster, Systeme, Varianten, Klassen, Grafiken, UnityPes, ColorConv;
Typ
TopticalFlowlk = Klasse
Privat
Imageold, ImageNew: TtriplelongintArray;
ImageGray, DX, DY, DXY: tdoublelongintArray;
Eigenwerte: TDoublextendedArray;
Wbpoint: tDoubleBooleanArray;
Höhe, Breite, L, Radiox, Radioy: Longint;
Verfahren CornerDetekte (Swidth, SHEIGHT: LONGINT; Qualität: erweitert);
Prozedur MakePyramid (var bilder: ttriplelongintArray; swidth, shight, sl: longint);
öffentlich
Rahmen: tbitmap;
Merkmale: TsinglePointinFoArray;
FeatureCount, SupportCount: Longint;
Geschwindigkeit, Radio: erweitert;
Verfahren init (swidth, sheight, sl: longint);
Verfahren initiert (Frames: TBitMap);
Verfahren CalopticalFlowlk (Frames: tbitmap);
Zerstörer zerstören; überschreiben;
Ende;
Durchführung
Verfahren topticalFlowlk.cornerDeTect (Swidth, Sheight: Longint; Qualität: erweitert);
var
Ich, J, Fi, FJ: Longint;
A, B, C, Sum, Minacecept, MaxiderenValue: erweitert;
beginnen
FeatureCount: = 0;
{
下面采用 Gute Funktion, um 介绍的方法 zu verfolgen 介绍的方法
J. Shi und C. Tomasi "Good Features to Track", CVPR 94
}
für i: = 1 bis swidth - 2 tun
für J: = 1 zu Shight - 2 Beginnen Sie an
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]);
Ende;
{求取 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;
für i: = 2 bis swidth - 3 tun
Für J: = 2 bis Sheight - 3 Beginnen Sie an
A: = 0; B: = 0; C: = 0;
für fi: = i - 1 bis i + 1 tun
Für FJ: = j - 1 bis j + 1 Beginnen Sie beginnen
A: = a + sqr (dx [fi, fj]);
B: = b + dxy [fi, fj];
C: = C + SQR (DY [fi, fj]);
Ende;
A: = a / 2; C: = C / 2;
Eigenwerte [i, j]: = (a + c - sqrt ((a - c) * (a - c) + b * b));
Wenn Eigenwerte [i, j]> MaxeigenValue dann MaxiderValue: = Eigenwerte [i, j];
Ende;
{求取矩阵
| ∑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 }
Minacececece: = MaxEinnenwert * Qualität;
{设置最小允许阀值}
für i: = 8 bis swidth - 9 tun
für J: = 8 bis SHEIGHT - 9 tun
Wenn Eigenwerte [i, j]> minaccececept dann beginnen
Wbpoint [i, j]: = true;
Inc (featureCount);
Ende sonst
Wbpoint [i, j]: = false;
für i: = 8 bis swidth - 9 tun
für J: = 8 bis SHEIGHT - 9 tun
Wenn Wbpoint [i, j] dann beginnen
Summe: = Eigenwerte [i, j];
für fi: = i - 8 bis i + 8 beginnen beginnen
für FJ: = j - 8 bis j + 8 tun
Wenn sqr (fi - i) + sqr (fj - j) <= 64 dann dann
if (Eigenvalues [fi, fj]> = sum) und ((fi <> i) oder (fj <> j)) und (wbpoint [fi, fj]) dann beginnen
Wbpoint [i, j]: = false;
Dec (featureCount);
brechen;
Ende;
Wenn nicht wbpoint [i, j], dann brechen Sie;
Ende;
Ende;
{用非最大化抑制来抑制假角点}
setLength (Funktionen, featureCount); fi: = 0;
für i: = 8 bis swidth - 9 tun
für J: = 8 bis SHEIGHT - 9 tun
Wenn Wbpoint [i, j] dann beginnen
Funktionen [fi] .info.x: = i;
Merkmale [fi] .info.y: = j;
Merkmale [fi] .Index: = 0;
Inc (Fi);
Ende;
{输出最终的点序列}
Ende;
Verfahren topticalFlowlk.init (Swidth, SHEIGHT, SL: LONGINT);
beginnen
Breite: = Swidth; Höhe: = SHEIGHT; L: = SL;
SetLength (Imageold, Breite, Höhe, L);
SetLength (ImageNew, Breite, Höhe, L);
Rahmen: = tbitmap.create;
Frame.Width: = Breite; Frame.Height: = Höhe;
Frame.Pixelformat: = PF24bit;
SetLength (ImageGray, Breite, Höhe);
SetLength (Eigenwerte, Breite, Höhe);
SetLength (DX, Breite, Höhe);
SetLength (DY, Breite, Höhe);
SetLength (DXY, Breite, Höhe);
SetLength (Wbpoint, Breite, Höhe);
FeatureCount: = 0;
Ende;
procedure topticalFlowlk.makepyramid (var bilder: ttriplelongintArray; swidth, shight, sl: longint);
var
I, J, K, II, JJ, Nwidth, Nheight, Owidth, Ohht: Longint;
beginnen
{生成金字塔图像}
Owidth: = Swidth; Oheight: = Shight;
für k: = 1 bis sl - 1 Beginnen Sie an
NWIDTH: = (OWIDTH + 1) SHR 1; nHeight: = (ohtheight + 1) SHR 1;
für i: = 1 bis nwidth - 2 Beginnen Sie an
II: = I Shl 1;
für J: = 1 bis nHeight - 2 Beginnen Sie beginnen
JJ: = J SHL 1;
Bilder [i, j, k]: = (Bilder [II, JJ, K - 1] Shl 2 +
Bilder [II - 1, JJ, K - 1] SHL 1 + Bilder [II + 1, JJ, K - 1] SHL 1 + Bilder [II, JJ - 1, K - 1] SHL 1 + Bilder [II, JJ + 1, k - 1] SHL 1 +
Bilder [II - 1, JJ - 1, K - 1] + Bilder [II + 1, JJ - 1, K - 1] + Bilder [II - 1, JJ + 1, K - 1] + Bilder [II + 1 , JJ + 1, k - 1]) SHR 4;
{高斯原则 , Shl 右移位 , Shr 左移位}
Ende;
Ende;
für i: = 1 bis nwidth - 2 Beginnen Sie an
II: = I Shl 1;
Bilder [i, 0, k]: = (Bilder [ii, 0, k - 1] Shl 2 +
Bilder [II - 1, 0, K - 1] SHL 1 + Bilder [II + 1, 0, K - 1] SHL 1 + Bilder [II, 0, K - 1] Shl 1 + Bilder [II, 1, K. - 1] Shl 1 +
Bilder [II - 1, 0, K - 1] + Bilder [II + 1, 0, K - 1] + Bilder [II - 1, 1, K - 1] + Bilder [II + 1, 1, K - 1 ]) SHR 4;
Bilder [i, nHeight - 1, k]: = (Bilder [ii, oheight - 1, k - 1] SHL 2 +
Bilder [ii - 1, oheight - 1, k - 1] SHL 1 + Bilder [II + 1, Oheight - 1, K - 1] SHL 1 + Bilder [II, OHTH - 2, K - 1] SHL 1 + Bilder [ii, oheight - 1, k - 1] SHL 1 +
Bilder [ii - 1, oheight - 2, k - 1] + Bilder [ii + 1, oheight - 2, k - 1] + Bilder [ii - 1, oheight - 1, k - 1] + Bilder [II + 1 , Oheight - 1, K - 1]) SHR 4;
Ende;
{处理上下边}
für J: = 1 bis nHeight - 2 Beginnen Sie beginnen
JJ: = J SHL 1;
Bilder [0, j, k]: = (Bilder [0, jj, k - 1] Shl 2 +
Bilder [0, JJ, K - 1] SHL 1 + Bilder [1, JJ, K - 1] SHL 1 + Bilder [0, JJ - 1, K - 1] SHL 1 + Bilder [0, JJ + 1, k - 1] Shl 1 +
Bilder [0, jj - 1, k - 1] + Bilder [1, jj - 1, k - 1] + Bilder [0, jj + 1, k - 1] + Bilder [1, jj + 1, k - 1 ]) SHR 4;
Bilder [NWIDTH - 1, J, K]: = (Bilder [Owidth - 1, JJ, K - 1] Shl 2 +
Bilder [Owidth - 2, JJ, K - 1] SHL 1 + Bilder [Owidth - 1, JJ, K - 1] SHL 1 + Bilder [Owidth - 1, JJ - 1, K - 1] SHL 1 + Bilder [Owidth - 1, jj + 1, k - 1] SHL 1 +
Bilder [Owidth - 2, jj - 1, k - 1] + Bilder [Owidth - 1, jj - 1, k - 1] + Bilder [Owidth - 2, JJ + 1, k - 1] + Bilder [Owidth - 1 , JJ + 1, k - 1]) SHR 4;
Ende;
{处理左右边}
Bilder [0, 0, k]: = (Bilder [0, 0, k - 1] Shl 2 +
Bilder [0, 0, k - 1] Shl 1 + Bilder [1, 0, k - 1] SHL 1 + Bilder [0, 0, k - 1] Shl 1 + Bilder [0, 1, k - 1] Shl 1 +
Bilder [0, 0, k - 1] + Bilder [1, 0, k - 1] + Bilder [0, 1, k - 1] + Bilder [1, 1, k - 1]) SHR 4;
{处理左上点}
Bilder [0, nHeight - 1, k]: = (Bilder [0, Oheight - 1, K - 1] SHL 2 +
Bilder [0, Oheight - 1, K - 1] SHL 1 + Bilder [1, Oheight - 1, K - 1] SHL 1 + Bilder [0, Oheight - 2, K - 1] SHL 1 + Bilder [0, Oheight - 1, k - 1] SHL 1 +
Bilder [0, Oheight - 2, K - 1] + Bilder [1, oheight - 2, k - 1] + Bilder [0, Oheight - 1, K - 1] + Bilder [1, ohtyht - 1, k - 1 ]) SHR 4;
{处理左下点}
Bilder [Nwidth - 1, 0, k]: = (Bilder [Owidth - 1, 0, k - 1] Shl 2 +
Bilder [owidth - 2, oheight - 1, k - 1] SHL 1 + Bilder [Owidth - 1, Oheight - 1, K - 1] SHL 1 + Bilder [Owidth - 1, Oheight - 1, K - 1] SHL 1 + Bilder [owidth - 1, oheight - 1, k - 1] SHL 1 +
Bilder [owidth - 2, oheight - 1, k - 1] + bilder [owidth - 1, oheight - 1, k - 1] + bilder [owidth - 2, oheight - 1, k - 1] + Bilder [Owidth - 1 , Oheight - 1, K - 1]) SHR 4;
{处理右上点}
Bilder [NWIDTH - 1, NHEIGHT - 1, K]: = (Bilder [owidth - 1, oheight - 1, k - 1] Shl 2 +
Bilder [owidth - 2, oheight - 1, k - 1] SHL 1 + Bilder [Owidth - 1, Oheight - 1, K - 1] SHL 1 + Bilder [Owidth - 1, oheight - 2, k - 1] SHL 1 + Bilder [owidth - 1, oheight - 1, k - 1] SHL 1 +
Bilder [owidth - 2, oheight - 2, k - 1] + bilder [owidth - 1, oheight - 2, k - 1] + bilder [owidth - 2, oheight - 1, k - 1] + Bilder [Owidth - 1 , Oheight - 1, K - 1]) SHR 4;
{处理右下点}
Ende;
Ende;
procedure topticalFlowlk.initFeatures (Frames: tbitmap);
var
Ich, J: Longint;
Zeile: prgbtrips;
beginnen
SetStRetchBltMode (Frame.Canvas.Handle, Stretch_Halftone);
StretchBlt (Frame.Canvas.Handle, 0, 0, Breite, Höhe, Frames.Canvas.Handle, 0, 0, Frames.width, Frames.Height, Srccopy);
Für i: = 0 bis zur Höhe - 1 Beginnen Sie an
Zeile: = Frame.Scanline [i];
für J: = 0 bis Breite - 1 Beginnen Sie an
Imageold [j, i, 0]: = (Zeile^.rgBtBlue * 11 + line^.rgbtgreen * 59 + line^.rgbTred * 30) div 100;
ImageGray [J, I]: = Imageold [j, i, 0];
Neigung);
Ende;
Ende;
{初始化金字塔图像第一层 Imageold [x, y, 0]}
Makepyramid (Imageold, Breite, Höhe, L);
{生成金字塔图像}
Cornerdetekte (Breite, Höhe, 0,01);
{进行强角点检测}
Ende;
procedure topticalflowlk.calopticalFlowlk (Frames: tbitmap);
var
I, J, Fi, FJ, K, LL, M, DX, Dy, GX, GY, PX, PY, KX, KY, ED, EDC, NWIDTH, NHELT: Longint;
NX, NY, VX, VY, A, B, C, D, E, F, IK: Extended;
Ix, iy: tDoublextendedArray;
Zeile: prgbtrips;
Änderung: Boolean;
beginnen
SetStRetchBltMode (Frame.Canvas.Handle, Stretch_Halftone);
StretchBlt (Frame.Canvas.Handle, 0, 0, Breite, Höhe, Frames.Canvas.Handle, 0, 0, Frames.width, Frames.Height, Srccopy);
Für i: = 0 bis zur Höhe - 1 Beginnen Sie an
Zeile: = Frame.Scanline [i];
für J: = 0 bis Breite - 1 Beginnen Sie an
ImageNew [j, i, 0]: = (Zeile^.rgBtBlue * 11 + line^.rgbtgreen * 59 + line^.rgbTred * 30) div 100;
Neigung);
Ende;
Ende;
{初始化金字塔图像第一层 ImageNew [x, y, 0]}
Makepyramid (ImageNew, Breite, Höhe, L);
{生成金字塔图像}
SetLength (IX, 15, 15); SetLength (iy, 15, 15);
{申请差分图像临时数组}
für m: = 0 zu featureCount - 1 Beginnen Sie an
{算法细节见:
Jean-yves Bouguet "Pyramiden-Implementierung der Lucas Kanade-Feature-Tracker Beschreibung des Algorithmus"}
GX: = 0; gy: = 0;
Für LL: = l - 1 nach unten beginnen beginnen
px: = features [m] .info.x Shr ll;
PY: = Merkmale [m] .info.y Shr ll;
{对应当前金字塔图像的 u 点 : u [l]: = u/2^l}
NWIDTH: = Breite SHR ll; nHeight: = Höhe SHR ll;
A: = 0; B: = 0; C: = 0;
für i: = max (px - 7, 1) bis min (px + 7, nwidth - 2) tun
für J: = max (py - 7, 1) bis min (py + 7, nHeight - 2) beginnen beginnen
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];
Ende;
{计算 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;
Wenn ABS (d)> 1E-8 dann beginnen
für k: = 1 bis 10 beginnen beginnen
E: = 0; F: = 0;
für i: = max (px - 7, 1) bis min (px + 7, nwidth - 2) tun
für J: = max (py - 7, 1) bis min (py + 7, nHeight - 2) beginnen beginnen
fi: = i - px + 7; FJ: = J - PY + 7;
kx: = i + gx + dx; ky: = j + gy + dy;
Wenn KX <0 dann KX: = 0; Wenn KX> NWIDTH - 1, dann KX: = NWIDTH - 1;
Wenn Ky <0 dann KY: = 0; Wenn KY> nHeight - 1 dann KY: = nHeight - 1;
IK: = Imageold [i, j, ll] - ImageNew [KX, KY, LL];
E: = e + ik * ix [fi, fj];
F: = f + ik * iy [fi, fj];
Ende;
{计算 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}
Ende;
Ende;
gx: = (gx + dx) SHL 1; gy: = (gy + dy) SHL 1;
{得到下一层的预测运动向量 g}
Ende;
gx: = gx div 2; gy: = gy div 2;
px: = px + gx; py: = py + gy;
Merkmale [m] .info.x: = px;
Merkmale [m] .info.y: = py;
Merkmale [m] .Vector.x: = gx;
Merkmale [m] .Vector.y: = gy;
if (px> breit - 1) oder (px <0) oder (py> Höhe - 1) oder (py <0), fasst [m] .Index: = 1;
{失去特征点处理}
Ende;
für k: = 0 bis l - 1 Beginnen Sie an
NWIDTH: = Breite SHR K; nHeight: = Höhe SHR K;
für i: = 0 bis nwidth - 1 tun
für J: = 0 bis nHeight - 1 tun
Imageold [i, j, k]: = imageNew [i, j, k];
Ende;
{复制 j 到 i}
wiederholen
Änderung: = falsch;
für i: = 0 zu featureCount - 1 beginnen
Wenn Funktionen [i] .Index = 1 dann
für J: = i + 1 zu featureCount - 1 Beginnen Sie an
Merkmale [j - 1]: = Features [j];
Änderung: = wahr;
Ende;
Wenn ändert, dann brechen Sie;
Ende;
Wenn Sie sich ändern, dann dec (featureCount);
bis sich nicht ändert;
setLength (Funktionen, featureCount);
{删除失去的特征点}
Geschwindigkeit: = 0; Radio: = 0; Radiox: = 0; Radioy: = 0;
Wenn featureCount> 0 dann beginnen
SupportCount: = 0;
für i: = 0 zu featureCount - 1 tun
if (features [i] .vector.x <> 0) und (features [i] .Vector.y <> 0) dann beginnen
Radiox: = Radiox + Funktionen [i] .Vector.x;
Radioy: = Radioy + Merkmale [i] .Vector.y;
Geschwindigkeit: = Geschwindigkeit + SQRT (SQR (Features [i] .Vector.x) + SQR (Funktionen [i] .Vector.y));
Inc (SupportCount);
Ende;
Wenn SupportCount> 0 dann beginnen
Geschwindigkeit: = Geschwindigkeit / Stützcount; //*0.0352778;
Radio: = Arctan2 (Radioy, Radiox);
Ende;
Ende;
{计算平均速度和整体方向}
Frame.Canvas.pen.Style: = psSolid;
Frame.Canvas.Pen.Width: = 1;
Frame.Canvas.pen.Color: = Cllim;
Frame.Canvas.brush.Style: = bsclear;
für i: = 0 zu featureCount - 1 tun
Frame.Canvas.ellipse (Funktionen [i] .info.x - 4, Funktionen [i] .info.y - 4, Funktionen [i] .info.x + 4, Funktionen [i] .info.y + 4) ;
{用绿圈标识特征点}
Frame.Canvas.pen.Color: = Clyellow;
für i: = 0 zu featureCount - 1 beginnen
Frame.Canvas.Moveto (Funktionen [i] .info.x - Funktionen [i] .Vector.x, Funktionen [i] .info.y - Funktionen [i] .Vector.y);
Frame.Canvas.lineto (Funktionen [i] .info.x, Funktionen [i] .info.y);
Ende;
{用黄色线条表示运动向量}
if (SupportCount> 0) und (Geschwindigkeit> 3) beginnen
Frame.Canvas.pen.Style: = psSolid;
Frame.Canvas.Pen.Width: = 6;
Frame.Canvas.pen.Color: = clwhite;
Frame.Canvas.ellipse (Breite Div 2 - 100, Höhe Div 2 - 100, Breite Div 2 + 100, Höhe Div 2 + 100);
Frame.Canvas.Moveto (Breite Div 2, Höhe Div 2);
Frame.Canvas.pen.Color: = clbblue;
Frame.Canvas.lineto (Breite Div 2 + Trunc (90 * cos (Radio)), Höhe Div 2 + Trunc (90 * sin (Radio)));
Frame.Canvas.Pen.Style: = PSClear;
Frame.Canvas.brush.Style: = bssolid;
Frame.Canvas.brush.Color: = Clred;
Frame.Canvas.ellipse (Breite Div 2 + Trun (90 * cos (Radio)) - 8, Höhe Div 2 + Trun (90 * sin (Radio)) - 8, Breite Div 2 + TRUNC (90 * cos (Radio) ) + 8, Höhe Div 2 + TRUNC (90 * sin (Radio)) + 8);
Ende;
Ende;
destructor topticalflowlk.destroy;
beginnen
setLength (Imageold, 0);
setLength (ImageNew, 0);
setLength (ImageGray, 0);
SetLength (Eigenwerte, 0);
setLength (dx, 0);
setLength (dy, 0);
setLength (dxy, 0);
setLength (wbpoint, 0);
geerbt;
Ende;
Ende.