本文主要对 TV-L1 光流算法的计算流程进行了梳理，不涉及算法原理的详细推理。在阅读论文时可参考以下的开源代码。
论文：Zach, C., Pock, T., & Bischof, H. (2007). A duality based approach for realtime tv-l 1 optical flow. In Pattern Recognition: 29th DAGM Symposium, Heidelberg, Germany, September 12-14, 2007. Proceedings 29 (pp. 214-223). Springer Berlin Heidelberg.
代码：https://github.com/pariasm/bwd-nlkalman/tree/master/lib/tvl1flow
1 TV-L1算法流程
记输入两帧 N N N 维单通道图像分别为 I 0 {I_0} I0 和 I 1 {I_1} I1， N N N 维坐标矢量为 x = ( x 1 , . . . , x N ) T \bm{x}={\left( {{x_1},...,{x_N}} \right)^T} x=(x1,...,xN)T，对应的 N N N 维光流矢量为 u = ( u 1 , . . . , u N ) T \bm{u} = {\left( {{u_1},...,{u_N}} \right)^T} u=(u1,...,uN)T，则带有 L 1 {L_1} L1 范数图像保真度惩罚项的总变分（Total Variation, TV）正则化代价函数的定义为
E = ∫ Ω { λ ∣ I 0 ( x ) − I 1 ( x + u ( x ) ) ∣ + ∣ ∇ u ∣ } d x , (1.1) E = \int_\Omega {\left\{ {\lambda \left| {{I_0}\left( \bm{x} \right) - {I_1}\left( {\bm{x} + \bm{u}\left( \bm{x} \right)} \right)} \right| + \left| {\nabla \bm{u}} \right|} \right\}d\bm{x}} ,\tag{1.1} E=∫Ω{λ∣I0(x)−I1(x+u(x))∣+∣∇u∣}dx,(1.1)
其中 Ω \Omega Ω 为全体图像坐标， λ \lambda λ 为 L 1 {L_1} L1 惩罚项的权重， ∇ \nabla ∇ 为梯度算子（读nabla），且有
∣ ∇ u ∣ = ∑ d = 1 N ∣ ∇ u d ∣ 2 = ∑ d = 1 N ∑ i = 1 N ( ∂ u d ∂ x i ) 2 . (1.2) \left| {\nabla \bm{u}} \right| = \sqrt {\sum\limits_{d = 1}^N {{{\left| {\nabla {u_d}} \right|}^2}} } = \sqrt {\sum\limits_{d = 1}^N {\sum\limits_{i = 1}^N {{{\left( {\frac{{\partial {u_d}}}{{\partial {x_i}}}} \right)}^2}} } } .\tag{1.2} ∣∇u∣=d=1∑N∣∇ud∣2=d=1∑Ni=1∑N(∂xi∂ud)2.(1.2)
式 (1.1) 积分项前半部分的意义是，基于所求得光流矢量 u \bm{u} u 对图像 I 0 {I_0} I0 进行运动补偿映射后所得图像应该与图像 I 1 {I_1} I1 尽可能保存一致，这是图像保真度方面的要求；而后半部分的意义是，在局部区域中不同像素位置的光流变化应该是相对平缓的，这是物体作为一个整体运动该有的性质。
1.1 讨论只有水平运动的情况
为了便于理解，我们以常见的二维图像为例，并假设只存在水平单一方向的运动，此时可简化地用标量 u ( x ) u\left( \bm{x} \right) u(x) 来表示光流，这里不太严谨地用 x + u ( x ) \bm{x} + u\left( \bm{x} \right) x+u(x) 代表 x + ( u ( x ) , 0 ) T \bm{x} + {\left( {u\left( \bm{x} \right),0} \right)^T} x+(u(x),0)T。但要注意，虽然我们假设只存在水平方向的运动，但不代表所有像素的水平运动速度都是一样，即光流函数 u ( x ) u\left( \bm{x} \right) u(x) 在水平或者垂直方向上都存在着大小变化，都可能具有非零的梯度值。在大多数光流算法中，我们常假设图像局部区域的像素值是低频或者线性变化的（如渐变等），这时我们可通过泰勒公式展开对 I 1 {I_1} I1 进行线性近似，即
I 1 ( x + u ( x ) ) = I 1 ( x ) + u ( x ) I 1 x ( x ) , (1.3) {I_1}\left( {\bm{x} + u\left( \bm{x} \right)} \right) = {I_1}\left( \bm{x} \right) + u\left( \bm{x} \right)I_1^x\left( \bm{x} \right),\tag{1.3} I1(x+u(x))=I1(x)+u(x)I1x(x),(1.3)
其中 I 1 x ( x ) I_1^x\left( \bm{x} \right) I1x(x) 为图像 I 1 {I_1} I1 水平方向的偏导数，这里我们常通过中心差分来估计，即
I 1 x ( x , y ) = 1 2 I 1 ( x + 1 , y ) − I 1 ( x − 1 , y ) . (1.4) I_1^x\left( {x,y} \right) = \frac{1}{2}{I_1}\left( {x + 1,y} \right) - {I_1}\left( {x - 1,y} \right).\tag{1.4} I1x(x,y)=21I1(x+1,y)−I1(x−1,y).(1.4)
对于式 (1.3) 的线性近似，当光流 u ( x ) u\left( \bm{x} \right) u(x) 较大时会存在较大的误差，这也是显而易见的，因为越大的邻域更有可能包含更多的异常区域。但问题是我们事先并不知道各像素的光流大小，也就无从保证线性近似的准确性。一种直观的解决办法是通过多尺度迭代来实现由粗到细（Coarse-to-Fine）的光流修正，即首先将原图像缩小到线性化误差较小的尺度，计算得到此尺度下的光流值，然后将所得光流值放大至原来的尺度，以此光流值作为引导计算出更精确的光流值；同样地，在同一尺度内，我们也可利用已求得的光流值作为引导，再进行一次或多次的迭代。关于多尺度的光流计算方法我们且留到后面再表，这里假设已得到了用于引导计算的初始光流 u 0 ( x ) {u_0}\left( \bm{x} \right) u0(x)，则式 (1.3) 的线性化可修改为
I 1 ( x + u ( x ) ) = I 1 ( x + u 0 ( x ) ) + ( u ( x ) − u 0 ( x ) ) I 1 x ( x + u 0 ( x ) ) . (1.5) {I_1}\left( {\bm{x} + u\left( \bm{x} \right)} \right) = {I_1}\left( {\bm{x + }{u_0}\left( \bm{x} \right)} \right) + \left( {u\left( \bm{x} \right) - {u_0}\left( \bm{x} \right)} \right)I_1^x\left( {\bm{x + }{u_0}\left( \bm{x} \right)} \right).\tag{1.5} I1(x+u(x))=I1(x+u0(x))+(u(x)−u0(x))I1x(x+u0(x)).(1.5)
为了保持公式的简洁，后面不妨省略 u ( x ) u\left( \bm{x} \right) u(x)的坐标符号。于是式 (1.1) 的代价函数可近似为
E = ∫ Ω { λ ∣ u I 1 x ( x + u 0 ) + I 1 ( x + u 0 ) − u 0 I 1 x ( x + u 0 ) − I 0 ( x ) ∣ + ∣ ∇ u ∣ } d x . (1.6) E = \int_\Omega {\left\{ {\lambda \left| {uI_1^x\left( {\bm{x} + {u_0}} \right) + {I_1}\left( {\bm{x + }{u_0}} \right) - {u_0}I_1^x\left( {\bm{x} + {u_0}} \right) - {I_0}\left( \bm{x} \right)} \right| + \left| {\nabla u} \right|} \right\}d\bm{x}} .\tag{1.6} E=∫Ω{λ∣uI1x(x+u0)+I1(x+u0)−u0I1x(x+u0)−I0(x)∣+∣∇u∣}dx.(1.6)
为了方便，我们把图像残差部分表示为
ρ ( u , u 0 , x ) = u I 1 x ( x + u 0 ) + I 1 ( x + u 0 ) − u 0 I 1 x ( x + u 0 ) − I 0 ( x ) . (1.7) \rho \left( {u,{u_0},\bm{x}} \right) = uI_1^x\left( {\bm{x} + {u_0}} \right) + {I_1}\left( {\bm{x + }{u_0}} \right) - {u_0}I_1^x\left( {\bm{x} + {u_0}} \right) - {I_0}\left( \bm{x} \right).\tag{1.7} ρ(u,u0,x)=uI1x(x+u0)+I1(x+u0)−u0I1x(x+u0)−I0(x).(1.7)
在后面的算法流程中，我们可发现 ρ ( u , u 0 , x ) \rho \left( {u,{u_0},\bm{x}} \right) ρ(u,u0,x) 只会被像素 x \bm{x} x 本身使用，所以为了公式简便，同时突出我们所要优化的光流目标，不妨直接写为 ρ ( u ) \rho \left( u \right) ρ(u)。
为了优化式 (1.6)，作者引入了一个辅助光流变量 v v v，注意 v v v 和 u u u 具有同样的意义，在这里也是表示只有水平方向的运动，并不是代表垂直方向的运动。于是式 (1.6) 可表示为以下的近似凸优化问题：
E θ = ∫ Ω { ∣ ∇ u ∣ + 1 2 θ ( u − v ) 2 + λ ∣ ρ ( v ) ∣ } d x , (1.8) {E_\theta } = \int_\Omega {\left\{ {\left| {\nabla u} \right| + \frac{1}{{2\theta }}{{\left( {u - v} \right)}^2} + \lambda \left| {\rho \left( v \right)} \right|} \right\}d\bm{x}} ,\tag{1.8} Eθ=∫Ω{∣∇u∣+2θ1(u−v)2+λ∣ρ(v)∣}dx,(1.8)
其中 θ \theta θ 为一个很小的常数，通常按经验配置，用于惩罚 u u u 和 v v v 相差较大的情况，保证式 (1.8) 能近似收敛于式 (1.6)。式 (1.8) 的优化可通过交替优化 u u u 和 v v v 来实现，即在每一次迭代中：
-
固定 v v v，优化 u u u 以实现以下最优化问题：
min u ∫ Ω { ∣ ∇ u ∣ + 1 2 θ ( u − v ) 2 } d x . (1.9) \mathop {\min }\limits_u \int_\Omega {\left\{ {\left| {\nabla u} \right| + \frac{1}{{2\theta }}{{\left( {u - v} \right)}^2}} \right\}d\bm{x}} .\tag{1.9} umin∫Ω{∣∇u∣+2θ1(u−v)2}dx.(1.9)
-
固定 u u u，优化 v v v 以实现以下最优化问题：
min v ∫ Ω { 1 2 θ ( u − v ) 2 + λ ∣ ρ ( v ) ∣ } d x . (1.10) \mathop {\min }\limits_v \int_\Omega {\left\{ {\frac{1}{{2\theta }}{{\left( {u - v} \right)}^2} + \lambda \left| {\rho \left( v \right)} \right|} \right\}d\bm{x}} .\tag{1.10} vmin∫Ω{2θ1(u−v)2+λ∣ρ(v)∣}dx.(1.10)
式 (1.9) 与式 (1.10) 的优化推导涉及变分法的内容，目前已超出我们的知识范围，因此这里仅给出论文的结论，深入的理解还需要一定时间学习消化。
对于式 (1.9) 的优化，其结果可表示为：
u = v − θ d i v p , (1.11) u = v - \theta \bf{div}\bm{p},\tag{1.11} u=v−θdivp,(1.11)
其中 p \bm{p} p 是一个跟随迭代过程不断更新的辅助二维向量 ( p x , p y ) T {\left( {{p_x},{p_y}} \right)^T} (px,py)T， d i v \bf{div} div 为散度算子。实际上 p \bm{p} p 的维度与图像的维数相同，其大概的意义是反映光流 u u u 在水平以及垂直方向上的变化，但不是精确对应。正如我们前面所说的， u u u 只是表示水平方向上的运动，并不涉及垂直方向的运动；但即便是水平方向的运动，不同位置像素的水平运动速度也不一定相同，也就具有了水平光流在水平以及垂直两个方向上的变化。因此，如果我们还需要讨论垂直方向上的运动，则应该再增加一个额外的二维辅助向量，用于反映垂直方向运动速度在水平以及垂直方向上的变化。同理，如果我们需要处理三维的图像，则每一个维度上的运动都需要一个三维的辅助向量，依此类推。关于多维图像的优化我们后面再表，这里继续以只有水平运动为例进行阐述。对于 p \bm{p} p 的散度，其计算如下：
d i v p = [ ∂ ∂ x , ∂ ∂ y ] [ p x , p y ] T = ∂ p x ∂ x + ∂ p y ∂ y , (1.12) \bf{div}\bm{p} = \left[ {\frac{\partial }{{\partial x}},\frac{\partial }{{\partial y}}} \right]{\left[ {{p_x},{p_y}} \right]^T} = \frac{{\partial {p_x}}}{{\partial x}} + \frac{{\partial {p_y}}}{{\partial y}},\tag{1.12} divp=[∂x∂,∂y∂][px,py]T=∂x∂px+∂y∂py,(1.12)
即散度是一个标量。要注意我们前面的公式都省略了像素坐标， p x {p_x} px 和 p y {p_y} py 实际上都可以表示为与图像同尺寸的二维阵列，因此在求各偏导数时只需使用数值差分近似即可。辅助向量 p \bm{p} p 的更新并不直接依赖于更新后的光流 u u u，其可直接由辅助变量 v v v 更新，即：
p k + 1 = p k + τ ∇ ( d i v p k − v / θ ) 1 + τ ∣ ∇ ( d i v p k − v / θ ) ∣ , (1.13) {\bm{p}^{k + 1}} = \frac{{{\bm{p}^k} + \tau \nabla \left( {\bf{div}{\bm{p}^k} - v/\theta } \right)}}{{1 + \tau \left| {\nabla \left( {\bf{div}{\bm{p}^k} - v/\theta } \right)} \right|}},\tag{1.13} pk+1=1+τ∣∇(divpk−v/θ)∣pk+τ∇(divpk−v/θ),(1.13)
其中 k k k 是迭代的序号， τ \tau τ 是更新步长，为了收敛通常要求 τ ≤ 1 / 8 \tau \le 1/8 τ≤1/8。常取 p 0 = 0 {\bm{p}^0} = {\bf{0}} p0=0。实际上，代入式 (1.11)， p k + 1 {\bm{p}^{k + 1}} pk+1 还可由 u k + 1 {u^{k + 1}} uk+1 来更新，即：
p k + 1 = p k − τ / τ θ θ ∇ u k + 1 1 + τ / τ θ θ ∣ ∇ u k + 1 ∣ , (1.14) {\bm{p}^{k + 1}} = \frac{{{\bm{p}^k} - {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\nabla {u^{k + 1}}}}{{1 + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left| {\nabla {u^{k + 1}}} \right|}},\tag{1.14} pk+1=1+τ/τθθ∣∇uk+1∣pk−τ/τθθ∇uk+1,(1.14)
虽然采用式 (1.13) 的更新方式可以实现与式 (1.11) 的并行化，但实际测试表明式 (1.14) 的串行式更新通常能够获得更好的效果。另外，由于辅助向量 p \bm{p} p 的正负方向并没有限制，我们还可以从式 (1.11) 和 (1.14) 的反方向进行更新，于是分别有
u = v − θ d i v ( − p ) = v + θ d i v p , (1.15) u = v - \theta \bf{div}\left( { - \bm{p}} \right) = v + \theta \bf{div}\bm{p},\tag{1.15} u=v−θdiv(−p)=v+θdivp,(1.15)
以及
− p k + 1 = − p k + τ ∇ ( − d i v p k − v / θ ) 1 + τ ∣ ∇ ( − d i v p k − v / θ ) ∣ ⇒ p k + 1 = p k + τ ∇ ( d i v p k + v / θ ) 1 + τ ∣ ∇ ( d i v p k + v / θ ) ∣ = p k + τ / τ θ θ ∇ u k + 1 1 + τ / τ θ θ ∣ ∇ u k + 1 ∣ . (1.16) \begin{aligned} & - {\bm{p}^{k + 1}} = \frac{{ - {\bm{p}^k} + \tau \nabla \left( { - \bf{div}{\bm{p}^k} - v/\theta } \right)}}{{1 + \tau \left| {\nabla \left( { - \bf{div}{\bm{p}^k} - v/\theta } \right)} \right|}}\\ & \Rightarrow {\bm{p}^{k + 1}} = \frac{{{\bm{p}^k} + \tau \nabla \left( {\bf{div}{\bm{p}^k} + v/\theta } \right)}}{{1 + \tau \left| {\nabla \left( {\bf{div}{\bm{p}^k} + v/\theta } \right)} \right|}} = \frac{{{\bm{p}^k} + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\nabla {u^{k + 1}}}}{{1 + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left| {\nabla {u^{k + 1}}} \right|}}. \end{aligned} \tag{1.16} −pk+1=1+τ∣∇(−divpk−v/θ)∣−pk+τ∇(−divpk−v/θ)⇒pk+1=1+τ∣∇(divpk+v/θ)∣pk+τ∇(divpk+v/θ)=1+τ/τθθ∣∇uk+1∣pk+τ/τθθ∇uk+1.(1.16)
以上两种更新方式的效果都是一样的，在算法实现中为了统一使用加法运算可以倾向于选择后者。
对于式 (1.10) 的优化，其结果可以表示为一个阈值化操作：
v = u + { λ θ I 1 x ( x + u 0 ) i f ρ ( u ) < − λ θ ( I 1 x ( x + u 0 ) ) 2 , − λ θ I 1 x ( x + u 0 ) i f ρ ( u ) > λ θ ( I 1 x ( x + u 0 ) ) 2 , − ρ ( u ) / ρ ( u ) I 1 x ( x + u 0 ) I 1 x ( x + u 0 ) i f ∣ ρ ( u ) ∣ ≤ λ θ ( I 1 x ( x + u 0 ) ) 2 . (1.17) v = u + \left\{ {\begin{array}{c} {\lambda \theta I_1^x\left( {\bm{x} + {u_0}} \right)}&{if}&{\rho \left( u \right) < - \lambda \theta {{\left( {I_1^x\left( {\bm{x} + {u_0}} \right)} \right)}^2},}\\ { - \lambda \theta I_1^x\left( {\bm{x} + {u_0}} \right)}&{if}&{\rho \left( u \right) > \lambda \theta {{\left( {I_1^x\left( {\bm{x} + {u_0}} \right)} \right)}^2},}\\ { - {{\rho \left( u \right)} \mathord{\left/ {\vphantom {{\rho \left( u \right)} {I_1^x\left( {\bm{x} + {u_0}} \right)}}} \right. } {I_1^x\left( {\bm{x} + {u_0}} \right)}}}&{if}&{\left| {\rho \left( u \right)} \right| \le \lambda \theta {{\left( {I_1^x\left( {\bm{x} + {u_0}} \right)} \right)}^2}.} \end{array}} \right.\tag{1.17} v=u+⎩ ⎨ ⎧λθI1x(x+u0)−λθI1x(x+u0)−ρ(u)/ρ(u)I1x(x+u0)I1x(x+u0)ifififρ(u)<−λθ(I1x(x+u0))2,ρ(u)>λθ(I1x(x+u0))2,∣ρ(u)∣≤λθ(I1x(x+u0))2.(1.17)
当残差 ρ ( u ) \rho \left( u \right) ρ(u) 足够小时，辅助变量 v v v 的更新会逐渐收敛到 u u u，从而使优化问题逐步回归到最初的代价函数；而当残差较大时，辅助变量 v v v 会偏离 u u u 一个固定的步长（注意每个像素的 u 0 {u_0} u0 是一个常量），以此搜索更好的结果。因此基于以上的阈值化操作，我们可以使 ρ ( u ) \rho \left( u \right) ρ(u) 逐渐消弭，从而保证光流计算的准确性。而从式 (1.17) 也可以看出，该阈值化操作是一个逐像素运算，非常适合并行化实现；而且我们可以在迭代开始之前就把 I 1 x ( x + u 0 ) I_1^x\left( {\bm{x} + {u_0}} \right) I1x(x+u0) 求出来，从而节省很大的计算量。
1.2 二维图像TV-L1算法流程
基于以上只有水平运动情况下的 TV-L1 算法讨论，我们可以给出二维图像在一般运动下水平与垂直二维度光流的计算流程。这套流程也可十分容易地扩展到更高维度的图像光流计算中。
-
给定二维输入图像 I 0 ( x ) {I_0}\left( \bm{x} \right) I0(x) 和 I 1 ( x ) {I_1}\left( \bm{x} \right) I1(x)，图像坐标为 x = ( x , y ) T \bm{x} = {\left( {x,y} \right)^T} x=(x,y)T，记每个像素的光流为 u ( x ) = ( u x ( x ) , u y ( x ) ) T \bm{u}\left( \bm{x} \right) = {\left( {{u_x}\left( \bm{x} \right),{u_y}\left( \bm{x} \right)} \right)^T} u(x)=(ux(x),uy(x))T。分配一个与图像同尺寸的二维阵列，用于存储水平光流 u x ( x ) {u_x}\left( \bm{x} \right) ux(x)，并由一个已知的初始水平光流二维阵列 u x , 0 ( x ) {u_{x,0}}\left( \bm{x} \right) ux,0(x) 进行初始化；同理还需另一个二维阵列存储垂直光流 u y ( x ) {u_y}\left( \bm{x} \right) uy(x)，并由一个已知的初始垂直光流二维阵列 u y , 0 ( x ) {u_{y,0}}\left( \bm{x} \right) uy,0(x)进行初始化。分配两个与图像同尺寸的二维阵列，用于存储水平光流的辅助向量 p u x ( x ) = ( p u x x ( x ) , p u x y ( x ) ) T {\bm{p}_{ux}}\left( \bm{x} \right) = {\left( {p_{ux}^x\left( \bm{x} \right),p_{ux}^y\left( \bm{x} \right)} \right)^T} pux(x)=(puxx(x),puxy(x))T 的两个分量 p u x x ( x ) p_{ux}^x\left( \bm{x} \right) puxx(x) 和 p u x y ( x ) p_{ux}^y\left( \bm{x} \right) puxy(x)；同理，还需分配两个与图像同尺寸的二维阵列，用于存储垂直光流的辅助向量 p u y ( x ) {\bm{p}_{uy}}\left( \bm{x} \right) puy(x) 的两个分量 p u y x ( x ) p_{uy}^x\left( \bm{x} \right) puyx(x) 和 p u y y ( x ) p_{uy}^y\left( \bm{x} \right) puyy(x)；这四个二维阵列都需要初始化为 0 \bf{0} 0。设置参数 λ , θ , τ \lambda ,\theta ,\tau λ,θ,τ 的值。设置一个迭代停止阈值 ε \varepsilon ε，当光流的更新幅度低于此阈值时迭代结束。设置一个标量变量 E d i f f {E_{diff}} Ediff，用于统计每次迭代光流更新的幅度。
-
分配两个与图像同尺寸的二维阵列，分别用于存储图像 I 1 ( x ) {I_1}\left( \bm{x} \right) I1(x) 的水平梯度 I 1 x ( x ) I_1^x\left( \bm{x} \right) I1x(x) 与垂直梯度 I 1 y ( x ) I_1^y\left( \bm{x} \right) I1y(x)，其中梯度通过中心差分近似，即
I 1 x ( x ) = 1 2 ( I 1 ( x + 1 , y ) − I 1 ( x − 1 , y ) ) , I 1 y ( x ) = 1 2 ( I 1 ( x , y + 1 ) − I 1 ( x , y − 1 ) ) . (1.18) \begin{array}{l} I_1^x\left( \bm{x} \right) = \frac{1}{2}\left( {{I_1}\left( {x + 1,y} \right) - {I_1}\left( {x - 1,y} \right)} \right),\\ I_1^y\left( \bm{x} \right) = \frac{1}{2}\left( {{I_1}\left( {x,y + 1} \right) - {I_1}\left( {x,y - 1} \right)} \right). \end{array} \tag{1.18} I1x(x)=21(I1(x+1,y)−I1(x−1,y)),I1y(x)=21(I1(x,y+1)−I1(x,y−1)).(1.18)
另外分配两个与图像同尺寸的二维阵列，基于初始水平与垂直光流，分别插值计算 I 1 x ( x + u 0 ( x ) ) I_1^x\left( {\bm{x} + {\bm{u}_0}\left( \bm{x} \right)} \right) I1x(x+u0(x)) 以及 I 1 y ( x + u 0 ( x ) ) I_1^y\left( {\bm{x} + {\bm{u}_0}\left( \bm{x} \right)} \right) I1y(x+u0(x))，为了简洁分别记为 I 1 , w x ( x ) I_{1,w}^x\left( \bm{x} \right) I1,wx(x) 以及 I 1 , w y ( x ) I_{1,w}^y\left( \bm{x} \right) I1,wy(x)。因为目标位置 ( x + u x , 0 ( x ) , y + u y , 0 ( x ) ) \left( {x + {u_{x,0}}\left( \bm{x} \right),y + {u_{y,0}}\left( \bm{x} \right)} \right) (x+ux,0(x),y+uy,0(x)) 可能为非整数，此时的插值可通过双线性或双三次插值来实现。同理，分配一个与图像尺寸相同的二维阵列，基于初始水平与垂直光流，插值计算 I 1 ( x + u 0 ( x ) ) {I_1}\left( {\bm{x} + {\bm{u}_0}\left( \bm{x} \right)} \right) I1(x+u0(x))，记为 I 1 , w ( x ) {I_{1,w}}\left( \bm{x} \right) I1,w(x)。注意以上的插值只需要在迭代开始前进行一次，后续迭代过程中无需更新。
-
分配两个与图像同尺寸的二维阵列，分别用于存储水平与垂直光流辅助变量 v x ( x ) {v_x}\left( \bm{x} \right) vx(x) 和 v y ( x ) {v_y}\left( \bm{x} \right) vy(x)，两者无需初始化。对于二维以及更高维的图像，式 (1.17) 的阈值化操作可改写为：
v = u + { λ θ ∇ I 1 ( x + u 0 ) i f ρ ( u ) < − λ θ ∣ ∇ I 1 ( x + u 0 ) ∣ 2 , − λ θ ∇ I 1 ( x + u 0 ) i f ρ ( u ) > λ θ ∣ ∇ I 1 ( x + u 0 ) ∣ 2 , − ρ ( u ) ∇ I 1 ( x + u 0 ) ∣ ∇ I 1 ( x + u 0 ) ∣ 2 i f ∣ ρ ( u ) ∣ ≤ λ θ ∣ ∇ I 1 ( x + u 0 ) ∣ 2 , (1.19) \bm{v} = \bm{u} + \left\{ {\begin{array}{c} {\lambda \theta \nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)}&{if}&{\rho \left( \bm{u} \right) < - \lambda \theta {{\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|}^2},}\\ { - \lambda \theta \nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)}&{if}&{\rho \left( \bm{u} \right) > \lambda \theta {{\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|}^2},}\\ { - \frac{{\rho \left( u \right)\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)}}{{{{\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|}^2}}}}&{if}&{\left| {\rho \left( \bm{u} \right)} \right| \le \lambda \theta {{\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|}^2},} \end{array}} \right.\tag{1.19} v=u+⎩ ⎨ ⎧λθ∇I1(x+u0)−λθ∇I1(x+u0)−∣∇I1(x+u0)∣2ρ(u)∇I1(x+u0)ifififρ(u)<−λθ∣∇I1(x+u0)∣2,ρ(u)>λθ∣∇I1(x+u0)∣2,∣ρ(u)∣≤λθ∣∇I1(x+u0)∣2,(1.19)
其中
∇ I 1 ( x + u 0 ) = [ I 1 , w x ( x ) , I 1 , w y ( x ) ] T , ∣ ∇ I 1 ( x + u 0 ) ∣ 2 = ( I 1 , w x ( x ) ) 2 + ( I 1 , w y ( x ) ) 2 . (1.20) \begin{array}{l} \nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right) = {\left[ {I_{1,w}^x\left( \bm{x} \right),I_{1,w}^y\left( \bm{x} \right)} \right]^T},\\ {\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|^2} = \sqrt {{{\left( {I_{1,w}^x\left( \bm{x} \right)} \right)}^2} + {{\left( {I_{1,w}^y\left( \bm{x} \right)} \right)}^2}} . \end{array}\tag{1.20} ∇I1(x+u0)=[I1,wx(x),I1,wy(x)]T,∣∇I1(x+u0)∣2=(I1,wx(x))2+(I1,wy(x))2.(1.20)
且有
ρ ( u ) = ⟨ u , ∇ I 1 ( x + u 0 ) ⟩ + I 1 ( x + u 0 ) − ⟨ u 0 , ∇ I 1 ( x + u 0 ) ⟩ − I 0 ( x ) = ⟨ u , ∇ I 1 ( x + u 0 ) ⟩ + ρ c ( x ) , (1.21) \begin{aligned} \rho \left( \bm{u} \right) &= \left\langle {\bm{u},\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right\rangle + {I_1}\left( {\bm{x + }{\bm{u}_0}} \right) - \left\langle {{\bm{u}_0},\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right\rangle - {I_0}\left( \bm{x} \right)\\ &= \left\langle {\bm{u},\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right\rangle + {\rho _c}\left( \bm{x} \right), \end{aligned}\tag{1.21} ρ(u)=⟨u,∇I1(x+u0)⟩+I1(x+u0)−⟨u0,∇I1(x+u0)⟩−I0(x)=⟨u,∇I1(x+u0)⟩+ρc(x),(1.21)
其中 ⟨ a , b ⟩ \left\langle {\bm{a},\bm{b}} \right\rangle ⟨a,b⟩ 代表向量的内积。因为 ∣ ∇ I 1 ( x + u 0 ) ∣ 2 {\left| {\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right|^2} ∣∇I1(x+u0)∣2 不随迭代而改变，如果想节省后续计算量可在迭代开始前分配存储空间并事先计算存储。同样，因为式 (1.21) 中可以提取出不跟随迭代过程改变的常量项 ρ c ( x ) {\rho _c}\left( \bm{x} \right) ρc(x)，可分配一个与图像同尺寸的二维阵列，并在迭代开始前进行初始化存储，即
ρ c ( x ) = I 1 ( x + u 0 ) − ⟨ u 0 , ∇ I 1 ( x + u 0 ) ⟩ − I 0 ( x ) = I 1 , w ( x ) − u x , 0 ( x ) I 1 , w x ( x ) − u y , 0 ( x ) I 1 , w y ( x ) − I 0 ( x ) . (1.22) \begin{array}{l} {\rho _c}\left( \bm{x} \right) = {I_1}\left( {\bm{x} + {\bm{u}_0}} \right) - \left\langle {{\bm{u}_0},\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right\rangle - {I_0}\left( \bm{x} \right)\\ = {I_{1,w}}\left( \bm{x} \right) - {u_{x,0}}\left( \bm{x} \right)I_{1,w}^x\left( \bm{x} \right) - {u_{y,0}}\left( \bm{x} \right)I_{1,w}^y\left( \bm{x} \right) - {I_0}\left( \bm{x} \right). \end{array}\tag{1.22} ρc(x)=I1(x+u0)−⟨u0,∇I1(x+u0)⟩−I0(x)=I1,w(x)−ux,0(x)I1,wx(x)−uy,0(x)I1,wy(x)−I0(x).(1.22)
完成以上初始化工作后，我们将进入循环迭代的流程。
-
根据当前迭代的光流 u x ( x ) {u_x}\left( \bm{x} \right) ux(x) 和 u y ( x ) {u_y}\left( \bm{x} \right) uy(x)，计算残差项
ρ ( u ) = ⟨ u , ∇ I 1 ( x + u 0 ) ⟩ + ρ c ( x ) = u x ( x ) I 1 , w x ( x ) + u y ( x ) I 1 , w y ( x ) + ρ c ( x ) . (1.23) \begin{aligned} \rho \left( \bm{u} \right) &= \left\langle {\bm{u},\nabla {I_1}\left( {\bm{x} + {\bm{u}_0}} \right)} \right\rangle + {\rho _c}\left( \bm{x} \right)\\ &= {u_x}\left( \bm{x} \right)I_{1,w}^x\left( \bm{x} \right) + {u_y}\left( \bm{x} \right)I_{1,w}^y\left( \bm{x} \right) + {\rho _c}\left( \bm{x} \right). \end{aligned}\tag{1.23} ρ(u)=⟨u,∇I1(x+u0)⟩+ρc(x)=ux(x)I1,wx(x)+uy(x)I1,wy(x)+ρc(x).(1.23)
根据式 (1.19)，更新 v x ( x ) {v_x}\left( \bm{x} \right) vx(x) 与 v y ( x ) {v_y}\left( \bm{x} \right) vy(x) 的值。
-
初始化 E d i f f = 0 {E_{diff}} = 0 Ediff=0。对于水平光流 u x ( x ) {u_x}\left( \bm{x} \right) ux(x)，首先计算其辅助向量 p u x ( x ) {\bm{p}_{ux}}\left( \bm{x} \right) pux(x) 的散度，其中在计算偏导数时常只采用后向差分近似，而不是中心差分，即有
d i v p u x ( x ) = ∂ p u x x ( x ) ∂ x + ∂ p u x y ( x ) ∂ y ≈ p u x x ( x , y ) − p u x x ( x − 1 , y ) + p u x y ( x , y ) − p u x y ( x , y − 1 ) . (1.24) \begin{aligned} & \bf{div}{\bm{p}_{ux}}\left( \bm{x} \right) = \frac{{\partial p_{ux}^x\left( \bm{x} \right)}}{{\partial x}} + \frac{{\partial p_{ux}^y\left( \bm{x} \right)}}{{\partial y}}\\ & \approx p_{ux}^x\left( {x,y} \right) - p_{ux}^x\left( {x - 1,y} \right) + p_{ux}^y\left( {x,y} \right) - p_{ux}^y\left( {x,y - 1} \right). \end{aligned}\tag{1.24} divpux(x)=∂x∂puxx(x)+∂y∂puxy(x)≈puxx(x,y)−puxx(x−1,y)+puxy(x,y)−puxy(x,y−1).(1.24)
统计水平光流更新幅度以及更新水平光流可得
E d i f f = E d i f f + ( v x ( x ) + θ d i v p u x ( x ) − u x ( x ) ) 2 , u x ( x ) = v x ( x ) + θ d i v p u x ( x ) . (1.25) \begin{gathered} {E_{diff}} = {E_{diff}} + {\left( {{v_x}\left( \bm{x} \right) + \theta \bf{div}{\bm{p}_{ux}}\left( \bm{x} \right) - {u_x}\left( \bm{x} \right)} \right)^2},\\ {u_x}\left( \bm{x} \right) = {v_x}\left( \bm{x} \right) + \theta \bf{div}{\bm{p}_{ux}}\left( \bm{x} \right). \end{gathered}\tag{1.25} Ediff=Ediff+(vx(x)+θdivpux(x)−ux(x))2,ux(x)=vx(x)+θdivpux(x).(1.25)
同理对于垂直光流，有
d i v p u y ( x ) ≈ p u y x ( x , y ) − p u y x ( x − 1 , y ) + p u y y ( x , y ) − p u y y ( x , y − 1 ) . (1.26) \bf{div}{\bm{p}_{uy}}\left( \bm{x} \right) \approx p_{uy}^x\left( {x,y} \right) - p_{uy}^x\left( {x - 1,y} \right) + p_{uy}^y\left( {x,y} \right) - p_{uy}^y\left( {x,y - 1} \right).\tag{1.26} divpuy(x)≈puyx(x,y)−puyx(x−1,y)+puyy(x,y)−puyy(x,y−1).(1.26)
以及
E d i f f = E d i f f + ( v y ( x ) + θ d i v p u y ( x ) − u y ( x ) ) 2 , u y ( x ) = v y ( x ) + θ d i v p u y ( x ) . (1.27) \begin{gathered} {E_{diff}} = {E_{diff}} + {\left( {{v_y}\left( \bm{x} \right) + \theta \bf{div}{\bm{p}_{uy}}\left( \bm{x} \right) - {u_y}\left( \bm{x} \right)} \right)^2},\\ {u_y}\left( \bm{x} \right) = {v_y}\left( \bm{x} \right) + \theta \bf{div}{\bm{p}_{uy}}\left( \bm{x} \right). \end{gathered}\tag{1.27} Ediff=Ediff+(vy(x)+θdivpuy(x)−uy(x))2,uy(x)=vy(x)+θdivpuy(x).(1.27)
-
基于以上已更新的水平光流 u x ( x ) {u_x}\left( \bm{x} \right) ux(x)，根据式 (1.16) 更新辅助向量 p u x ( x ) {\bm{p}_{ux}}\left( \bm{x} \right) pux(x)，这里 u x ( x ) {u_x}\left( \bm{x} \right) ux(x) 的梯度我们采用前向差分近似，从而与前面求散度时所用的后向差分相互弥补误差，即
p u x x ( x ) = p u x x ( x ) + τ / τ θ θ ( u x ( x + 1 , y ) − u x ( x , y ) ) 1 + τ / τ θ θ ∣ ∇ u x ( x ) ∣ , p u x y ( x ) = p u x y ( x ) + τ / τ θ θ ( u x ( x , y + 1 ) − u x ( x , y ) ) 1 + τ / τ θ θ ∣ ∇ u x ( x ) ∣ , (1.28) \begin{gathered} p_{ux}^x\left( \bm{x} \right) = \frac{{p_{ux}^x\left( \bm{x} \right) + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left( {{u_x}\left( {x + 1,y} \right) - {u_x}\left( {x,y} \right)} \right)}}{{1 + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left| {\nabla {u_x}\left( \bm{x} \right)} \right|}},\\ p_{ux}^y\left( \bm{x} \right) = \frac{{p_{ux}^y\left( \bm{x} \right) + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left( {{u_x}\left( {x,y + 1} \right) - {u_x}\left( {x,y} \right)} \right)}}{{1 + {\tau \mathord{\left/ {\vphantom {\tau \theta }} \right. } \theta }\left| {\nabla {u_x}\left( \bm{x} \right)} \right|}}, \end{gathered}\tag{1.28} puxx(x)=1+τ/τθθ∣∇ux(x)∣puxx(x)+τ/τθθ(ux(x+1,y)−ux(x,y)),puxy(x)=1+τ/τθθ∣∇ux(x)∣puxy(x)+τ/τθθ(ux(x,y+1)−ux(x,y)),(1.28)
其中
∣ ∇ u x ( x ) ∣ = ( ∂ u x ( x ) ∂ x ) 2 + ( ∂ u x ( x ) ∂ y ) 2 ≈ ( u x ( x + 1 , y ) − u x ( x , y ) ) 2 + ( u x ( x , y + 1 ) − u x ( x , y ) ) 2 . (1.29) \begin{aligned} & \left| {\nabla {u_x}\left( \bm{x} \right)} \right| = \sqrt {{{\left( {\frac{{\partial {u_x}\left( \bm{x} \right)}}{{\partial x}}} \right)}^2} + {{\left( {\frac{{\partial {u_x}\left( \bm{x} \right)}}{{\partial y}}} \right)}^2}} \\ & \approx \sqrt {{{\left( {{u_x}\left( {x + 1,y} \right) - {u_x}\left( {x,y} \right)} \right)}^2} + {{\left( {{u_x}\left( {x,y + 1} \right) - {u_x}\left( {x,y} \right)} \right)}^2}} . \end{aligned}\tag{1.29} ∣∇ux(x)∣=(∂x∂ux(x))2+(∂y∂ux(x))2≈(ux(x+1,y)−ux(x,y))2+(ux(x,y+1)−ux(x,y))2.(1.29)
对于垂直光流的辅助向量两个分量 p u y x ( x ) p_{uy}^x\left( \bm{x} \right) puyx(x) 和 p u y y ( x ) p_{uy}^y\left( \bm{x} \right) puyy(x) 的更新亦同理，只需把 u x {u_x} ux 替代为 u y {u_y} uy 即可，这里不再赘述。
-
检查光流的平均更新幅度，判断迭代是否以及停滞。假设图像的总像素个数为 N p i x {N_{pix}} Npix，如果满足
E d i f f / E d i f f N p i x N p i x ≤ ε , (1.30) \sqrt {{{{E_{diff}}} \mathord{\left/ {\vphantom {{{E_{diff}}} {{N_{pix}}}}} \right. } {{N_{pix}}}}} \le \varepsilon ,\tag{1.30} Ediff/EdiffNpixNpix≤ε,(1.30)
则终止迭代程序，将此次更新所得光流 u x ( x ) {u_x}\left( \bm{x} \right) ux(x) 和 u y ( x ) {u_y}\left( \bm{x} \right) uy(x) 作为输出；否则返回第 (4) 步继续以上流程。我们也可设置一个最大迭代次数来避免收敛过慢的情况。
以上便是二维图像的 TV-L1 光流计算流程。由于算法保留了初始光流 u 0 ( x ) {\bm{u}_0}\left( \bm{x} \right) u0(x) 这一可配置项，我们事实上可以设计一个更大的迭代过程。例如，在一开始我们可以将初始光流 u 0 ( x ) {\bm{u}_0}\left( \bm{x} \right) u0(x) 全部设置为 0 {\bf{0}} 0，这是十分符合直觉的；基于以上的算法流程，我们可以得到一个可能不那么准确的光流输出 u ( x ) \bm{u}\left( \bm{x} \right) u(x)，这是因为大多数光流法都不太擅长处理特别大幅度的运动；如果以所得的光流 u ( x ) \bm{u}\left( \bm{x} \right) u(x) 作为初始光流再执行一次光流算法，即便该初始光流中存在一些不那么准确的部分，但在一个大的邻域内依然我们仍能找到很多正确的参照对象，TV-L1 算法也具有误差修正的能力，于是理论上我们能够获得更加准确的光流输出。除此以外，我们也可以在多尺度上获取更准确的初始光流，因为通过缩小图像，我们一方面可以使大幅度的运动变为小幅度的运动，另一方面还可以降低原始尺度图像中可能存在的噪声影响。因此，多尺度是大多数光流法的有效改善手段。
2 多尺度的TV-L1算法流程
这里我们将给出 TV-L1 算法与多尺度方法结合的算法流程，实际上这种多尺度方法可以轻松移植到其他光流算法中。
-
给定输入图像 I 0 , s = 0 {I_{0,s = 0}} I0,s=0 和 I 1 , s = 0 {I_{1,s = 0}} I1,s=0，其中 s = 0 s = 0 s=0 表示图像在第 0 尺度下。设置图像缩小倍率 f d n < 1 {f_{dn}} < 1 fdn<1，以及最大缩小次数 N d n {N_{dn}} Ndn。设置在每一尺度下 TV-L1 算法的优化次数 N w a r p {N_{warp}} Nwarp。TV-L1 算法所需的参数这里不作为重点，可参考第 1 节内容。
-
分别对 I 0 , s = 0 {I_{0,s = 0}} I0,s=0 和 I 1 , s = 0 {I_{1,s = 0}} I1,s=0 根据所设置的缩小倍率 f d n {f_{dn}} fdn 进行缩小，得到第 1 尺度下的图像 I 0 , s = 1 {I_{0,s = 1}} I0,s=1 和 I 1 , s = 1 {I_{1,s = 1}} I1,s=1。具体缩小插值的方法不限，如双线性或双三次插值方法都可。继续对 I 0 , s = 1 {I_{0,s = 1}} I0,s=1 和 I 1 , s = 1 {I_{1,s = 1}} I1,s=1 进行 f d n {f_{dn}} fdn 倍率的缩小，以此类推，直至获得第 N d n {N_{dn}} Ndn 尺度下的图像 I 0 , s = N d n {I_{0,s = {N_{dn}}}} I0,s=Ndn 和 I 1 , s = N d n {I_{1,s = {N_{dn}}}} I1,s=Ndn。
-
在第 N d n {N_{dn}} Ndn 尺度下，将初始光流 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn 初始化为 0 {\bf{0}} 0，注意初始光流的水平光流和垂直光流分别由与第 N d n {N_{dn}} Ndn 尺度下图像同尺寸的二维阵列存储。基于 I 0 , s = N d n {I_{0,s = {N_{dn}}}} I0,s=Ndn 和 I 1 , s = N d n {I_{1,s = {N_{dn}}}} I1,s=Ndn 以及初始光流 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn，根据第 1 节的 TV-L1 算法流程，计算得到输出光流 u s = N d n {\bm{u}_{s = {N_{dn}}}} us=Ndn，将 u s = N d n {\bm{u}_{s = {N_{dn}}}} us=Ndn 复制到 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn。在第 N d n {N_{dn}} Ndn 尺度下基于 I 0 , s = N d n {I_{0,s = {N_{dn}}}} I0,s=Ndn 和 I 1 , s = N d n {I_{1,s = {N_{dn}}}} I1,s=Ndn 以及更新的初始光流 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn 重新执行 TV-L1 算法流程，同样将输出光流 u s = N d n {\bm{u}_{s = {N_{dn}}}} us=Ndn 复制到 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn；重复此操作，直至在第 N d n {N_{dn}} Ndn 尺度下总共执行了 N w a r p {N_{warp}} Nwarp 次 TV-L1 算法流程。
-
在第 N d n − 1 {N_{dn}} - 1 Ndn−1 尺度下，首先将第 N d n {N_{dn}} Ndn 尺度下最后更新所得的初始光流 u 0 , s = N d n {\bm{u}_{0,s = {N_{dn}}}} u0,s=Ndn 的每个值除以缩小倍率 f d n {f_{dn}} fdn，即图像尺度缩放时其光流大小也会等倍率缩放；然后分别将其两个分量即水平光流图和垂直光流图通过双线性或双三次等插值方法插值放大至第 N d n − 1 {N_{dn}} - 1 Ndn−1 尺度下的图像尺寸，得到第 N d n − 1 {N_{dn}} - 1 Ndn−1 尺度下的初始光流 u 0 , s = N d n − 1 {\bm{u}_{0,s = {N_{dn}} - 1}} u0,s=Ndn−1。基于 I 0 , s = N d n − 1 {I_{0,s = {N_{dn}} - 1}} I0,s=Ndn−1 和 I 1 , s = N d n − 1 {I_{1,s = {N_{dn}} - 1}} I1,s=Ndn−1 以及初始光流 u 0 , s = N d n − 1 {\bm{u}_{0,s = {N_{dn}} - 1}} u0,s=Ndn−1，执行和第 (3) 步一样的（除了最初的初始光流不一样外）总共 N w a r p {N_{warp}} Nwarp 次的 TV-L1 循环流程。
-
对于第 N d n − 2 {N_{dn}} - 2 Ndn−2 尺度，取第 N d n − 1 {N_{dn}} - 1 Ndn−1 尺度下最后更新所得的初始光流 u 0 , s = N d n − 1 {\bm{u}_{0,s = {N_{dn}} - 1}} u0,s=Ndn−1 ，执行和第 (4) 步相似的流程。以此类推，直至执行完毕第 0 尺度即原图尺度下的算法流程。将第 0 尺度下最后更新所得初始光流 u 0 , s = 0 {\bm{u}_{0,s = 0}} u0,s=0 作为最终的光流输出，算法结束。
扫一扫
1440
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
