Study of optical flow

记录几种常用的光流方法的原理和思想

参考了

Motion field

考虑在某个场景下，物体朝着三维空间的某个方向上移动，我们通过相机将物体的运动投影到图像屏幕的平面上形成视频，我们会从视频看到物体的运动方向。当物体在相机前移动时，图像中会有相应的更改。因此，如果对象上的点 $p_{o}$ 以速度 $v_{o}$ 移动，则可以将成像点 $p_{i}$ 看成以速度 $v_{i}$ 在图像平面上的运动。所有这些三维空间运动速度投影在二维图像平面的向量的集合形成了运动场 (motion field)。

然而我们想通过物体在视频中的移动速度来推导物体在三维空间中的移动速度是不现实的。原因如下，在一个很短的时间段 $δ t$ 内，场景中物体的速度为 $v_{0} = \frac{d r_{0}}{d t}$ 以及图像的速度 $v_{i} = \frac{d r_{i}}{d t}$ ，我们想要关联这两者的关系。根据位置关系我们能得到以下信息：

\frac{r_{i}}{f} = \frac{r_{0}}{r_{0} \cdot z}

其中 $r_{0} \cdot z$ 表示 $r_{0}$ 在 $z$ 轴上的投影。

图像 (运动场) 的速度为 (https://www.cse.psu.edu/~rtc12/CSE486/lecture22_6pp.pdf 有一些更详细的)：

v_{i} = \frac{d r_{i}}{d t} = \frac{\frac{f \cdot r_{0}}{r_{0} \cdot z}}{d t} = f \frac{(r_{0} \cdot z) v_{0} - (v_{0} \cdot z) r_{0}}{(r_{0} \cdot z)^{2}} = f \frac{(r_{0} \times v_{0}) \times z}{(r_{0} \cdot z)^{2}}

然而在绝大部分情况下，我们并不能直接通过上式来计算物体的运动场。退而求其次，我们只能通过计算视频光流（optical flow）对物体的运动场进行估计。在理想情况下，光流估计 = 运动场。然而，也有一些场景光流估计不等于运动场，比如下图中左图是一个单侧光源照射在一个表面质地相同的旋转的球体上，无论球体旋转到什么位置，其反射的光线在相机上成像都是相同的，此时物体的运动场存在，然而没有光流。下图右图中，存在一个围绕静止球体旋转的单侧光源，此时物体并没有旋转，而且其反射的光线在相机上的成像随着光源的旋转而改变，此时物体没有运动场存在，但是却有光流。

又如下图中，旋转灯的运动方向是横向朝右的，而其呈现出来的光流方向却是纵向向下的

同样，在下面两幅图中，图像的纹理并没有移动，而人眼却能观察到光流。

Optical flow

所谓的光流 (optical flow) 指的是物体或场景相对观测者（相机）运动引起的，往往将在连续帧图像上的特定坐标点上的亮度变化定义为该坐标点的光流矢量，对图像上每个坐标点都计算其光流矢量，则构成图像的光流场，我们希望用光流场来估计物体的运动场，计算光流的算法主要分为稀疏光流算法 (sparse optical flow) 和密集光流 (dense optical flow) 两类，常用的光流算法主要有以下几种：

Lucas-Kanade (LK) 光流法 (稀疏光流)
Horn-Schunck 算法

重要假设

基于光流的算法研究往往围绕着以下假设开始建模：

1. 光强恒定假设

2. 瞬时移动假设

考虑 $t$ 时刻，图像某点在坐标 $(x, y)$ 上的强度 (intensity) 为 $I (x, y, t)$ ，在一个很短的时间间隔 $δ t$ 后，该点移动到不远处坐标 $(x + δ x, y + δ y)$ ，其中 $δ$ 表示一个很小的值，其强度为 $I (x + δ x, y + δ y, t + δ t)$ ，在该光强恒定假设中，该点在不同时刻的强度恒定不变：

I (x, y, t) = I (x + δ x, y + δ y, t + δ t)

将右式进行一阶泰勒展开可得：

I (x + δ x, y + δ y, t + δ t) = I (x, y, t) + \frac{\partial I}{\partial x} δ x + \frac{\partial I}{\partial y} δ y + \frac{\partial I}{\partial t} δ t

PS：泰勒展开的原式为：

f (x) = \frac{x_{0}}{0!} + \frac{f^{'} (x_{0})}{1!} (x - x_{0}) + \frac{f^{″} (x_{0})}{2!} (x - x_{0})^{2} + . . . + \frac{f^{n} (x_{0})}{n!} (x - x_{0})^{n} + o (x - x_{0})

将一阶泰勒展开代入上式可以得到：

I (x, y, t) = I (x, y, t) + \frac{\partial I}{\partial x} δ x + \frac{\partial I}{\partial y} δ y + \frac{\partial I}{\partial t} δ t => \frac{\partial I}{\partial x} δ x + \frac{\partial I}{\partial y} δ y + \frac{\partial I}{\partial t} δ t = 0 => \frac{\partial I}{\partial x} \frac{δ x}{δ t} + \frac{\partial I}{\partial y} \frac{δ y}{δ t} + \frac{\partial I}{\partial t} \frac{δ t}{δ t} = 0

令 $u = \frac{δ x}{δ t}, v = \frac{δ y}{δ t}$ 分别表示光流分别为沿 $x$ 轴与 $y$ 轴的速度矢量， $I_{x} = \frac{\partial I}{\partial x}, I_{y} = \frac{\partial I}{\partial y}, I_{t} = \frac{\partial I}{\partial t}$ 分别表示为图像中像素点的灰度沿 $x, y, t$ 方向的偏导数。在具体操作中， $(I_{x}, I_{y}, I_{t})$ 可以很容易的从连续的两帧图像中计算出来。如下图所示，图像中像素点的灰度沿着 $x$ 方向的偏导数 $I_{x}$ 为右边四个方块减去左边四个方块的均值，同理 $I_{y}$ 为上面四个方块减去下面四个方块的均值， $I_{t}$ 为后面四个方块减去正面四个方块的均值。

我们现在有了这样的约束方程：

I_{x} u + I_{y} v + I_{t} = 0

如下图所示，由于 $I_{x}, I_{y}, I_{t}$ 已知，而 $u, v$ 未知，对于某点坐标的光流矢量 $u (u, v)$ 我们有下图二维坐标的黄线的约束，当然我们可以将该矢量分解为垂直于约束直线方向和平行于约束直线的方向的两个矢量 $u_{p}$ 和 $u_{n}$ 。

当然我们也可以轻易计算出垂直于约束直线方向的矢量的方向和大小：

因此问题的难点在于我们无法计算平行于约束直线的光流分量 $u_{p}$ 的大小，当然对于约束直线来说，存在一个方程和两个未知数，因此，仅有该约束，光流问题的求解是欠定的（under constrained），即该存在不同满足单一约束的解。从物理意义上，由于我们观察（假设）的是一个很小的区域的图像像素的灰度变化，如下图孔径问题（aperture problem）所示，在孔径中我们观察到的光流的方向是朝着右下角的方向，然而从整体观察来看，其物体实际运动可能是向下或者是向右的。

因此我们需要额外的约束来解决光流估计问题，也就是不同光流估计算法提出的出发点和需要解决的问题。

Lucas-Kanade (LK) 光流法

LK 光流法，进一步假设，图像中一个小的窗口中的每个像素，其运动场也就是光流是相同的。用数学的语言表达，在一个大小为 $n \times n$ 的窗口 $W$ 中的所有像素点 $(k, l) \in W$ ，其满足：

I_{x} (k, l) u + I_{y} (k, l) v + I_{t} (k, l) = 0

使用矩阵的形式表达如下：

[\begin{matrix} I_{x} (1, 1) & I_{y} (1, 1) \\ I_{x} (k, l) & I_{y} (k, l) \\ ⋮ & ⋮ \\ I_{x} (n, n) & I_{y} (n, n) \end{matrix}] [\begin{matrix} u \\ v \end{matrix}] = [\begin{matrix} I_{t} (1, 1) \\ I_{t} (k, l) \\ ⋮ \\ I_{t} (n, n) \end{matrix}] A u = B

此时约束方程数 ( $n \times n$ ) 大于未知数 (2)，该问题会变成一个超定问题，解决该超定问题一种方式看成是解线性最小二乘拟合，使得偏差平方和最小。一种快速和简单的方法如下，将上式两边同时左乘上 $A$ 的转置 $A^{T}$ 得：

[\begin{matrix} \sum_{W} I_{x} I_{x} & \sum_{W} I_{x} I_{y} \\ \sum_{W} I_{x} I_{y} & \sum_{W} I_{y} I_{y} \end{matrix}] [\begin{matrix} u \\ v \end{matrix}] = [\begin{matrix} - \sum_{W} I_{x} I_{t} \\ - \sum_{W} I_{y} I_{t} \end{matrix}] A^{T} A u = A^{T} B

可以解得：

u = (A^{T} A)^{- 1} A^{T} B

其中， $A^{T} A$ 是一个 $2 \times 2$ 的矩阵，从拟合效果来说，上式成立存在两个条件：

$A^{T} A $ $ 必须是可逆的 (i n v e r t i b l e) ，即其行列式 $ $ | A^{T} A | \neq 0$
$A^{T} A $ $ 必须是适定的 (w e l l - c o n d i t i o n e d) ，如果 $ $ λ_{1} $ $ 和 $ $ λ_{2} $ $ 是 $ $ A^{T} A $ $ 的特征值，需要满足： - $ $ λ_{1} > ϵ $ $ 并且 $ $ λ_{2} > ϵ$
- $λ_{1} \geq λ_{2} $ $ 并且 $ $ λ_{1} ≫̸ λ_{2}$

如果 $λ_{1}$ 和 $λ_{2}$ 都很小，这说明矩阵的基向量的重要程度很低。如果其中一个特征值 $λ_{1}$ 远大于另一个特征值 $λ_{2}$ ，则说明矩阵在一个方向上的基向量很重要，而另一个方向上的基向量重要程度很低。如果这两个条件之一不满足，其求解出来的光流是不准确的。

从实际的物理意义来看，如果窗口定位在光滑区域上，由于光滑区域的各方向的梯度都很小，则其矩阵的 $λ_{1}$ 和 $λ_{2}$ 都很小，因此估计的光流是不可靠的。

当另一个条件不满足时，其对应的实际物理意义是，窗口定位在边缘的区域上，则垂直于边的梯度很大，平行于边的梯度很小，这样会产生矩阵一个特征值 $λ_{1}$ 远大于另一个特征值 $λ_{2}$ 的结果，从上面的图可以得知，这样估计得到的光流也是不可靠的。说明了即使是适定问题，其任然有可能是个病态问题，体现在其输入端极小的误差，也可能导致输出值的巨大误差。

那么什么样的区域对于光流的估计是更为准确的呢？结果是窗口定位在富含纹理的区域，这种区域各方向的梯度都比较大而且多样，得到的特征值都比较大，估计的光流较为准确。

当然物理意义上也可以认为是 Harris 角点检测成立的条件，本质上是用 $A^{T} A$ 的特征值 $λ_{1} 、 λ_{2}$ 来分类图像中平滑、边缘、角点的区域。

因此 LK 算法在平滑和边缘的区域表现较差，而在角点或者梯度丰富的区域表现得较好。

Horn-Schunck(HS) 光流法

More details: http://www.ipol.im/pub/art/2013/20/

HS 光流法是最早提出来估计光流的方法（之一？），他把光强恒定的假设写成需要优化的损失函数，并加上一个平滑约束，新的约束的出发点是运动场在图像上趋于平稳变化，即速度的变化率趋于零，在该假设下，相邻帧的图像中光流分别为沿 $x$ 轴与 $y$ 轴的速度矢量 $u$ 和 $v$ 随着像素点移动而发生的改变是缓慢的，如果变化过大则给于惩罚， $λ$ 是一个用于权衡重要性的超参数：

L = \iint_{i m a g e} (I_{x} u + I_{y} v + I_{t})^{2} + λ (u_{x}^{2} + u_{y}^{2} + v_{x}^{2} + v_{y}^{2}) d x d y L = \iint_{i m a g e} E_{d a t a} (u, v) + λ E_{s m o o t h n e s s} (u, v) d x d y

目标是为了最小化上述损失函数，其中 $u_{x}$ 是速度矢量 $u$ 对 $x$ 的偏导 $\frac{\partial u}{\partial x}$ ，同理 $u_{y} = \frac{\partial u}{\partial v}, u_{x} = \frac{\partial v}{\partial x}, u_{y} = \frac{\partial u}{\partial y}$ ，由此看出 HS 光流法给出的是一个全局约束。

我们可以通过多维欧拉 - 拉格朗日方程

对其进行求解有：

\frac{\partial L}{\partial u} - \frac{\partial}{\partial x} \frac{\partial L}{\partial u_{x}} - \frac{\partial}{\partial y} \frac{\partial L}{\partial u_{y}} = 0 \frac{\partial L}{\partial v} - \frac{\partial}{\partial x} \frac{\partial L}{\partial v_{x}} - \frac{\partial}{\partial y} \frac{\partial L}{\partial v_{y}} = 0

其中 $L = (I_{x} u + I_{y} v + I_{t})^{2} + λ (u_{x}^{2} + u_{y}^{2} + v_{x}^{2} + v_{y}^{2})$ ，代入上式后得：

I_{x} (I_{x} u + I_{y} v + I_{t}) - α^{2} Δ u = 0 I_{y} (I_{x} u + I_{y} v + I_{t}) - α^{2} Δ v = 0

其中 $Δ = \frac{\partial^{2}}{\partial x^{2}} + \frac{\partial^{2}}{\partial y^{2}}$ 是拉普拉斯算子，在实际中，拉普拉斯算子是用有限差分进行数值逼近的，常被写为 $Δ u (x, y) = 4 (\overset{―}{u} (x, y) - u (x, y))$ ，其中 $\overset{―}{u} (x, y)$ 是一个加权平均模板：

\frac{1}{12} (\begin{matrix} 1 & 2 & 1 \\ 2 & - 12 & 2 \\ 1 & 2 & 1 \end{matrix})

用该式替换上式后得：

(I_{x}^{2} + 4 α^{2}) u + I_{x} I_{y} v = 4 α^{2} \overset{―}{u} - I_{x} I_{t} I_{x} I_{y} u + (I_{y}^{2} + 4 α^{2}) v = 4 α^{2} \overset{―}{v} - I_{y} I_{t}

这是一个二元一次线性方程组，可以对图像中的每个像素求解。然而，由于解依赖于光流场的邻近值，所以一旦邻近值被更新，解就必须迭代更新，利用 Cramer 规则可以推导出一下迭代公式：

u^{k + 1} = u^{- k} - \frac{I_{x} (I_{x} u^{- k} + I_{y} v^{- k} + I_{t})}{4 α^{2} + I_{x}^{2} + I_{y}^{2}} v^{k + 1} = v^{- k} - \frac{I_{y} (I_{x} u^{- k} + I_{y} v^{- k} + I_{t})}{4 α^{2} + I_{x}^{2} + I_{y}^{2}}

Hierarchical method (金字塔改进）

上述的方法中存在的问题是：

鲁棒性较差，无法准确处理原本图象梯度较大的边界区域；
并没有解决我们先前提出的针对变化幅度比较大情况下如何准确求解光流的问题。于是乎有了 Coarse-to-Fine 方法：

针对光流估计的算法还提出了一个改进算法 (Coarse-to-Fine)，其要解决的问题是，当物体运动速度很快的时候，即连续两帧物体位置差异较大的情况，那么，在完整的分辨率下，光强恒定和瞬时移动的假设将不成立。这个改进方法的想法是，通过降分辨率的方式来使假设成立，因为在低分辨率情况下，物体运动的位置变化变小了，打个比方，物体在原始分辨率图像下运动 4px 的像素，在 1/4 分辨率大小的图像上只运动 1px 的像素，因此在分辨率降低的情况下可以使得原来的光强恒定和瞬时移动的假设成立。为了计算原始分辨率图像下的光流场，如下图所示，先在低分辨率情况下使用光流算法得到图像的光流场，使用仿射变换将光流场转换到中等分辨率的情况下，利用中等分辨率大小的光流场和前一帧的图像计算得到一个中间插值 (interpolation) 图像，利用插值帧和后一帧的图像计算中等分辨率下的光流场。以此类推，得到原始分辨率下的光流场。