金字塔LK光流法

金字塔LK光流法

最近看的一篇论文中有金字塔LK光流法，于是看了些东西，整理一下。

光流法

光流（optical flow）是空间运动物体在观察成像平面上的像素运动的瞬时速度。光流法是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。通常将二维图像平面特定坐标点上的灰度瞬时变化率定义为光流矢量。
光流计算基于物体的运动的光学特性，提出了两个假设：

1. 运动物体的灰度值在很短的时间间隔内保持不变；
2. 给定邻域内的速度向量变化缓慢；

假设一个像素点为$(x,y)$，$I(x,y,t)$为该像素点$t$时刻的的亮度，用$u(x,y)$和$v(x,y)$来表示该像素点在水平和竖直方向上的速度分量：
$$u=\frac{dx}{dt};v=\frac{dy}{dt}$$
那么，在时间间隔很短的情况下（比如相邻两帧），$u$和$v$便可以看作各自方向上的位移。
我们假设时间间隔为$\Delta t$，那么此时的亮度即为$I(x+\Delta x,y+\Delta y,t+\Delta t)$。在假设$\Delta t$很小的前提下，将其泰勒展开，并舍弃掉二阶无穷小项：$$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\text{\hspace{1em}(1)}$$
根据假设（1）：运动物体的灰度值在很短的时间间隔内保持不变，可得$I(x+dx,y+dy,t+dt)=I(x,y,t)$，所以公式（1）中的$\frac{\partial I}{\partial x}\Delta x+\frac{\partial I}{\partial y}\Delta y+\frac{\partial I}{\partial t}\Delta t=0$，即$$\begin{array}{c}-\frac{\partial I}{\partial t}=\frac{\partial I}{\partial x}\frac{dx}{dt}+\frac{\partial I}{\partial y}\frac{dy}{dt}=\frac{\partial I}{\partial x}u+\frac{\partial I}{\partial y}v\\ -I_t=I_{x}u+I_{y}v\\ -I_t=\left[\begin{array}{ll}{I}_x& I_y\end{array}\right]\left[\begin{array}{l}u\\ v\end{array}\right]\end{array}\text{\hspace{1em}(2)}$$
其中，$I_x$和$I_y$是图像在$(x,y)$处的梯度，$I_x=I(x+1,y)-I(x,y)$，$I_y=I(x,y+1)-I(x,y)$,$I_t$是图像在$(x,y)$处关于时间$t$的导数，$I_t\approx (I(x,y,t)-I(x,y,t-1))$。
这就是基本的光流约束方程。但此时要求出$u$和$v$，却只有一个方程。一个方程两个未知量，这是没法求解的。于是便有了LK光流法。

LK光流法

LK光流法提出了第三个假设：邻域内光流一致，即一个场景中的同一表面的局部邻域内具有相似的运动，在图像平面上的投影也在邻近区域，且邻近点速度一致（邻域内所有像素点光流一样）。
这时，邻域$w\times w$内所有像素点（假设为$k$个）运动一致，那么公式（2）便有$k$个：$$\{\begin{array}{l}{I}_{1x}u+I_{1y}v=-I_{1t}\\ I_{2x}u+I_{2y}v=-I_{2t}\\ \cdots \\ I_{kx}u+I_{ky}v=-I_{kt}\end{array}$$
但$k$个方程，两个未知量，也是无法求解的。所以只能找出该方程的最优解。
记：
$$A=\left[\begin{array}{c}{\left[I_x,I_y\right]}_1\\ \cdots \\ {\left[I_x,I_y\right]}_k\end{array}\right],b=\left[\begin{array}{c}{I}_{t1}\\ \cdots \\ I_{tk}\end{array}\right]$$
则方程可以写为：$A\left[\begin{array}{l}u\\ v\end{array}\right]=-b$，这里写最小二乘法的解法，还可以用牛顿迭代法等，便不做推导。将方程写成：$A\vec{x}=\vec{z}$，得到：
$$\begin{array}{l}{A}^{T}A\vec{x}=A^T\vec{z}\\ \vec{x}={\left(A^{T}A\right)}^{-1}{A}^T\vec{z}\end{array}$$
即：$$\vec{x}=\left[\begin{array}{c}u\\ v\end{array}\right]={\left[\begin{array}{cc}{\sum }_{i=1}^k{I}_{ix}^2& {\sum }_{i=1}^k{I}_{ix}{I}_{iy}\\ {\sum }_{i=1}^k{I}_{ix}{I}_{iy}& {\sum }_{i=1}^k{I}_{iy}^2\end{array}\right]}^{-1}\left[\begin{array}{c}-{\sum }_{i=1}^k{I}_{ix}{I}_t\\ -{\sum }_{i=1}^k{I}_{iy}{I}_t\end{array}\right]$$
但只有在$A^{T}A$可逆的时候，才能得出答案。在图像中，沿着两个方向都有像素发生变化的区域，$A^{T}A$才可逆；反之，在灰度变化很小的区域，$A^{T}A$一般不可逆。这限制了LK光流法的应用范围，也是其被称为“稀疏光流法”的主要原因。

迭代求解LK法

整理一下求解过程。
相邻的两帧图像$I$和$J$，对于$I$中的像素点$u={\left[\begin{array}{ll}{u}_x& u_y\end{array}\right]}^T$，需要在图像$J$中找到像素点$v=u+d={\left[\begin{array}{ll}{u}_x+d_x& u_y+d_y\end{array}\right]}^T$使其与前一个像素点最为相似。我们把$d={\left[\begin{array}{ll}{d}_x& d_y\end{array}\right]}^T$称为$u$点的光流。为了求解，假设邻域内的点具有相同的光流，邻域采用$w_x$和$w_y$两个参数表示。则求解$d$便成了使以下目标函数最小的优化问题：$$\epsilon (d)=\epsilon \left(d_x,d_y\right)=\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}{\left(I(x,y)-J\left(x+d_x,y+d_y\right)\right)}^2$$
最优解的导数为$0$，则有：
$$\frac{\partial \epsilon (d)}{\partial d}=-2\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left(I(x,y)-J\left(x+d_x,y+d_y\right)\right)\left[\frac{\partial J}{\partial x}\frac{\partial J}{\partial y}\right]$$
一阶泰勒展开后：
$$\frac{\partial \epsilon (d)}{\partial d}\approx -2\sum _{x=u_x-w_{x}y=u_y-w_y}^{x=u_x+w_{x}y=u_y+w_y}\left(I(x,y)-J(x,y)-\left[\frac{\partial J}{\partial x}\frac{\partial J}{\partial y}\right]\left[\begin{array}{l}{d}_x\\ d_y\end{array}\right]\right)\left[\begin{array}{ll}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]$$
由于时间间隔很短，所以$d={\left[\begin{array}{ll}{d}_x& d_y\end{array}\right]}^T$足够小，就可以用$\left[\begin{array}{ll}\frac{\partial I}{\partial x}& \frac{\partial I}{\partial y}\end{array}\right]$来代替$\left[\begin{array}{ll}\frac{\partial J}{\partial x}& \frac{\partial J}{\partial y}\end{array}\right]$。同时，为了使式子看上去简单一些，便定义$\delta I(x,y)=I(x,y)-J(x,y)$和$\nabla I={\left[\begin{array}{ll}\frac{\partial I}{\partial x}& \frac{\partial I}{\partial y}\end{array}\right]}^T$，所以上述公式可以转换为：$$\frac{1}{2}\frac{\partial \epsilon (d)}{\partial d}=\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left(\nabla I^{T}d-\delta I\right)\nabla I^T$$
其中，$\nabla I^{T}d-\delta I$是标量，所以可以将式子中的$\nabla I^T$换成$\nabla I$，这样变换之后，虽然最后的结果从$1\times 2$的向量变成了$2\times 1$的向量，但值是一样的。所以：$$\begin{array}{rl}\frac{1}{2}{\left[\frac{\partial \epsilon (d)}{\partial d}\right]}^T& \approx \sum _{x=u_x-w_y}^{x=u_x+w_{x}y=u_y-u_y-w_y}\left(\nabla I^{T}d-\delta I\right)\nabla I\\ & =\sum _{x=u_x-w_y}^{x=u_x+w_{x}y=u_y-u_y-w_y}\nabla I^{T}d\nabla I-\delta I\nabla I\\ & =\sum _{x=u_x-w_y}^{x=u_x+w_{x}y=u_y-u_y-w_y}\nabla I\nabla I^{T}d-\delta I\nabla I\\ & =\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]d-\sum _{x=u_x-w_x}^{x=u_x+w_y}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left[\begin{array}{l}\delta I\cdot I_x\\ \delta I\cdot I_y\end{array}\right]\end{array}$$
令$$\begin{array}{rl}G& =\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right]d\\ b& =\sum _{x=u_x-w_x}^{x=u_x+w_y}\sum _{y=u_y-w_y}^{y=u_y+w_y}\left[\begin{array}{l}\delta I\cdot I_x\\ \delta I\cdot I_y\end{array}\right]\end{array}$$
则公式最终可以写成$d=G^{-1}b$的形式。有时候为了求得更加精确的$d$，可以采用迭代求解的方式求$d$。得到的位置$(x+d_x,y+d_y)$往往不是整数点像素，需要采用类似双线性插值的方法来计算得到最终的值。
用$K(K\ge 1)$来表示迭代次数，对于第$k(K\ge k\ge 1)$次迭代，前面的$k-1$次迭代提供了初始值$d^{k-1}={\left[\begin{array}{ll}{d}_x^{k-1}& d_y^{k-1}\end{array}\right]}^T$，移动后得到的图像为$J_k(x,y)=J\left(x+d_x^{k-1},y+d_y^{k-1}\right)$，第k次迭代需要求解的问题就是计算优化变量${\eta }^k={\left[\begin{array}{ll}{\eta }_x^k& {\eta }_y^k\end{array}\right]}^T$，使得下面的目标函数最小：$${\epsilon }^k\left({\eta }^k\right)=\epsilon \left({\eta }_x^k,{\eta }_y^k\right)=\sum _{x=u_x-w_x}^{x=u_x+w_x}\sum _{y=u_y-w_y}^{y=u_y+w_y}{\left(I(x,y)-J_k\left(x+{\eta }_x^k,y+{\eta }_y^k\right)\right)}^2$$
也可以写成${\eta }^k=G^{-1}{b}_k$的形式。这里的$G$只需要计算一次，而$b_k$每一次迭代都不一样，需要重新计算。之后第$k$次迭代的结果$d_k=d_{k-1}+{\eta }_k$

金字塔LK光流法

之前的三条假设分别为：

1. 运动物体的灰度值在很短的时间间隔内保持不变；
2. 给定邻域内的速度向量变化缓慢；
3. 邻域内光流一致；

这些假设在实际场景中很难满足，尤其是假设2。我们在求解过程中应用了泰勒展开，泰勒展开只有在变量变化很小的情况下才能使用，而如果帧之间的像素运动比较大，泰勒展开便不怎么合适了。所以就有大神提出了金字塔LK光流法，既然害怕像素运动太快，无法使用泰勒展开，那么就将整张图片进行缩小，降低其分辨率。对于运动较快的像素点，总能在图像分辨率降到一定程度时，其运动变得足够小，满足泰勒展开的条件。
算法将原始图像作为第0层，宽和高缩小$2^L$倍的图像作为第$L$层，形成金字塔的样子。一些小细节：

1. 对图像进行降采样之前通常采用低通滤波器进行滤波，防止降采样后发生锯齿现象；
2. 对图像添加额外像素圈，将额外像素圈内的像素值填充为图像的真实边界值，这样就算在计算边缘点邻域时，只计算该点邻域的有效部分；
3. 原始图像的宽和高不严格满足可以整除$2^L$，可以微调原图像，也可以按照$n_x^L\le \frac{n_x^{L-1}+1}{2},n_y^L\le \frac{n_y^{L-1}+1}{2}$来取宽和高；

算法大致流程其实就是先计算最顶层的光流大小，定义为$g^{Lm}={\left[\begin{array}{ll}0& 0\end{array}\right]}^T$，传入到下一层作为初始值，最后计算出原始图像上的光流大小。假设第$L+1$层计算到的第$L$层的光流大小为$g^L={\left[\begin{array}{ll}{g}_x^L& g_y^L\end{array}\right]}^{\prime }$，需要求得$d^L=\left[\begin{array}{ll}{d}_x^L& d_y^L\end{array}\right]$使以下函数最小：
$${\epsilon }^L\left(d^L\right)={\epsilon }^L\left(d_x^L,d_y^L\right)=\sum _{x=u_x^L-w_{x}y=u_y^L-w_y}^{x=u_x^L+w_{x}y=u_y^L+w_y}{\left(I^L(x,y)-J^L\left(x+g_x^L+d_x^L,y+g_y^L+d_y^L\right)\right)}^2$$
由于泰勒展开是在变化微小的情况下才能使用，所以此处需要$x+g_x^L$才能接近真实值，才能泰勒展开，相对于把邻域窗口进行了平移。后面便是和之前LK光流法类似，不再叙述。
计算得到$d^L=\left[\begin{array}{ll}{d}_x^L& d_y^L\end{array}\right]$后，定义$g^{L-1}=2\left(g^L+d^L\right)$作为第$L-1$层的光流初始值大小，重复上述过程，最终计算得到的结果为$d={\sum }_{L=0}^{Lm}{2}^L{d}^L$。虽然对于原始图像来说，位移可能偏大，但对于每一层的图像来说，其对应的$d^L$其实都很小，因此可以使用泰勒展开求解。

整理一下整体流程

将图像$I$和$J$建立为${\left\{I^L\right\}}_{L=0,\dots L_m}$和${\left\{J^L\right\}}_{L=0,\dots L_m}$；
初始化最高层的光流为$g^{Lm}={\left[\begin{array}{ll}0& 0\end{array}\right]}^T$；
从$L=L_m$开始进行迭代：
   计算图像$I^L$中像素点$u$的位置：$u^L=\left[\begin{array}{ll}{u}_x^L& u_y^L\end{array}\right]=u/2^L$
   计算图像$I^L$在$x$方向上的梯度：$I_x^L(x,y)=\frac{I^L(x+1,y)-I^L(x-1,y)}{2}$（中心差分法）
   计算图像$I^L$在$y$方向上的梯度：$I_y^L(x,y)=\frac{I^L(x,y+1)-I^L(x,y-1)}{2}$（中心差分法）
   计算矩阵G：$$G=\sum _{x=u^{L}x-w_x}^{x=u^{L}x+w_x}\sum _{y=u^{L}y-w_y}^{y=u^{L}y+w_y}\left[\begin{array}{cc}{I}_x^{L^2}& I_x^L{I}_y^L\\ I_x^L{I}_y^L& I_y^{L^2}\end{array}\right]$$
   初始化迭代$k$次的初始值：$d^0={\left[\begin{array}{ll}0& 0\end{array}\right]}^T$
   进行k次迭代
     计算图像差异：$\delta I^k(x,y)=I^L(x,y)-J^L\left(x+g_x^L+d_x^{k-1},y+g_y^L+d_y^{k-1}\right)$
     计算向量$b_k$：$b_k={\sum }_{x=u^{L}x-w_x}^{x=u^{L}x+w_x}{\sum }_{y=u^{L}y-w_y}^{y=u^{L}y+w_y}\left[\begin{array}{l}\delta I^k(x,y)I_x(x,y)\\ \delta I^k(x,y)I_y(x,y)\end{array}\right]$
     计算光流：${\eta }_k=G^{-1}{b}_k$
     计算下一次迭代的初始值：$d^k=d^{k-1}+{\eta }^k$
   结束$k$上的迭代
   得到$L$层图像上的光流：$d^L=d^k$
  为$L-1$层提供光流的初始值：$g^{L-1}=2\left(g^L+d^L\right)$
结束$L$上的迭代
计算最终的光流$d=g^0+d^0$
找到图像$J$中对应图像$I$中像素点$u$的位置：$v=u+d$