Lucas-Kanade 算法原理以及应用，正向、反向、additive、Compositional光流法

先祭出一片神级总结性的文章：Lucas-Kanade 20 Years On: A Unifying Framework

Lucas-Kanade 算法原理以及应用

一 算法原理

1.1 目标函数

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

∑x[I(W(x;p+Δp))−T(x)]2        （1）∑x[I(W(x;p+Δp))−T(x)]2        （1）

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

2D平移：

W(x;p)=(x+p1y+p2),p=[p1p2]TW(x;p)=(x+p1y+p2),p=[p1p2]T

3D仿射变化

W(x;p)=(1+p1p2p31+p4p5p6)⎛⎝⎜xy1⎞⎠⎟W(x;p)=(1+p1p3p5p21+p4p6)(xy1)

其中

[p=(p1p2p3p4p5p6)T][p=(p1p2p3p4p5p6)T]

1.2 一阶泰勒公式展开

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

\[∑x[I(W(x;p))+∇I∂W∂pΔp−T(x)]2\]       （2）\[∑x[I(W(x;p))+∇I∂W∂pΔp−T(x)]2\]       （2）

其中

[∇I=(IxIy)][∇I=(IxIy)]

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

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

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

2∑x[∇I∂W∂p]T[I(W(x;p))+∇I∂W∂pΔp−T(x)]      （3）2∑x[∇I∂W∂p]T[I(W(x;p))+∇I∂W∂pΔp−T(x)]      （3）

让公式3等于0，则

Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]      （4）Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]      （4）

其中

H=∑x[∇I∂W∂p]T[∇I∂W∂p]H=∑x[∇I∂W∂p]T[∇I∂W∂p]

二 LK算在跟踪的应用

这部分将LK算法应用到具体的目标跟踪中，假设跟踪目标用一个角度、尺度可变的矩形进行描述
将矩形框的位移、角度和尺度参数代入公式W(x;p)W(x;p)、xx和pp求得

2.1 平移、角度尺度版本

\[W(x;p)=(xScosθ−ySsinθ+ΔxxSsinθ+yScosθ+Δy)\]\[W(x;p)=(xScos⁡θ−ySsin⁡θ+ΔxxSsin⁡θ+yScos⁡θ+Δy)\]

变换参数p=(Δx,Δy,θ,S)Tp=(Δx,Δy,θ,S)T，顺时针方向为正方向
则有以下推导

∂W∂p=(1001−xSsinθ−yScosθxScosθ−ySsinθxcosθ−ysinθxsinθ+ycosθ)∂W∂p=(10−xSsin⁡θ−yScos⁡θxcos⁡θ−ysin⁡θ01xScos⁡θ−ySsin⁡θxsin⁡θ+ycos⁡θ)

∇I∂W∂p=(∂I∂x∂I∂y)(1001−xSsinθ−yScosθxScosθ−ySsinθxcosθ−ysinθxsinθ+ycosθ)∇I∂W∂p=(∂I∂x∂I∂y)(10−xSsin⁡θ−yScos⁡θxcos⁡θ−ysin⁡θ01xScos⁡θ−ySsin⁡θxsin⁡θ+ycos⁡θ)
=∂I∂x∂I∂y(−(xsinθ+ycosθ)∂I∂x+(xcosθ−ysinθ)∂I∂y)S(xcosθ−ysinθ)∂I∂x+(xsinθ+ycosθ)∂I∂y=∂I∂x∂I∂y(−(xsin⁡θ+ycos⁡θ)∂I∂x+(xcos⁡θ−ysin⁡θ)∂I∂y)S(xcos⁡θ−ysin⁡θ)∂I∂x+(xsin⁡θ+ycos⁡θ)∂I∂y

所以∇I∂W(x;0)∂p=(∂I∂x∂I∂y−y∂I∂x+x∂I∂yx∂I∂x+y∂I∂y)∇I∂W(x;0)∂p=(∂I∂x∂I∂y−y∂I∂x+x∂I∂yx∂I∂x+y∂I∂y)

2.2 平移版本

[W(x;p)=(x+Δx y+Δy)][W(x;p)=(x+Δx y+Δy)]，其中p=(Δx,Δy)Tp=(Δx,Δy)T则有：

∂W∂p=(1001)∂W∂p=(1001)

\[∇I∂W∂p=(∂I∂x∂I∂y)(1001)=(∂I∂x∂I∂y)\]\[∇I∂W∂p=(∂I∂x∂I∂y)(1001)=(∂I∂x∂I∂y)\]

所以可得：[∇I∂W(x;0)∂p=(∂I∂x∂I∂y)][∇I∂W(x;0)∂p=(∂I∂x∂I∂y)]

2.3 平移、尺度版本

[W(x;p)=(Sx+Δx Sy+Δy)][W(x;p)=(Sx+Δx Sy+Δy)]，其中[p=(Δx,Δy,S)T][p=(Δx,Δy,S)T]，可得

\[∇I∂W∂p=(∂I∂x∂I∂y)(1001xy)=(∂I∂x∂I∂y)\]\[∇I∂W∂p=(∂I∂x∂I∂y)(10x01y)=(∂I∂x∂I∂y)\]

所以

\[∇I∂W(x;0)∂p=(∂I∂x∂I∂yx∂I∂x+y∂I∂y)\]\[∇I∂W(x;0)∂p=(∂I∂x∂I∂yx∂I∂x+y∂I∂y)\]

2.4 算法流程

1. 根据pp按WW截取图像 II，并对II做归一化；
2. 生成模板的下标矩阵xx，yy；
3. 计算模板的梯度∇I∇I
4. 计算∇I∂W(x;0)∂p∇I∂W(x;0)∂p
5. 计算H=∑x[∇I∂W∂p]T[∇I∂W∂p]H=∑x[∇I∂W∂p]T[∇I∂W∂p]
6. 计算误差[T(x)−I(W(x;p))][T(x)−I(W(x;p))]
7. 计算\[∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]\[∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]
8. 根据公式\[Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]\[Δp=H−1∑x[∇I∂W∂p]T[T(x)−I(W(x;p))]\]，求出状态参数的变化
9. 更新目标状态：\[p=p+Δp\]\[p=p+Δ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）：基于光流法的目标跟踪的改进算法

光流法的目标是完成图像点的跟踪, 因此这里假设存在一个输入图像I, 以及要跟踪的点x, 存在另外一个图像块T, 我们的目标是完成图像块T到输入图像I的匹配.

文章针对Lucas-Canade光流法做了一个总结,
文章对lucas-canade (Forward Additive, FA)算法做了简介, 引入了Compositional算法以及Inverse方法. 因此对应组合形成4种方法分别是Forward Additive(FA), Forward Compositional(FC)以及新提出的Inverse Compositional(IC)算法, Inverse Additive(IA)算法. 对这4种方法文章分别从算法的目标, 算法的推导, 算法对wrap矩阵的要求, 算法计算复杂度以及算法与其他方法的等效情况.

除了四种方法的原理, 文章还介绍了几种优化方法如何实现上述光流的计算. 从高斯牛顿开始介绍, 文章介绍了牛顿法的原理以及对牛顿法的改进. 高斯牛顿法是对牛顿法的一种近似. 当牛顿法中的Hessian矩阵使用雅克比近似时, 牛顿法变为高斯牛顿. 当Hessian矩阵使用c*I(对角为c, 其他元素为0)的矩阵, 牛顿法变为梯度下降法.

Lucas-Kanade 光流法简介

直接法使用光度值不变约束，可以表示为

T(x)=I(W(x;p))T(x)=I(W(x;p))

II表示被匹配的图像， TT表示模板图像。 WW表示关于 xx和 pp的一个函数。如果该函数建模一个平移，那么

W=[x+p1y+p2]W=[x+p1y+p2]

如果是一个仿射变换

W=[(1+p1)x+p3y+p5p2x+(1+p4)y+p6]=[1+p1p2p31+p4p5p6]⎡⎣⎢xy1⎤⎦⎥W=[(1+p1)x+p3y+p5p2x+(1+p4)y+p6]=[1+p1p3p5p21+p4p6][xy1]

这里有6个参数表示了一个仿射变换。
为了让上述问题可解。Lucas等引入FAIA（Forward Additional Image Alignment）方法使用单一的运动模型代替独立像素位移差，其中运动模型的额参数依赖于运模型的建立。之后研究者在其上做了拓展分别是FCIA（Forward Composition Image Alignment），ICIA(Inverse Compositional Image Alignment)和IAIA（Inverse Additional Image Alignment）\cite{baker2004lucas}。其中ICIA使用于直接法SVO中块匹配方法，
看到上面一堆名词，有点晕了，可以看到他们都是围绕Image Alignment来的，那么为什么使用LK光流法呢，LK是这四个中的FAIA。那是因为原始的光流法算起来太慢，然后我们吧光流用在了运动估计上，因此可以根据运动模型降低光流法的计算量。

增量方式\更新方式	forward	inverse
additive	FAIA	IAIA
compositional	FCIA	ICIA

前向与后向的对比

前向方法对于输入图像进行参数化(包括仿射变换及放射增量). 后向方法则同时参数输入图像和模板图像, 其中输入图像参数化仿射变换, 模板图像参数参数化仿射增量. 因此后向方法的计算量显著降低. 由于图像灰度值和运动参数非线性, 整个优化过程为非线性的.

参数化过程主要要计算: 图像的梯度, 位置对运动参数导数, 运动参数增量. 前向方法中Hessian是运动参数的函数. 提高效率的主要思想是交换模板图像和输入图像的角色.
后向方法在迭代中Hessian是固定的.

前向方法和后向方法在目标函数上不太一样，一个是把运动向量pp都是跟着I（被匹配图像），但是前向方法中的迭代的微小量ΔpΔp使用I计算的，后巷方法中的ΔpΔp使用T计算的。因此计算雅克比矩阵的时候，一个的微分在ΔpΔp处，而另外一个在0处。所以如果使用雅克比矩阵计算Hessian矩阵，后者计算的结果是固定的。
举例:
FAIA的目标函数（前向方法）

∑x[I(W(x;p+Δp))−T(x)]2∑x[I(W(x;p+Δp))−T(x)]2

ICIA的目标函数(后向方法)

∑x[T(x;Δp)−I(W(x;p))]2∑x[T(x;Δp)−I(W(x;p))]2

如果使用一阶泰勒展开FAIA(前向方法)的目标函数变为

∑x[I(W(x;p))+∇I∂W∂pΔp−T(x)]2∑x[I(W(x;p))+∇I∂W∂pΔp−T(x)]2

ICIA(后向方法)的泰勒展开为

∑x[T(W(x;0))−I(W(x;p))+∇T∂W∂pΔp]2∑x[T(W(x;0))−I(W(x;p))+∇T∂W∂pΔp]2

而雅克比矩阵为

∂W∂p=[x,0,y,0,1,00,x,0,y,0,1]∂W∂p=[x,0,y,0,1,00,x,0,y,0,1]

只和xy有关，图像I的梯度是要在 W(x;p)W(x;p)处计算的，后向方法中图像T的梯度在 W(x;0)W(x;0) 处计算，因此Hessian矩阵不依赖与处计算，因此Hessian矩阵不依赖与 pp。 后向方法中对模板图像参数化, Hessian矩阵只需要计算一次. 因为模板是在迭代过程中(优化 pp)的每一步固定的。因此会减小计算量。而对输入图像参数化, 由于输入图像的位置是运动的函数, 因此运动参数变化后, 梯度需要重新求解.

Compositional 与 Additive对比

通过增量的表示方式来区分方法. 迭代更新运动参数的时候，如果迭代的结果是在原始的值(6个运动参数)上增加一个微小量，那么称之为Additive，如果在仿射矩阵上乘以一个矩阵(增量运动参数形成的增量仿射矩阵)，这方法称之为Compositional。两者理论上是等效的，而且计算量也差不多。

算法的目标

FAIA:

∑x[I(W(x;p+Δp))−T(x)]2∑x[I(W(x;p+Δp))−T(x)]2

FCIA:

∑x[I(W(W(x;Δp);p))−T(x)]2∑x[I(W(W(x;Δp);p))−T(x)]2

ICIA以及,IAIA:

∑x[T(W(x;Δp))−I(W(x;p))]2∑x[T(W(x;Δp))−I(W(x;p))]2

对于warp的要求

FAIA: W(x,p)W(x,p)对于pp可微.
FCIA: warp集合包含identitywarp, warp集合包含在Composition操作上是闭的(semi-group), 其中包括Homograph, 3D rotation等.
ICIA: semi-group, 另外要求增量warp可逆, 其中包括Homograph, 3D rotation等, 但不包括piece wise affine.
IAIA: 适用于2D平移, 2D相似, 2D仿射等.

算法简介

FAIA

目标函数为

∑x[I(W(x;p+Δp))−T(x)]2∑x[I(W(x;p+Δp))−T(x)]2

更新的方式为

p←p+Δpp←p+Δp

步进的计算方法为

Δp=H−1∑x[I(W)(x;p)−T(x)]Δp=H−1∑x[I(W)(x;p)−T(x)]

算法每个步骤中的时间复杂度

伪代码

FCIA

目标函数为

∑x[I(W(W(x;Δp);p))−T(x)]2∑x[I(W(W(x;Δp);p))−T(x)]2

更新方式为

W(x;p)←W(x;p)∘W(x;Δp)W(x;p)←W(x;p)∘W(x;Δp)

步进的计算

Δp=H−1∑x[I(W)(x;p)−T(x)]Δp=H−1∑x[I(W)(x;p)−T(x)]

ICIA

为了避免花费很多时间来计算hessian矩阵,如果该矩阵是恒定的,那么只需要计算一次.然后事实上Hessian矩阵是关于pp的函数,很多研究给出了该矩阵的近似计算方法,然而很难估计近似的效果,有的时候近似不是很完善.提出该方法的出发点是交换图像和模板,
文章给出了前向和反向的方法是等效的,并给出了证明.

对比IAIA发现ICIA的迭代中不需要对图像梯度进行wrap, 另外计算Hessian中同样如此.

IAIA

仍然是交换I和T. 这样可以避免每个迭代中计算梯度图像.

∑x[I(W(x)+∇T(∂W∂x)−1∂W∂pp−T(x)]2∑x[I(W(x)+∇T(∂W∂x)−1∂W∂pp−T(x)]2

update

p←p+Δpp←p+Δp

实际上这种方法能够使用的运动很少, 对于warp的要求很高, 因此不常用. 文中之后给出了the Inverse Additive 和 Compositional Algorithms 方法在 Affine Warps中的等效性.

总结

• 两个前向方法的计算复杂度相似,后向方法几乎相等.后向方法的速度远比前向方法要快.
• 前向additive可以用于任何变形(warp),
• 反向compositional只能用于warps that form groups.
• 反向additive 可以用于simple 2D linear warps such as translations and affine warps.
  如果不进考虑效率的话可以使用两种前向方法.前向compositional的方法中Jacobian是常量,因此有一定的优势.

如果效率很重要的话,那么后向compositional方法是首选,推导很简单,很容易确定.

雅克比矩阵和残差计算的方式有关, 由于 compositional 计算误差的方式会使得雅克比矩阵为常数,通常采用compositional的形式

梯度下降方法的近似

文章介绍了4种方法分别是高斯牛顿, LM, 梯度下降和Hessian矩阵对角近似. 对这些方法文章分别进行了分step性能, iteration性能等的测试.

牛顿法中通过对Hessian使用雅克比矩阵近似可以得到高斯牛顿.

Algorithm	order of the Taylor approximations	Hessian	work better at
The Gauss-Newton Algorithm	first order		approximations to the Hessian
The Newton Algorithm	sendond order
Steepest Descent	-	indentity matrix	furhter away from lcoal local minima
The Levenberg-Marquardt Algorithm	combine diagonal and full Hessian	error gets worse

结论

1. Gauss-Newton 和 Levenberg-Marquardt的收敛性能类似, 但跟另外2种方法稍好.
2. Levenberg-Marquardt实现的效率和Gauss-Newton接近, 并且不好好于高斯牛顿.

算法的选择

• 噪声，如果图像的噪声比较大，那么最好使用反向算法，反之使用前向方法．
• 关于效率的已经讨论过了．

refer

1. project mainpage https://www.ri.cmu.edu/research_project_detail.html?project_id=515&menu_id=261
   2.https://www.cs.cmu.edu/afs/cs/academic/class/15385-s12/www/lec_slides/Baker&Matthews.pdf