Lucas-Kanade 算法原理以及应用

Lucas-Kanade 算法原理以及应用

一 算法原理

1.1 目标函数

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

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

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

2D平移：

$$W\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\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 = {\left( {\begin{array}{*{20}{c}} {{p_1}}&{{p_2}}&{{p_3}}&{{p_4}}&{{p_5}}&{{p_6}} \end{array}} \right)^T}]$$

1.2 一阶泰勒公式展开

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

$$\[{\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）$$
其中
$$[\nabla I = \left( {\begin{array}{*{20}{c}} {{I_x}}&{{I_y}} \end{array}} \right)]$$

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

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

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

$$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，则

$$\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 = \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\left( {x;p} \right)$、$x$和$p$求得

2.1 平移、角度尺度版本

$$\[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)\]$$

变换参数${\bf{p}} = {\left( {\Delta x,\Delta y,\theta ,S} \right)^T}$，顺时针方向为正方向
则有以下推导

$$\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)$$

$\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)$
$=\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}$

所以$\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\left( {{\bf{x}};} \right.\left. {\bf{p}} \right) = \left( {\begin{array}{*{20}{c}} {x + \Delta x}\ {y + \Delta y} \end{array}} \right)]$，其中${\bf{p}} = {\left( {\Delta x,\Delta y} \right)^T}$则有：

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

$$\[\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)\]$$

所以可得：$[\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\left( {{\bf{x}};} \right.\left. {\bf{p}} \right) = \left( {\begin{array}{*{20}{c}} {Sx + \Delta x}\ {Sy + \Delta y} \end{array}} \right)]$，其中$[{\bf{p}} = {\left( {\Delta x,\Delta y,S} \right)^T}]$，可得

$$\[\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)\]$$
所以 $$\[\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 算法流程

1. 根据$p$按$W$截取图像 $I$，并对$I$做归一化；
2. 生成模板的下标矩阵$x$，$y$；
3. 计算模板的梯度$\nabla I$
4. 计算$\nabla I\frac{{\partial W\left( {{\bf{x}};} \right.\left. 0 \right)}}{{\partial {\bf{p}}}}$
5. 计算$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]}$
6. 计算误差$\left[ {T\left( {\bf{x}} \right) - I\left( {W\left( {{\bf{x}};{\bf{p}}} \right)} \right)} \right]$
7. 计算$\[\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]} \]$
8. 根据公式$\[\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]} \]$，求出状态参数的变化
9. 更新目标状态：$\[{\bf{p}} = {\bf{p}} + \Delta {\bf{p}}\]$
10. 判断结果是否收敛，若不收敛，则返回步骤1

三 小结

1、本文的方法本质上是通过梯度下降方法来寻找局部最优解，因此需要初始位置要在最优解的领域内。也就是说，前后两帧目标状态不发发生明显变化的情况。
2、克服目标的大范围运动，可以通过图像金字塔的方法进行跟踪
3、在opencv中，其LK光流算法是实现图像子块的位置跟踪，子块大小一般为5*5.opencv的这个函数是结合图像金字塔，实现对子块的大范围跟踪。但这个函数不能直接得到子块的尺度和角度变化。
4、opencv的LK目标跟踪算法不能应用于光照突变的情况。但是如果目标函数的I和T如果是经过标准化的，相信能提高对光照变化的抗干扰能力

四 参考文献：

1. Lucas-Kanade 20 Years On: A Unifying Framework IJCV 2004
2. 基于Lucas-kanade目标跟踪算法（本文算法实现代码）
3. 基于光流法的目标跟踪（代码）：使用opencv的稀疏光流法实现的跟踪算法，是一个基于点跟踪的目标跟踪算法
4. 基于前向后向光流的目标跟踪（代码）（Forward-Backward Error: Automatic Detection of Tracking Failures）：基于光流法的目标跟踪的改进算法