{
作者: 刘留
参考文献为
Jean-Yves Bouguet "Implementação piramidal do Lucas Kanade Recursos Rastreador Descrição do algoritmo"
http://www.aivisoft.net/
Geo.cra [at] gmail [dot] com
}
Unidade Opticalflowlk;
interface
usos
Matemática, Windows, Sysutils, Variantes, Classes, Gráficos, Unitypes, ColorConv;
tipo
Topticalflowlk = classe
Privado
Imageold, imagenew: ttriplelongintary;
ImageGray, DX, DY, DXY: TDOUBLELONGINTARRAY;
Valores próprios: TdoubleExtendedArray;
WBPoint: TdoubleBooleanArray;
Altura, largura, L, Radiox, Radioy: Longint;
Procedimento CornerDetect (Swidth, Sheight: Longint; Qualidade: Extended);
Procedimento makepyramid (Var Images: ttriplelongintary; swidth, sheight, sl: longnt);
público
Quadro: tbitmap;
Recursos: TsinglePointInfoarray;
FeatureCount, SupportCount: Longint;
Velocidade, rádio: estendido;
Procedimento init (Swidth, Sheight, SL: Longint);
Procedimento initFeatures (quadros: tbitmap);
procedimento calopticticlowlk (quadros: tbitmap);
destruidor destruir; substituir;
fim;
implementação
procedimento Topticalflowlk.CornerDetect (Swidth, Sheight: Longint; Qualidade: Extended);
var
I, J, Fi, FJ: Longint;
A, B, C, Soma, Minacept, MaxeigenValue: estendido;
começar
FeatureCount: = 0;
{
下面采用 Bom recurso para rastrear 介绍的方法
J. Shi e C. Tomasi "Bons Recursos para Rastrear", CVPR 94
}
para i: = 1 para swidth - 2 fazer
para j: = 1 para sheight - 2 começa
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]);
fim;
{求取 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;
para i: = 2 para swidth - 3
para j: = 2 para Sheight - 3 começa
a: = 0; b: = 0; C: = 0;
para fi: = i - 1 a i + 1 fazer
para fj: = j - 1 a j + 1 começa
a: = a + sqr (dx [fi, fj]);
b: = b + dxy [fi, fj];
c: = c + sqr (dy [fi, fj]);
fim;
a: = a / 2; c: = c / 2;
Autovalores [i, j]: = (a + c - sqrt ((a - c) * (a - c) + b * b));
Se autovalores [i, j]> maxeigenValue, então maxeigenValue: = eigenValues [i, j];
fim;
{求取矩阵
| ∑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 }
MinAcept: = maxeigenValue * Qualidade;
{设置最小允许阀值}
para i: = 8 para swidth - 9 faça
para j: = 8 para Sheight - 9 faça
Se os valores próprios [i, j]> MinAcept, então comece
Wbpoint [i, j]: = true;
Inc (featureCount);
fim mais
Wbpoint [i, j]: = false;
para i: = 8 para swidth - 9 faça
para j: = 8 para Sheight - 9 faça
Se wbpoint [i, j] então comece
soma: = autovalores [i, j];
para fi: = i - 8 a i + 8 começa
para fj: = j - 8 a j + 8 do
Se sqr (fi - i) + sqr (fj - j) <= 64 então
if (autovalores [fi, fj]> = soma) e ((fi <> i) ou (fj <> j)) e (wbpoint [fi, fj]) então comece
Wbpoint [i, j]: = false;
DEC (featureCount);
quebrar;
fim;
se não wbpoint [i, j] então quebre;
fim;
fim;
{用非最大化抑制来抑制假角点}
setLength (recursos, featureCount); fi: = 0;
para i: = 8 para swidth - 9 faça
para j: = 8 para Sheight - 9 faça
Se wbpoint [i, j] então comece
Apresenta [fi] .info.x: = i;
Apresenta [fi] .info.y: = j;
Apresenta [fi] .Index: = 0;
Inc (fi);
fim;
{输出最终的点序列}
fim;
procedimento Topticalflowlk.init (Swidth, Sheight, SL: Longint);
começar
Largura: = swidth; Altura: = Sheight; L: = SL;
SetLength (Imageold, Largura, Altura, L);
setLength (imagenew, largura, altura, l);
Quadro: = tbitmap.create;
Frame.width: = width; Frame.Height: = altura;
Frame.pixelFormat: = pf24bit;
setLength (ImageGray, largura, altura);
SetLength (autovalores, largura, altura);
setLength (dx, largura, altura);
SetLength (DY, largura, altura);
SetLength (dxy, largura, altura);
setLength (wbpoint, largura, altura);
FeatureCount: = 0;
fim;
procedimento Topticalflowlk.makepiramid (Var Images: ttriplelongintary; Swidth, Sheight, SL: Longint);
var
i, j, k, ii, jj, nwidth, nheight, owidth, ohight: longnt;
começar
{生成金字塔图像}
Owidth: = Swidth; OHEUTH: = Sheight;
para k: = 1 a sl - 1 começa
nWidth: = (Owidth + 1) SHR 1; NHEight: = (OHEIXE + 1) SHR 1;
para i: = 1 a Nwidth - 2 começa
ii: = i shl 1;
para J: = 1 a NHeight - 2 começa
jj: = j shl 1;
Imagens [i, j, k]: = (imagens [ii, jj, k - 1] shl 2 +
Imagens [ii - 1, jj, k - 1] shl 1 + imagens [ii + 1, jj, k - 1] shl 1 + imagens [ii, jj - 1, k - 1] shl 1 + imagens [ii, jj + 1, k - 1] shl 1 +
Imagens [II - 1, JJ - 1, K - 1] + Imagens [II + 1, JJ - 1, K - 1] + Imagens [II - 1, JJ + 1, K - 1] + Imagens [II + 1 , jj + 1, k - 1]) shr 4;
{高斯原则 ， shl 右移位 ， shr 左移位}
fim;
fim;
para i: = 1 a Nwidth - 2 começa
ii: = i shl 1;
Imagens [i, 0, k]: = (imagens [ii, 0, k - 1] shl 2 +
Imagens [ii - 1, 0, k - 1] shl 1 + imagens [ii + 1, 0, k - 1] shl 1 + imagens [ii, 0, k - 1] shl 1 + imagens [ii, 1, k - 1] shl 1 +
Imagens [ii - 1, 0, k - 1] + imagens [ii + 1, 0, k - 1] + imagens [ii - 1, 1, k - 1] + imagens [ii + 1, 1, k - 1 ]) shr 4;
Imagens [i, nheight - 1, k]: = (imagens [ii, ohight - 1, k - 1] shl 2 +
Imagens [ii - 1, ohight - 1, k - 1] shl 1 + imagens [ii + 1, ohight - 1, k - 1] shl 1 + imagens [ii, oheight - 2, k - 1] shl 1 + imagens [ii, ohight - 1, k - 1] shl 1 +
Imagens [ii - 1, ohight - 2, k - 1] + imagens [ii + 1, ohight - 2, k - 1] + imagens [ii - 1, ohight - 1, k - 1] + imagens [ii + 1 , ohight - 1, k - 1]) shr 4;
fim;
{处理上下边}
para J: = 1 a NHeight - 2 começa
jj: = j shl 1;
Imagens [0, J, K]: = (imagens [0, JJ, K - 1] Shl 2 +
Imagens [0, jj, k - 1] shl 1 + imagens [1, jj, k - 1] shl 1 + imagens [0, jj - 1, k - 1] shl 1 + imagens [0, jj + 1, k - 1] shl 1 +
Imagens [0, JJ - 1, K - 1] + Imagens [1, JJ - 1, K - 1] + Imagens [0, JJ + 1, K - 1] + Imagens [1, JJ + 1, K - 1 ]) shr 4;
Imagens [NWidth - 1, J, K]: = (imagens [Owidth - 1, JJ, K - 1] Shl 2 +
Imagens [Owidth - 2, JJ, K - 1] shl 1 + imagens [Owidth - 1, JJ, K - 1] shl 1 + imagens [Owidth - 1, JJ - 1, K - 1] Shl 1 + Imagens [Owidthth - 1, jj + 1, k - 1] shl 1 +
Imagens [Owidth - 2, JJ - 1, K - 1] + Imagens [Owidth - 1, JJ - 1, K - 1] + Imagens [Owidth - 2, JJ + 1, K - 1] + Imagens [Owidth - 1 , jj + 1, k - 1]) shr 4;
fim;
{处理左右边}
Imagens [0, 0, k]: = (imagens [0, 0, k - 1] shl 2 +
Imagens [0, 0, k - 1] shl 1 + imagens [1, 0, k - 1] shl 1 + imagens [0, 0, k - 1] shl 1 + imagens [0, 1, k - 1] shl 1 +
Imagens [0, 0, k - 1] + imagens [1, 0, k - 1] + imagens [0, 1, k - 1] + imagens [1, 1, k - 1]) shr 4;
{处理左上点}
Imagens [0, NHEight - 1, k]: = (imagens [0, ohight - 1, k - 1] shl 2 +
Imagens [0, ohight - 1, k - 1] shl 1 + imagens [1, oheight - 1, k - 1] shl 1 + imagens [0, oheight - 2, k - 1] shl 1 + imagens [0, ohight - 1, k - 1] shl 1 +
Imagens [0, ohight - 2, k - 1] + imagens [1, ohight - 2, k - 1] + imagens [0, oheight - 1, k - 1] + imagens [1, oheight - 1, k - 1 ]) shr 4;
{处理左下点}
Imagens [nWidth - 1, 0, k]: = (imagens [Owidth - 1, 0, k - 1] shl 2 +
Imagens [Owidth - 2, ohight - 1, k - 1] shl 1 + imagens [owidth - 1, ohight - 1, k - 1] shl 1 + imagens [owidth - 1, ohight - 1, k - 1] shl 1 + Imagens [Owidth - 1, ohight - 1, k - 1] shl 1 +
Imagens [Owidth - 2, ohight - 1, k - 1] + imagens [Owidth - 1, ohight - 1, k - 1] + imagens [Owidth - 2, ohight - 1, k - 1] + imagens [Owidth - 1 , ohight - 1, k - 1]) shr 4;
{处理右上点}
Imagens [NWidth - 1, NHEight - 1, K]: = (imagens [Owidth - 1, OHEIXO - 1, K - 1] Shl 2 +
Imagens [Owidth - 2, ohight - 1, k - 1] shl 1 + imagens [owidth - 1, ohight - 1, k - 1] shl 1 + imagens [owidth - 1, ohight - 2, k - 1] shl 1 + Imagens [Owidth - 1, ohight - 1, k - 1] shl 1 +
Imagens [Owidth - 2, ohight - 2, k - 1] + imagens [Owidth - 1, ohight - 2, k - 1] + imagens [Owidth - 2, ohight - 1, k - 1] + imagens [Owidth - 1 , ohight - 1, k - 1]) shr 4;
{处理右下点}
fim;
fim;
Procedimento Topticalflowlk.initFeatures (quadros: TBITMAP);
var
Eu, J: Longint;
Linha: prgbtriple;
começar
SetStretchBltMode (Frame.Canvas.Handle, Stretch_HalftOne);
StretchBlT (ardem.Canvas.Handle, 0, 0, largura, altura, Frames.Canvas.Handle, 0, 0, Frames.Width, Frames.Height, srcCopy);
para i: = 0 a altura - 1 Comece
Linha: = frame.scanline [i];
para j: = 0 para largura - 1 Comece
Imageold [j, i, 0]: = (linha^.rgbtBlue * 11 + linha^.rgbtGreen * 59 + linha^.rgbtred * 30) div 100;
ImageGray [j, i]: = imageold [j, i, 0];
Inclinação);
fim;
fim;
{初始化金字塔图像第一层 Imageold [x, y, 0]}
Makepyramid (imagem, largura, altura, l);
{生成金字塔图像}
CornerDetect (largura, altura, 0,01);
{进行强角点检测}
fim;
Procedimento Topticalflowlk.Calopticalflowlk (quadros: Tbitmap);
var
i, j, fi, fj, k, ll, m, dx, dy, gx, gy, px, py, kx, ky, ed, edc, nwidth, nheight: longnt;
nx, ny, vx, vy, a, b, c, d, e, f, ik: estendido;
Ix, iy: tdoubleExtendedArray;
Linha: prgbtriple;
Mudança: booleano;
começar
SetStretchBltMode (Frame.Canvas.Handle, Stretch_HalftOne);
StretchBlT (ardem.Canvas.Handle, 0, 0, largura, altura, Frames.Canvas.Handle, 0, 0, Frames.Width, Frames.Height, srcCopy);
para i: = 0 a altura - 1 Comece
Linha: = frame.scanline [i];
para j: = 0 para largura - 1 Comece
Imagenew [j, i, 0]: = (linha^.rgbtBlue * 11 + linha^.rgbtGreen * 59 + linha^.rgbtred * 30) div 100;
Inclinação);
fim;
fim;
{初始化金字塔图像第一层 imagenew [x, y, 0]}
Makepyramid (imagenew, largura, altura, l);
{生成金字塔图像}
SetLength (ix, 15, 15); SetLength (IY, 15, 15);
{申请差分图像临时数组}
para m: = 0 para apresentar -se - 1 Comece
{算法细节见:
Jean-Yves Bouguet "Implementação piramidal do Recurso Lucas Kanade Descrição do algoritmo"}
gx: = 0; gy: = 0;
para ll: = l - 1 down 0
px: = recursos [m] .info.x shr ll;
py: = recursos [m] .info.y shr ll;
{对应当前金字塔图像的 u 点: u [l]: = u/2^l}
nWidth: = largura shr ll; nheight: = altura shr ll;
A: = 0; B: = 0; C: = 0;
para i: = max (px - 7, 1) a min (px + 7, nwidth - 2)
para j: = max (py - 7, 1) a min (py + 7, nheight - 2) começa
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];
fim;
{计算 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;
Se abs (d)> 1e-8, então comece
para k: = 1 a 10 começa
E: = 0; F: = 0;
para i: = max (px - 7, 1) a min (px + 7, nwidth - 2)
para j: = max (py - 7, 1) a min (py + 7, nheight - 2) começa
fi: = i - px + 7; fj: = j - py + 7;
kx: = i + gx + dx; KY: = j + gy + dy;
se kx <0 então kx: = 0; se kx> nwidth - 1 então kx: = nwidth - 1;
Se ky <0 então ky: = 0; Se ky> nheight - 1 então ky: = nheight - 1;
Ik: = imageold [i, j, ll] - imagenew [kx, ky, ll];
E: = e + ik * ix [fi, fj];
F: = f + ik * iy [fi, fj];
fim;
{计算 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}
fim;
fim;
gx: = (gx + dx) shl 1; gy: = (gy + dy) shl 1;
{得到下一层的预测运动向量 g}
fim;
gx: = gx div 2; gy: = gy div 2;
px: = px + gx; py: = py + gy;
Apresenta [m] .info.x: = px;
Apresenta [m] .info.y: = py;
Apresenta [m] .vector.x: = gx;
Apresenta [m] .Vector.y: = gy;
if (px> largura - 1) ou (px <0) ou (py> altura - 1) ou (py <0) depois apresenta [m] .index: = 1;
{失去特征点处理}
fim;
para k: = 0 a L - 1 Comece
nWidth: = largura shr k; nheight: = altura shr k;
para i: = 0 a Nwidth - 1 do
para j: = 0 a NHEight - 1 do
Imageold [i, j, k]: = imagenew [i, j, k];
fim;
{复制 j 到 i}
repita
Mudança: = false;
para i: = 0 para apresentar -se - 1 começo
Se recursos [i] .Index = 1 então
para j: = i + 1 para apresentar -se - 1 começo
Recursos [j - 1]: = recursos [j];
Mudança: = true;
fim;
se mudar então quebre;
fim;
se mudar então DEC (featureCount);
até que não mude;
setLength (recursos, featureCount);
{删除失去的特征点}
Velocidade: = 0; Rádio: = 0; Radiox: = 0; Radioy: = 0;
Se featureCount> 0, então comece
Suporte de suporte: = 0;
para i: = 0 para apresentar -se - 1 fazer
if (apresenta [i] .Vector.x <> 0) e (recursos [i] .vector.y <> 0) então comece
Radiox: = Radiox + recursos [i] .vector.x;
Radioy: = radioy + apresenta [i] .vector.y;
Velocidade: = velocidade + sqrt (sqr (recursos [i] .vector.x) + sqr (recursos [i] .vector.y));
Inc (SupportCount);
fim;
Se supportCount> 0, então comece
Velocidade: = velocidade / suporte; //*0.0352778;
Rádio: = arctan2 (radioy, radiox);
fim;
fim;
{计算平均速度和整体方向}
Frame.cacanvas.pen.style: = PSSolid;
Frame.cacanvas.pen.width: = 1;
Frame.canvas.pen.color: = cllime;
Frame.canvas.brush.style: = bsclear;
para i: = 0 para apresentar -se - 1 fazer
Frame.Canvas.ellipse (Recursos [i] .info.x - 4, apresenta [i] .info.y - 4, apresenta [i] .info.x + 4, características [i] .info.y + 4) ;
{用绿圈标识特征点}
Frame.canvas.pen.color: = Clyellow;
para i: = 0 para apresentar -se - 1 começo
Frame.Canvas.moveto (recursos [i] .info.x - recursos [i] .vector.x, recursos [i] .info.y - recursos [i] .vetor.y);
Frame.Canvas.Lineto (apresenta [i] .info.x, recursos [i] .info.y);
fim;
{用黄色线条表示运动向量}
if (supportCount> 0) e (velocidade> 3), então comece
Frame.cacanvas.pen.style: = PSSolid;
Frame.cacanvas.pen.width: = 6;
Frame.canvas.pen.color: = clwhite;
Frame.canvas.ellipse (largura div 2 - 100, altura div 2 - 100, largura div 2 + 100, altura div 2 + 100);
Frame.Canvas.moveto (largura div 2, altura div 2);
Frame.canvas.pen.color: = clblue;
Frame.canvas.lineto (largura div 2 + trunc (90 * cos (rádio)), altura div 2 + trunc (90 * sin (rádio)));
Frame.cacanvas.pen.style: = psclear;
Frame.canvas.brush.style: = bssolid;
Frame.Canvas.brush.color: = clred;
Frame.canvas.ellipse (largura div 2 + trunc (90 * cos (rádio)) - 8, altura div 2 + trunc (90 * sin (rádio)) - 8, largura div 2 + trunc (90 * cos (rádio) ) + 8, altura div 2 + trunc (90 * sin (rádio)) + 8);
fim;
fim;
Destructor Topticalflowlk.Destroy;
começar
setLength (imageold, 0);
setLength (imagenew, 0);
SetLength (ImageGray, 0);
setLength (autovalores, 0);
setLength (dx, 0);
SetLength (dy, 0);
setLength (dxy, 0);
SetLength (WBPoint, 0);
herdado;
fim;
fim.