Lucas-Kanade 算法原理以及应用

Lucas-Kanade 算法原理以及应用

Lucas-Kanade 算法原理以及应用

一 算法原理
1 目标函数

2 一阶泰勒公式展开

3 最小化目标函数条件下的pDelta p

二 LK算在跟踪的应用
1 平移角度尺度版本

2 平移版本

3 平移尺度版本

4 算法流程

三 小结

四 参考文献

一 算法原理

1.1 目标函数

Lucas-Kanade Algorithm本质上是为了最小化目标函数：

∑x[I(W(x;p+Δp))−T(x)]2 （1）
{\sum\limits_x {\left[ {I\left( {W\left( {x;p + \Delta p} \right)} \right) - T\left( x \right)} \right]} ^2}~~~~~~~~（1）

类似高斯牛顿法。其中xx为图像下标，可以是二维（对应图像像素的坐标），也可以是一维（此时为图像展成一维数组时对应的下标）；pp为目标状态变量，W W为仿射变化函数，具体范例如下:

2D平移：

W(x;p)=(x+p1y+p2),p=[p1p2]TW\left( {x;p} \right) = \left( {\begin{array}{*{20}{c}}
{x + {p_1}}\\
{y + {p_2}}
\end{array}} \right),p = {\left[ {\begin{array}{*{20}{c}}
{{p_1}}&{{p_2}}
\end{array}} \right]^T}

3D仿射变化

W(x;p)=(1+p1p2p31+p4p5p6)⎛⎝⎜xy1⎞⎠⎟
W\left( {x;p} \right) = \left( {\begin{array}{*{20}{c}}
{1 + {p_1}}&{{p_3}}&{{p_5}}\\
{{p_2}}&{1 + {p_4}}&{{p_6}}
\end{array}} \right)\left( {\begin{array}{*{20}{c}}
x\\
y\\
1
\end{array}} \right)

其中[p=(p1p2p3p4p5p6)T][p = {\left( {\begin{array}{*{20}{c}}
{{p_1}}&{{p_2}}&{{p_3}}&{{p_4}}&{{p_5}}&{{p_6}}
\end{array}} \right)^T}]

1.2 一阶泰勒公式展开

对公式1进行一阶泰勒公式展开可得

\[∑x[I(W(x;p))+∇I∂W∂pΔp−T(x)]2\] （2）\[{\sum\limits_x {\left[ {I\left( {W\left( {x;p} \right)} \right) + \nabla I\frac{{\partial W}}{{\partial p}}\Delta p - T\left( x \right)} \right]} ^2}\]~~~~~~~（2）

其中

[∇I=(IxIy)][\nabla I = \left( {\begin{array}{*{20}{c}}
{{I_x}}&{{I_y}}
\end{array}} \right)]

设I I已经展开成一列n维向量，则∇I\nabla I为I I在W(x;p)W\left( {x;p} \right)的梯度；

1.3 最小化目标函数条件下的Δp\Delta p

求公式2关于Δp\Delta p的偏导数

2∑x[∇I∂W∂p]T[I(W(x;p))+∇I∂W∂pΔp−T(x)] （3）2\sum\limits_x {{{\left[ {\nabla I\frac{{\partial W}}{{\partial p}}} \right]}^T}\left[ {I\left( {W\left( {x;p} \right)} \right) + \nabla I\frac{{\partial W}}{{\partial p}}\Delta p - T\left( x \right)} \right]} ~~~~~~（3）

让公式3等于0，则

Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))] （4）\Delta p = {H^{{\rm{ - }}1}}\sum\limits_x {{{\left[ {\nabla I\frac{{\partial W}}{{\partial p}}} \right]}^T}\left[ {T\left( x \right) - I\left( {W\left( {x;p} \right)} \right)} \right]}~~~~~~（4）

其中

H=∑x[∇I∂W∂p]T[∇I∂W∂p]H = \sum\limits_x {{{\left[ {\nabla I\frac{{\partial W}}{{\partial p}}} \right]}^T}} \left[ {\nabla I\frac{{\partial W}}{{\partial p}}} \right]

二 LK算在跟踪的应用

这部分将LK算法应用到具体的目标跟踪中，假设跟踪目标用一个角度、尺度可变的矩形进行描述

将矩形框的位移、角度和尺度参数代入公式W(x;p)W\left( {x;p} \right)、xx和pp求得

2.1 平移、角度尺度版本

\[W(x;p)=(xScosθ−ySsinθ+ΔxxSsinθ+yScosθ+Δy)\]\[W\left( {{\bf{x}};} \right.\left. {\bf{p}} \right) = \left( {\begin{array}{*{20}{c}}
{xS\cos \theta - yS\sin \theta + \Delta x}\\
{xS\sin \theta + yS\cos \theta + \Delta y}
\end{array}} \right)\]

变换参数p=(Δx,Δy,θ,S)T{\bf{p}} = {\left( {\Delta x,\Delta y,\theta ,S} \right)^T}，顺时针方向为正方向

则有以下推导

∂W∂p=(1001−xSsinθ−yScosθxScosθ−ySsinθxcosθ−ysinθxsinθ+ycosθ)\frac{{\partial W}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
1&0&{ - xS\sin \theta - yS\cos \theta }&{x\cos \theta - y\sin \theta }\\
0&1&{xS\cos \theta - yS\sin \theta }&{x\sin \theta + y\cos \theta }
\end{array}} \right)

∇I∂W∂p=(∂I∂x∂I∂y)(1001−xSsinθ−yScosθxScosθ−ySsinθxcosθ−ysinθxsinθ+ycosθ)\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\left( {\begin{array}{*{20}{c}}
1&0&{ - xS\sin \theta - yS\cos \theta }&{x\cos \theta - y\sin \theta }\\
0&1&{xS\cos \theta - yS\sin \theta }&{x\sin \theta + y\cos \theta }
\end{array}} \right)

=∂I∂x∂I∂y(−(xsinθ+ycosθ)∂I∂x+(xcosθ−ysinθ)∂I∂y)S(xcosθ−ysinθ)∂I∂x+(xsinθ+ycosθ)∂I∂y=\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}&{\left( { - \left( {x\sin \theta + y\cos \theta } \right)\frac{{\partial I}}{{\partial x}} + \left( {x\cos \theta - y\sin \theta } \right)\frac{{\partial I}}{{\partial y}}} \right)S}&{\left( {x\cos \theta - y\sin \theta } \right)\frac{{\partial I}}{{\partial x}} + \left( {x\sin \theta + y\cos \theta } \right)\frac{{\partial I}}{{\partial y}}}
\end{array}

所以∇I∂W(x;0)∂p=(∂I∂x∂I∂y−y∂I∂x+x∂I∂yx∂I∂x+y∂I∂y)\nabla I\frac{{\partial W\left( {{\bf{x}};} \right.\left. 0 \right)}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}&{ - y\frac{{\partial I}}{{\partial x}} + x\frac{{\partial I}}{{\partial y}}}&{x\frac{{\partial I}}{{\partial x}} + y\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)

2.2 平移版本

[W(x;p)=(x+Δx y+Δy)][W\left( {{\bf{x}};} \right.\left. {\bf{p}} \right) = \left( {\begin{array}{*{20}{c}}
{x + \Delta x}\
{y + \Delta y}
\end{array}} \right)]，其中p=(Δx,Δy)T{\bf{p}} = {\left( {\Delta x,\Delta y} \right)^T}则有：

∂W∂p=(1001)\frac{{\partial W}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
1&0\\
0&1
\end{array}} \right)

\[∇I∂W∂p=(∂I∂x∂I∂y)(1001)=(∂I∂x∂I∂y)\]\[\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\left( {\begin{array}{*{20}{c}}
1&0\\
0&1
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\]

所以可得：[∇I∂W(x;0)∂p=(∂I∂x∂I∂y)][\nabla I\frac{{\partial W\left( {{\bf{x}};} \right.\left. 0 \right)}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)]

2.3 平移、尺度版本

[W(x;p)=(Sx+Δx Sy+Δy)][W\left( {{\bf{x}};} \right.\left. {\bf{p}} \right) = \left( {\begin{array}{*{20}{c}}
{Sx + \Delta x}\
{Sy + \Delta y}
\end{array}} \right)]，其中[p=(Δx,Δy,S)T][{\bf{p}} = {\left( {\Delta x,\Delta y,S} \right)^T}]，可得

\[∇I∂W∂p=(∂I∂x∂I∂y)(1001xy)=(∂I∂x∂I∂y)\]\[\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\left( {\begin{array}{*{20}{c}}
1&0&x\\
0&1&y
\end{array}} \right) = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\]

所以\[∇I∂W(x;0)∂p=(∂I∂x∂I∂yx∂I∂x+y∂I∂y)\]\[\nabla I\frac{{\partial W\left( {{\bf{x}};} \right.\left. 0 \right)}}{{\partial {\bf{p}}}} = \left( {\begin{array}{*{20}{c}}
{\frac{{\partial I}}{{\partial x}}}&{\frac{{\partial I}}{{\partial y}}}&{x\frac{{\partial I}}{{\partial x}} + y\frac{{\partial I}}{{\partial y}}}
\end{array}} \right)\]

2.4 算法流程

根据pp按WW截取图像 II，并对II做归一化；

生成模板的下标矩阵xx，yy；

计算模板的梯度∇I\nabla I

计算∇I∂W(x;0)∂p\nabla I\frac{{\partial W\left( {{\bf{x}};} \right.\left. 0 \right)}}{{\partial {\bf{p}}}}

计算H=∑x[∇I∂W∂p]T[∇I∂W∂p]H = \sum\limits_{\bf{x}} {{{\left[ {\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}}} \right]}^T}\left[ {\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}}} \right]}

计算误差[T(x)−I(W(x;p))]\left[ {T\left( {\bf{x}} \right) - I\left( {W\left( {{\bf{x}};{\bf{p}}} \right)} \right)} \right]

计算\[∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]\[\sum\limits_{\bf{x}} {{{\left[ {\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}}} \right]}^T}\left[ {T\left( {\bf{x}} \right) - I\left( {W\left( {{\bf{x}};{\bf{p}}} \right)} \right)} \right]} \]

根据公式\[Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]\[\Delta {\bf{p}} = {H^{{\rm{ - }}1}}\sum\limits_{\bf{x}} {{{\left[ {\nabla I\frac{{\partial W}}{{\partial {\bf{p}}}}} \right]}^T}\left[ {T\left( {\bf{x}} \right) - I\left( {W\left( {{\bf{x}};{\bf{p}}} \right)} \right)} \right]} \]，求出状态参数的变化

更新目标状态：\[p=p+Δp\]\[{\bf{p}} = {\bf{p}} + \Delta {\bf{p}}\]

判断结果是否收敛，若不收敛，则返回步骤1

三 小结

1、本文的方法本质上是通过梯度下降方法来寻找局部最优解，因此需要初始位置要在最优解的领域内。也就是说，前后两帧目标状态不发发生明显变化的情况。

2、克服目标的大范围运动，可以通过图像金字塔的方法进行跟踪

3、在opencv中，其LK光流算法是实现图像子块的位置跟踪，子块大小一般为5*5.opencv的这个函数是结合图像金字塔，实现对子块的大范围跟踪。但这个函数不能直接得到子块的尺度和角度变化。

4、opencv的LK目标跟踪算法不能应用于光照突变的情况。但是如果目标函数的I和T如果是经过标准化的，相信能提高对光照变化的抗干扰能力

四 参考文献：

Lucas-Kanade 20 Years On: A Unifying Framework IJCV 2004

基于Lucas-kanade目标跟踪算法（本文算法实现代码）

基于光流法的目标跟踪（代码）：使用opencv的稀疏光流法实现的跟踪算法，是一个基于点跟踪的目标跟踪算法

基于前向后向光流的目标跟踪（代码）（Forward-Backward Error: Automatic Detection of Tracking Failures）：基于光流法的目标跟踪的改进算法