Lucas Kanade 光流法（来自wiki 百科）

来源：互联网发布：淘宝处方药药师问什么编辑：程序博客网时间：2024/05/20 07:36

小伙伴们开始正式玩起了APM的PX4flow，加上课题方向也要用到光流法，因此从哪个角度来说，都是十分必要的。

光流法最常用的是用于机器视觉的跟踪算法，一是可以跟踪目标物体，而是求解目标的运动学参数（平移速度与旋转速度）。

光流法用到了一个很重要的概念，叫做场。图像像素点亮度的变化，形成了与空域有关的一个场，由于云台上的相机是运动的，它取景的位置是与时间有关，这构成了图像的时域。

本科学过热扩散方程，温度场也是一种时空结合的场，热扩散方程的建立的前提是能量守恒定律。如果所研究的温度场与外界绝缘，那么温度场内的热量是守恒的，这样就建立起了一种热扩散方程；如果所研究的温度场受到其他能量的作用，考虑到能量守恒，就能得到广义热扩散方程。

既然光流法研究的也是时空结合的场，要想建立类似的方程，必须做出类似的假设：研究的特征点领域的像素的亮度值在dt时间段前后是恒定的。

提取图像特征点的方法和标准不再本文讨论方法之内，我们主要看光流法方程如何建立的。

某t时刻p特征点图像附近区域的图像记作 $I(x,y,t)$ ，其中x，y是相对于p特征点的坐标。过了之后，p特征点在视野中的位置发生了改变，x,y方向的位移为 $\Delta x$ , $\Delta y$ ，此时时刻为 $t+\Delta t$ ，根据所做的假设，p特征点的领域图像时保持恒定的，那么有

$I(x,y,t) = I(x+\Delta x, y + \Delta y, t + \Delta t)$

我们倘若在新图像中搜索与原图像亮度值一致的区域，是可以做到的，这就是光流法中的SAD（sum of absolute difference）法，但会遇到一个很可怕的问题，我们不知道这段时间内图像有没有发生旋转，如果发生了旋转，又不知道旋转的方向大小，这个匹配起来太困难了，因为无法确定搜索区域的方向。

当然PX4等飞行器是能解决这个问题的，因为它们配备陀螺仪，能够精确的记录相机坐标系的旋转，矩阵旋转就能得到图像的旋转大小和方向。

图像算法本身能够解决这个问题，这里就要借助数学工具求解出了。

将等式右边用Taylor公式展开，得到

$I(x+\Delta x,y+\Delta y,t+\Delta t) = I(x,y,t) + \frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t+$ H.O.T.

略去高阶无穷小和利用 $I(x,y,t) = I(x+\Delta x, y + \Delta y, t + \Delta t)$ ，有

\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t = 0

即

\frac{\partial I}{\partial x}\frac{\Delta x}{\Delta t}+\frac{\partial I}{\partial y}\frac{\Delta y}{\Delta t}+\frac{\partial I}{\partial t}\frac{\Delta t}{\Delta t} = 0

那么有

\frac{\partial I}{\partial x}V_x+\frac{\partial I}{\partial y}V_y+\frac{\partial I}{\partial t} = 0

这样就能得到相机的平移速度啦！（这个和热扩散方程极像，连思想方法都一样。）

上面提到 $I(x,y,t)$ 为某t时刻p特征点图像附近区域的图像，该图像区域里有众多像素，可将光流法方程写的详细点，得到下式：

I_x(q_1) V_x + I_y (q_1) V_y = -I_t(q_1)

I_x(q_2) V_x + I_y (q_2) V_y = -I_t(q_2)

\vdots

I_x(q_n) V_x + I_y (q_n) V_y = -I_t(q_n)

这里 $q_1,q_2,\dots,q_n$ 是该图像区域内的像素点， $I_x(q_i),I_y(q_i),I_t(q_i)$ 为图像 $I$ 在此时在 $q_i$ 点领域里的多元函数分别对位置 x, y 和时间 t的偏微分估计。

线性代数这是就派上用场了，上面这么多关系式可以统一写成 $A v = b$ , 这里

A = \begin{bmatrix}I_x(q_1) & I_y(q_1) \\[10pt]I_x(q_2) & I_y(q_2) \\[10pt]\vdots & \vdots \\[10pt]I_x(q_n) & I_y(q_n) \end{bmatrix},\quad\quadv = \begin{bmatrix}V_x\\[10pt]V_y\end{bmatrix},\quad \mbox{and}\quadb = \begin{bmatrix}-I_t(q_1) \\[10pt]-I_t(q_2) \\[10pt]\vdots \\[10pt]-I_t(q_n)\end{bmatrix}

这是一个超定方程，因为未知数只有Vx,Vy两个，而方程个数往往远远多于这些，在矩阵分析中，A矩阵m*n维，且不是方阵，求解V时使用的是A矩阵的伪逆，伪逆的形式有很多种，具体就不展开了。Lucas-Kanade法在这里选取了线性最小二乘法的伪逆。即：线性方程组两边同时乘上A的转置，则

A^T A v=A^T b

这样就可以求逆了：

\mathrm{v}=(A^T A)^{-1}A^T b

注意：伪逆不是真正的逆，使用上述伪逆求解之后，很可能上述n个等式都不严格相等，而是近似相等，每个等式都存在偏差。线性最小二乘法的伪逆能保证所有的偏差的平方之和最小，这就是最小二乘法的含义。

题外话：wiki百科中一句话真是醍醐灌顶：使用最小二乘法的依据是我们认为偏差e是服从于正态分布，准确的表达式为：。作为一名工科生，做了这么多年实验，当然实验内容有的实在是不敢恭维，但是基本数据处理还是自己做的多，几乎无一例外地使用了最小二乘法（MATLAB拟合工具箱或者是polyfit函数），却不知道这背后的意义。很重要的原因是，概率、统计、随机这一块数学学得少而且抱着得过且过的态度，自然，还包括在之前的文章里提到我对着方面理解能力有限的原因。

这个题外话恰恰从物理原理上揭示了为什么使用线性最小二乘法的伪逆来求解光流速度，因为光流法遇到的噪声多是随机噪声，服从正态分布。

我们还可以从另一个数学问题诉求物理原理，我们知道Taylor展开时近似求导的时候，要求 $\Delta x$ , $\Delta y$ ，都足够小，这样的近似才更有意义，那么，就要求光流速度足够小，所以Lucas-Kanade法求解光流一个很重要的前提是相邻两帧之间的图像的特征点p的偏移很小，否则求解结果是不准的，这就诞生了许多方法：如金字塔式的LK法，LKT法等等，不再展开

后续（翻阅了刘未鹏和July的博客，终于找到了直线拟合的二乘法解释，参见从决策树学习谈到贝叶斯分类算法、EM、HMM和数学之美番外篇：平凡而又神奇的贝叶斯方法）：

2.4.3、最大似然与最小二乘

学过线性代数的大概都知道经典的最小二乘方法来做线性回归。问题描述是：给定平面上 N 个点，（这里不妨假设我们想用一条直线来拟合这些点——回归可以看作是拟合的特例，即允许误差的拟合），找出一条最佳描述了这些点的直线。

一个接踵而来的问题就是，我们如何定义最佳？我们设每个点的坐标为 (Xi, Yi) 。如果直线为 y = f(x) 。那么 (Xi, Yi) 跟直线对这个点的“预测”：(Xi, f(Xi)) 就相差了一个 ΔYi = |Yi – f(Xi)| 。最小二乘就是说寻找直线使得 (ΔY1)^2 + (ΔY2)^2 + .. （即误差的平方和）最小，至于为什么是误差的平方和而不是误差的绝对值和，统计学上也没有什么好的解释。然而贝叶斯方法却能对此提供一个完美的解释。

我们假设直线对于坐标 Xi 给出的预测 f(Xi) 是最靠谱的预测，所有纵坐标偏离 f(Xi) 的那些数据点都含有噪音，是噪音使得它们偏离了完美的一条直线，一个合理的假设就是偏离路线越远的概率越小，具体小多少，可以用一个正态分布曲线来模拟，这个分布曲线以直线对 Xi 给出的预测 f(Xi) 为中心，实际纵坐标为 Yi 的点 (Xi, Yi) 发生的概率就正比于 EXP[-(ΔYi)^2]。（EXP(..) 代表以常数 e 为底的多少次方）。

现在我们回到问题的贝叶斯方面，我们要想最大化的后验概率是：

P(h|D) ∝ P(h) * P(D|h)

又见贝叶斯！这里 h 就是指一条特定的直线，D 就是指这 N 个数据点。我们需要寻找一条直线 h 使得 P(h) * P(D|h) 最大。很显然，P(h) 这个先验概率是均匀的，因为哪条直线也不比另一条更优越。所以我们只需要看 P(D|h) 这一项，这一项是指这条直线生成这些数据点的概率，刚才说过了，生成数据点 (Xi, Yi) 的概率为 EXP[-(ΔYi)^2] 乘以一个常数。而 P(D|h) = P(d1|h) * P(d2|h) * .. 即假设各个数据点是独立生成的，所以可以把每个概率乘起来。于是生成 N 个数据点的概率为 EXP[-(ΔY1)^2] * EXP[-(ΔY2)^2] * EXP[-(ΔY3)^2] * .. = EXP{-[(ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + ..]} 最大化这个概率就是要最小化 (ΔY1)^2 + (ΔY2)^2 + (ΔY3)^2 + .. 。熟悉这个式子吗？

0 0