大家好，又见面了，我是你们的朋友全栈君。

光流法：Farnback

现实世界中，万物都在在运动，且运动的速度和方向可能均不同，这就构成了运动场。物体的运动投影在图像上反应的是像素的移动。这种像素的瞬时移动速度就是光流。光流法是利用图像序列中的像素在时间域上的变化、相邻帧之间的相关性来找到的上一帧跟当前帧间存在的对应关系，计算出相邻帧之间物体的运动信息的一种方法。

光流法按照不同的实现方式可以分为：基于梯度的方法、基于匹配的方法、基于能量的方法、基于相位的方法等。本文介绍的是一种基于梯度的经典光流方法：Farnback法。光流法的前提假设包括：相邻帧之间亮度恒定；相邻帧之间取时间连续或者运动变化微小；同一子图像中像素点具有相同的运动。

基本假设

假定图像序列记作$I(x,y,t)$,其中$X=[x,y]$。视频中的每个前后帧提取出来之后就是一个图像序列。假设图像亮度恒定，即图像亮度没有变化，则导数为0：

$$\frac{dI(X, t)}{d t}=\frac{\partial I}{\partial X}\frac{\partial X}{\partial t}+\frac{\partial I}{\partial t}=0$$

或者根据泰勒展开来得出上述式子：

$$I(X , t)=I(X + \Delta X, t + \Delta t) \approx I(X , t) + \frac{\partial I}{\partial X}\Delta X+\frac{\partial I}{\partial t}\Delta t \\ \frac{\partial I}{\partial X}\Delta X+\frac{\partial I}{\partial t}\Delta t=0\\ \frac{\partial I}{\partial X}\frac{\Delta X}{\Delta t}+\frac{\partial I}{\partial t}=0\\ \frac{\partial I}{\partial x}\frac{\partial x}{\partial t}+\frac{\partial I}{\partial y}\frac{\partial y}{\partial t}+\frac{\partial I}{\partial t}=0$$

其中，在微小时间内

∂X∂t ∂ X ∂ t
$\frac{\partial X}{\partial t}$或者

ΔXΔt Δ X Δ t
$\frac{\Delta X}{\Delta t}$表示速度，可以记为：

$$\frac{\partial X}{\partial t}=[\frac{\partial x}{\partial t}, \frac{\partial y}{\partial t}]=[u, v]=d$$

则有：

$$I_xu+I_yv+I_t=0\\\begin{bmatrix}I_x \ \ I_y\end{bmatrix}\begin{bmatrix}u \\ v\end{bmatrix}=-I_t$$

Farneback光流法

Farneback是一种基于梯度的方法，假设图像梯度恒定且假设局部光流恒定。局部光流恒定，即对于任意的$y\in N(x), d=\frac{\partial X}{\partial t} 不变$。梯度恒定即：

$$\frac{d}{dt}\nabla I(X, t)=\frac{\partial \nabla I}{\partial X}\frac{\partial X}{\partial t}+\frac{\partial \nabla I}{\partial t}=H(I) \cdot d +(\nabla I)_t=0$$

假设：

$$E(X, d)=||(H(I)\cdot d +(\nabla I)_t)||^2$$

上式在最优值处有导数为0：

$$\frac{\partial E}{\partial d}=0 \\ d=-\left(H^T(I)H(I)\right)^{-1}(H^T(I)(\nabla I)_t)$$

若对时间离散化：
（后向差分）

$$(\nabla I)_t(X, t) \approx \nabla I(X,t) - \nabla I(X,t-1)$$

（时间中心差分）

$$(\nabla I)_t(X, t-1/2) \approx \nabla I(X,t) - \nabla I(X,t-1)\\H(I)(X,t-1/2) \approx \frac{1}{2}(H(I)(X,t)+H(I)(X,t-1))$$

图像模型

图像一般是二维的（灰度图像），那么图像像素点的灰度值可以看成是一个二维变量的函数$f(x,y)$。假设以感兴趣的像素点为中心，构建一个局部坐标系（并不是针对整张图像）。对函数进行二项展开，可以近似为：

$$\begin{align}f(x,y) &\approx r_1+r_2x+r_3y+r_4x^2+r5y^2+r_6xy\\&=\left(\begin{matrix} x & y \end{matrix} \right)^T \left(\begin{matrix} r_{4} & r_{6}/2 \\ r_{6}/2 & r_{5} \end{matrix} \right) \left(\begin{matrix} x \\ y \end{matrix} \right)+ \left(\begin{matrix} r_{2} \\ r_{3} \end{matrix} \right)^T \left(\begin{matrix} x \\ y \end{matrix} \right)+r_{1}\\ &=\mathbf{x^TAx+b^Tx+}c \end{align}$$

其中，$\bf x$为二维列向量， $\bf A$为$2\times 2$的对称矩阵，$\bf b$为$2\times 1$的矩阵。注意，此处的系数确定后只针对在确定点$(x,y)$而言，对于其他点可能并不适用，也就是说，每个像素点对应一组系数。

取该像素点的一个邻域（通常以该像素为中心，大小为$2n+1的方形区域$），利用这些像素点的值和坐标来进行系数的估计，估计的算法可以使用加权最小二乘法。加权是因为在邻域内，距中心越近的像素点与中心像素具有更大的相关性，而越远的点提供的信息越少。其实可以将邻域以外的像素点的权重都看成是0。

位移估计

考虑多项式扩展是在一个像素的邻域内，如果像素经过移动$d$后，则整个多项式应该会发生变化。
原始位置：

$$f_1({\bf x})={\bf x^TA_1x+b_1^Tx}+c_1$$

像素移动后：

$$\begin{aligned}f_2({\bf x}) &=f_1({\bf x-d})\\ &={\bf (x-d)^TA_1(x-d)+b_1^T(x-d)}+c_1\\ &={\bf x^TA_1x+(b_1-2A_1d)^Tx + d^TA_1d-b_1^Td}+c_1\\ &={\bf x^TA_2x+b_2^Tx}+c_2 \end{aligned}$$

其中，

$${\bf A_2=A_1 \\b_2=b_1-2A_1d} \\c_2={\bf d^TA_1d-b_1^Td}+c_1$$

如果

A1 A 1
$\bf A_1$非奇异，则有上述的第二个式子可以得到：

$${\bf d}=-\frac{1}{2}\bf A_1^{-1}(b_2 - b_1)$$

按照理论推导，其中必定有

A1=A2 A 1 = A 2
$\bf A_1=A_2$,但实际情况中未必能满足这一项要求，因此可以通过来求平均来近似真实值。如果令：

$${\bf A(x)}=\frac{{\bf A_1(x) + A_2(x)}}{2}\\\Delta {\bf b(x)}=-\frac{1}{2} \bf(b_2 - b_1)$$

那么：

$$\bf A(x)d(x)=\Delta {\bf b(x)} \\d=(A^TA)^{-1}(A^T\Delta b)$$

可以构建目标函数来进行优化求得位移：

$$e({\bf x})=||{\bf Ad - \Delta b}||^2$$

实际情况中，这种方法求得的结果中噪声太多，因此可以使用兴趣像素点的邻域，然后使用加权的目标函数：

$$e({\bf x})=\sum_{\Delta x \in I} w(\Delta x)||{\bf A(x + \Delta x)d - \Delta b(x+\Delta x)}||^2$$

Reference

[1] 图像分析之光流之经典
[2] Farneback 光流算法详解与 calcOpticalFlowFarneback 源码分析
[3] 光流Optical Flow介绍与OpenCV的实现
[4] 光流法简单介绍
[5] Farneback, 2003, Two-Frame Motion Estimation Based on Polynomial Expansion