光流算法：Horn–Schunck光流讲解

光流理论

光流（optical flow）是空间运动物体在观察成像平面上的像素运动的瞬时速度。

光流法是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系，从而计算出相邻帧之间物体的运动信息的一种方法。

光流场是一个二维矢量场，它反映了图像上每一点灰度的变化趋势，可看成是带有灰度的像素点在图像平面上运动而产生的瞬时速度场。它包含的信息即是各像点的瞬时运动速度矢量信息。光流的集合。

基本假设

1. 亮度恒定不变。同一目标帧间变化亮度不发生改变。
2. 时间连续或运动为“小运动”。短时间内或者相邻帧间目标位置的变化不会剧烈。

基本假设对应的约束方程

像素点亮度在(x,y)位置以及(t)时刻表示为(I(x,y,t))，像素运动表示为(dx,dy,dt)，那么像素变化前后亮度相同表示为

[I(x,y,t) = I(x+dx,y+dy,t+dt) ]

泰勒展开：

[I(x+dx,y+dy,t+dt) = I(x,y,t) + frac{partial I}{partial x}dx + frac{partial I}{partial y}dy + frac{partial I}{partial t}dt + epsilon ]

其中(epsilon)是极小值，综合上式，可得

[frac{partial I}{partial x}dx + frac{partial I}{partial y}dy + frac{partial I}{partial t}dt = 0 ]

同除以(dt)，有

[frac{partial I}{partial x} frac{dx}{dt} + frac{partial I}{partial y} frac{dy}{dt} + frac{partial I}{partial t} = 0 ]

其中(u=frac{dx}{dt},v=frac{dy}{dt})为像素速度，也就是光流。令(I_x = frac{partial I}{partial x},I_y=frac{partial I}{partial y},I_t=frac{partial I}{partial t})，则

[I_xu+I_yv+I_t = 0 ]

((u,v))为所求光流。
但是求解((u,v))只有一个约束方程，是无法确定光流的，需要增加新的约束。不同的光流方法引入不同的约束条件，包括基于梯度的方法，基于匹配的方法，基于能量（频率）的方法，基于相位的方法和基于神经动力学的方法。

Horn–Schunck光流

Horn-Schunck光流算法通过引入全局平滑约束来做图像中的运动估计。Horn和Schunck设定图像中像素的运动速度和其临近像素的速度相似或相同，且光流场中的每处的速度变化是平滑的，不会突变。平滑约束表示为：

[ abla^2u = frac{partial^2u}{partial^2x} + frac{partial^2u}{partial^2y}, abla^2v = frac{partial^2v}{partial^2x} + frac{partial^2v}{partial^2y} ]

对于( abla^2u和 abla^2v)的处理，进行近似表示有

[ abla^2u simeq kappa (ar u_{i,j,k}-u_{i,j,k}), abla^2v simeq kappa (ar v_{i,j,k}-v_{i,j,k}) ]

其中领域均值表示为

[ar 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} } ]

[ar 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} } ]

约束写成：

[xi_c^2 = (ar u-u)^2+(ar v-v)^2 ]

综合计算

上面得到的两组约束为：

[xi_b = I_xu+I_yv+I_t ]

[xi_c^2 = (ar u-u)^2+(ar v-v)^2 ]

综合为：

[xi^2 = alpha^2xi_c^2+xi_b^2 ]

也就是最小化(xi^2)，其对应偏导为0即为所求

[frac{partial xi^2}{partial u} = -2alpha(ar u-u)+2(I_xu+I_yv+I_t)I_x ]

[frac{partial xi^2}{partial v} = -2alpha(ar v-v)+2(I_xu+I_yv+I_t)I_y ]

化简

[(alpha^2+I_x^2)u+I_xI_yv = (alpha^2 ar u-I_xI_t) \ I_xI_yu+(alpha^2+I_y^2)v= (alpha^2 ar u-I_yI_t) ]

得

[(alpha^2+I_x^2+I_y^2)u=+(alpha^2I_y^2)ar u -I_xI_yar v -I_xI_t \ (alpha^2+I_x^2+I_y^2)v=+(alpha^2I_x^2)ar v -I_xI_yar u -I_yI_t ]

或

[(alpha^2+I_x^2+I_y^2)(u-ar u) = -I_x[I_x ar u + I_y ar v + I_t] \ (alpha^2+I_x^2+I_y^2)(v-ar v) = -I_y[I_x ar u + I_y ar v + I_t] ]