Horn–Schunck光流算法 — Horn–Schunck Method

Horn-Schunck光流算法通过引入全局平滑约束来做图像中的运动估计。Horn和Schunck设定图像中像素的运动速度和其临近像素的速度相似或相同，且光流场中的每处的速度变化是平滑的，不会突变【学习本文内容前，需要先了解光流的基础知识】。有两种方法来表示平滑约束（最小化以下式子）：
$$(\frac{\partial u}{\partial x}{)}^2+(\frac{\partial u}{\partial y}{)}^2和(\frac{\partial v}{\partial x}{)}^2+(\frac{\partial v}{\partial y}{)}^2$$

或者：
$${\nabla }^{2}u=\frac{{\partial }^{2}u}{{\partial }^{2}x}+\frac{{\partial }^{2}u}{{\partial }^{2}y}和{\nabla }^{2}v=\frac{{\partial }^{2}v}{{\partial }^{2}x}+\frac{{\partial }^{2}v}{{\partial }^{2}y}$$

对于${\nabla }^{2}u$和${\nabla }^{2}v$的计算，可以参考Horn和Schunck在1981年发表的论文$DeterminingOpticalFlow$中的计算方法。如下：
$${\nabla }^{2}u\approx {}_3({\bar{u}}_{i,j,k}-u_{i,j,k})和{\nabla }^{2}v\approx {}_3({\bar{v}}_{i,j,k}-v_{i,j,k})$$

其中，左下标$3$表示像素的临近范围，${\bar{u}}_{i,j,k}$和${\bar{v}}_{i,j,k}$定义如下：
u ˉ i , j , k = 1 6 { u i − 1 , j , k + u i , j + 1 , k + u i + 1 , j , k + u i , j − 1 , k } + 1 12 { u i − 1 , j − 1 , k + u i − 1 , j + 1 , k + u i + 1 , j + 1 , k + u i + 1 , j − 1 , k } , v ˉ i , j , k = 1 6 { v i − 1 , j , k + v i , j + 1 , k + v i + 1 , j , k + v i , j − 1 , k } + 1 12 { v i − 1 , j − 1 , k + v i − 1 , j + 1 , k + v i + 1 , j + 1 , k + v i + 1 , j − 1 , k } , \begin{aligned} \bar{u}_{i,j,k} = & \frac{1}{6}\{ u_{i-1,j,k} + u_{i,j+1,k} + u_{i+1,j,k} + u_{i,j-1,k}\} \\ & + \frac{1}{12}\{ u_{i-1,j-1,k} + u_{i-1,j+1,k} + u_{i+1,j+1,k} + u_{i+1,j-1,k}\}, \\ \bar{v}_{i,j,k} = & \frac{1}{6}\{ v_{i-1,j,k} + v_{i,j+1,k} + v_{i+1,j,k} + v_{i,j-1,k}\} \\ & + \frac{1}{12}\{ v_{i-1,j-1,k} + v_{i-1,j+1,k} + v_{i+1,j+1,k} + v_{i+1,j-1,k}\}, \\ \end{aligned} uˉi,j,k​=vˉi,j,k​=​61​{ui−1,j,k​+ui,j+1,k​+ui+1,j,k​+ui,j−1,k​}+121​{ui−1,j−1,k​+ui−1,j+1,k​+ui+1,j+1,k​+ui+1,j−1,k​},61​{vi−1,j,k​+vi,j+1,k​+vi+1,j,k​+vi,j−1,k​}+121​{vi−1,j−1,k​+vi−1,j+1,k​+vi+1,j+1,k​+vi+1,j−1,k​},​

如图：

至此，Horn-Schunck光流算法已经建立了两个约束条件，可以完成对光流中的$u$和$v$的计算。
$$\begin{array}{rl}& {\xi }_b=u{I}_x+v{I}_y+I_t\\ & {\xi }_c^2={\nabla }^{2}u+{\nabla }^{2}v\end{array}$$

基于上两式，Horn和Schunck构建误差函数${\xi }^2$，如下：
$${\xi }^2=\int \int ({\alpha }^2{\xi }_c^2+{\xi }_b^2)dxdy$$

其中，$\alpha$值越大，光流越平滑。目标是计算误差函数${\xi }^2$取得最小值时，$u$和$v$的取值。通过变分法（博主还不懂变分法，以后慢慢学～～），误差函数${\xi }^2$可变换为：
$$\begin{array}{r}{I}_x^{2}u+I_x{I}_{y}v={\alpha }^2{\nabla }^{2}u-I_x{I}_t\\ I_x{I}_{y}u+I_y^{2}v={\alpha }^2{\nabla }^{2}v-I_y{I}_t\end{array}$$

代入${\nabla }^{2}u={\bar{u}}_{i,j,k}-u_{i,j,k}$和${\nabla }^{2}v={\bar{u}}_{i,j,k}-u_{i,j,k}$，得：
$$\begin{gathered} ({\alpha }^2+I_x^2)u+I_x{I}_{y}v=({\alpha }^2\bar{u}-I_x{I}_t) \\ {I}_x{I}_{y}u+({\alpha }^2+I_y^2)v=({\alpha }^2\bar{v}-I_y{I}_t) \end{gathered}$$

上式的矩阵形式为：
$$\left[\begin{array}{cc}{\alpha }^2+I_x^2& I_x{I}_y\\ I_x{I}_y& {\alpha }^2+I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]=\left[\begin{array}{c}{\alpha }^2\bar{u}-I_x{I}_t\\ {\alpha }^2\bar{v}-I_y{I}_t\end{array}\right]$$

通过对上式进行行列式计算，得到$u$和$v$的表达式：
$$\begin{array}{rl}& ({\alpha }^2+I_x^2+I_y^2)u=({\alpha }^2+I_y^2)\bar{u}-I_x{I}_y\bar{v}-I_x{I}_t\\ & ({\alpha }^2+I_x^2+I_y^2)v=-I_x{I}_y\bar{u}+({\alpha }^2+I_x^2)\bar{v}-I_y{I}_t\end{array}$$

上式的另一种形式可以是：
$$\begin{gathered} ({\alpha }^2+I_x^2+I_y^2)(u-\bar{u})=-I_x[I_x\bar{u}+I_y\bar{v}+I_t] \\ ({\alpha }^2+I_x^2+I_y^2)(v-\bar{v})=-I_y[I_x\bar{u}+I_y\bar{v}+I_t] \end{gathered}$$

由于一开始，$u$和$v$的取值未知，因此$\bar{u}$和$\bar{v}$的取值也未知。无法直接通过上式计算$u$和$v$的取值。不过，我们可以采用迭代的方式来计算：
$$\begin{array}{r}{u}^{n+1}={\bar{u}}^n-\frac{I_x[I_x{\bar{u}}^n+I_y{\bar{v}}^n+I_t]}{{\alpha }^2+I_x^2+I_y^2}\\ v^{n+1}={\bar{v}}^n-\frac{I_y[I_x{\bar{u}}^n+I_y{\bar{v}}^n+I_t]}{{\alpha }^2+I_x^2+I_y^2}\end{array}$$

至此，Horn–Schunck光流算法原理已介绍完毕。