光流不仅包含了被观察物体的运动信息，还携带有被观察物体丰富的3维结构等信息，被广泛地应用于目标识别、分割、跟踪、目标3D结构和形状恢复、视频检测、图像配准、运动估计等计算机视觉和图像处理领域。

在Horn等人^[1]和Lucas等人^[2]分别提出了全局和局部变分光流计算的两种基本方法后，各种针对变分光流计算的新理论和新方法不断涌现，在计算精度及鲁棒性等方面得到了较大提高。一般光流计算中的变分光流能量函数由数据项和平滑项组成，其中数据项主要包括各种基于图像数据先验知识的守恒假设，是决定光流计算的主要因素。Horn算法中采用了灰度守恒假设，并得到了相应的光流基本约束方程。Zimmer等人^[3]在Horn算法的模型中增加了惩罚函数形式的数据项和平滑项，且数据项增加了梯度守恒假设。Papenberg等人^[4]证明了光流计算基本约束方程与梯度守恒假设相结合计算精度较高、计算复杂度相对较低、总体性能最优。能量函数中平滑项主要由各种平滑策略构成，决定了在计算时光流扩散的方向和程度。常用的平滑项有各向同性扩散的非线性平滑策略^[5]、基于图像结构张量的各向异性扩散策略^[6]、将各项同性扩散与各向异性扩散策略相结合的策略^[7]和马尔可夫随机场(MRF)策略^[8]等。

马尔可夫随机场模型(MRF)假设图像的邻域空间分布具有局部性、一致性，使得MRF光流计算模型具有优良的抗噪能力。早期的MRF能量函数优化算法(ICM算法、模拟退火算法和动态规划等)速度慢，容易陷入局部最优。随着图割(GC)^[9]、置信传播(BP)^[10]等新的强有力的优化算法的出现，MRF能量函数的优化速度有了质的提升。

虽然图像序列变分光流计算技术在运动检测方面已取得了很大的进展，但由于变分光流计算的时间消耗通常较大，限制了其在实时运动检测中的应用。采用并行计算是提高光流计算效率的一个解决途径。2003年以来，GPU的运算能力就远远超过了CPU。在存储器带宽方面，GPU也优于CPU。由于在并行计算方面的优越性，除了进行传统的图形图像显示外，GPU还被用来进行诸如图像分割^[11]、配准^[12]、深度学习^[13]等通用计算(GPGPU)。为了解决GPU可编程和存储问题，NVIDIA公司于2007年正式发布一种新的统一计算架构CUDA。经过近年来的发展，GPU性能显著提高，CUDA的通用计算能力也在显著提高，而且功能也在不断完善。Cuomo等人^[14]在CUDA平台上实现了一个GPU加速的光流场计算方法，该方法将光流变分计算变为两个线性方程的迭代解的问题，其计算速度仍然不够快。Hossen等人^[15]在CUDA平台上实现了一种GPU加速的Lucas光流场计算方法。Grauer-Gray等人^[16]在CUDA平台上提出了基于马尔可夫能量函数的变分光流并行计算方法，但是该方法只是简单将串行方法并行化，没有针对光流计算特点进行优化，导致单线程计算负荷重，计算效率依然不够高。

本文提出了一种优化的基于马尔可夫模型的光流并行计算方法，该优化方法采用了一种4层的3维网络结构进行并行BP计算，每层负责4邻域MRF中一个方向的BP信息计算及传递，同时为每层中的每个像素点分配多个线程并行计算，大大降低单线程计算负荷，大幅度提高计算效率。

一般光流计算中的变分光流能量函数由数据项和平滑项组成，其能量函数如下

$ E\left( {u,v} \right) = {E_{{\rm{smooth}}}}\left( {u,v} \right) + \lambda {E_{{\rm{data}}}}\left( {u,v} \right) $

(1)

式中，$u$, $v$分别代表两个方向上的光流，$\lambda $为比例系数。式(1)可以通过变分方式求解，常用的方法是最速下降法，速度快，但是容易陷入局部最小，对噪声敏感。把图像中$p$点的光流值离散为取值范围为1到$L$的光流标记$u_p$和$v_p$，当$L$为偶数时，其对应的光流取值范围为$\left({ - \left({\alpha L/2} \right)\sim \left({\left({\alpha L/2} \right)} \right. - 1} \right)$，其中$\alpha $为采样间隔，本文取$L$=16。令

$ {E_{{\rm{smooth}}}}\left( {u,v} \right) = \sum\limits_{q \in {\mathit{\boldsymbol{\mathcal{N}}}_p}} {V\left( {{u_p} - {u_q},{v_p} - {v_q}} \right)} $

(2)

$ {E_{{\rm{data}}}}\left( {u,v} \right) = \sum\limits_{p \in {\cal P}} {{D_p}\left( {{u_p},{v_p}} \right)} $

(3)

则变分光流能量模型转化为马尔可夫能量模型

$ \begin{array}{*{20}{c}} {E\left( {u,v} \right) = }\\ {\sum\limits_{q \in {\mathit{\boldsymbol{\mathcal{N}}}_p}} {V\left( {{u_p} - {u_q},{v_p} - {v_q}} \right)} + \sum\limits_{p \in \mathit{\boldsymbol{\mathcal{P}}}} {{D_p}\left( {{u_p},{v_p}} \right)} } \end{array} $

(4)

式中，等式右边第1项对应平滑项，$V$为代价函数，表示把光流标记($u_p$, $v_p$)和($u_q$, $v_q$)同时赋给两个相邻节点$p$和$q$的代价，$\mathit{\boldsymbol{\mathcal{N}}}p$表示像素点$p$的邻域，一般2维MRF中取上下左右4邻域，$\mathit{\boldsymbol{\mathcal{P}}}$代表图像域。

第2项对应数据项，${D_p}\left({{u_p}, {v_p}} \right)$表示将标记($u_p$, $v_p$)赋给$p$点的代价。光流场应该是连续而稠密的矢量集合，但是在物体边界处可能会出现因计算不准确而造成的误差较大的情况，采用非平方惩罚函数平滑运动矢量，可以大大减弱边界不连续带来的不利影响。本文在数据项中采用$\psi \left({{D_p}} \right) = \sqrt {D_p^2 + {L^2}} $作为非平方惩罚函数，其中$L$为采样范围。则本文采用的数据项可表示为

$ {E_{{\rm{data}}}}\left( {u,v} \right) = \sum\limits_{p \in {\cal P}} {\psi \left( {{D_p}\left( {{u_p},{v_p}} \right)} \right)} $

(5)

在计算代价函数${D_p}\left({{u_p}, {v_p}} \right)$时，Felzenszwalb等人^[8]基于灰度守恒采用了方程

$ \begin{array}{*{20}{c}} {{D_p}\left( {{u_p},{v_p}} \right) = }\\ {\left\| {{I_1}\left( {x + \alpha {u_p},y + \alpha {v_p}} \right) - {I_2}\left( {x,y} \right)} \right\|} \end{array} $

(6)

采用该数据项的不足之处在于坐标($x + \alpha {u_p}, y + \alpha {v_p}$)多数情况下不是整数，这样需要采用插值方法计算式(6)，如果插值不合理会导致光流计算错误。为此本文在式(7)中根据Horn灰度守恒假设^[1]推导出的光流基本约束方程，增加了一个新的数据项，该数据项不涉及数据插值，提高了光流计算精度，即

$ \begin{array}{*{20}{c}} {{D_p}\left( {{u_p},{v_p}} \right) = \left\| {{I_1}\left( {x + \alpha {u_p},y + \alpha {v_p}} \right) - {I_2}\left( {x,y} \right)} \right\| + }\\ {\gamma \left\| \begin{array}{l} I_1^x\left( {x,y} \right)\alpha {u_p} + I_1^y\left( {x,y} \right)\alpha {v_p} + \\ {I_2}\left( {x,y} \right) - {I_1}\left( {x,y} \right) \end{array} \right\|} \end{array} $

(7)

式中，$I_i^x$和$I_i^y$为第$i$幅输入图像在$X$和$Y$方向上的梯度，$\gamma $为比例系数。

置信传播(BP)算法的主要思想是：对于马尔可夫随机场中的每一个节点，通过信息传播，把该节点的能量状态传递给相邻的节点，从而影响相邻节点的能量状态，经过一定次数的迭代，每个节点的能量将收敛于一个稳态。采用BP算法中的Max-Product方法，节点之间传送的信息定义为

$ \begin{array}{*{20}{c}} {m_{pq}^t\left( {{f_q}} \right) = \mathop {\min }\limits_{{f_p}} \left( {V\left( {{f_p},{f_q}} \right) + {D_p}\left( {{f_p}} \right) + } \right.}\\ {\left. {\sum\limits_{s \in \mathit{\boldsymbol{\mathcal{N}}}\left( p \right)\backslash q} {m_{sp}^{t - 1}\left( {{f_p}} \right)} } \right)} \end{array} $

(8)

式中，$m_{pq}^t\left({{f_q}} \right)$表示在$t$时刻节点$p$向邻域点$q$的标记$f_q$传播的信息，$s \in \mathit{\boldsymbol{\mathcal{N}}}\left(p \right)\backslash p$表示$p$点排除$q$点的邻域，其传递过程如图 1所示。在信息经过$T$遍传播后，得到节点$q$的置信度为

$ {b_q}\left( {{f_q}} \right) = {D_q}\left( {{f_q}} \right) + \sum\limits_{p \in \mathit{\boldsymbol{\mathcal{N}}}\left( q \right)} {m_{sp}^T\left( {{f_q}} \right)} $

(9)

对节点$q$，使置信度${b_q}\left({{f_q}} \right)$最小的$f_q$就是节点$s \in \mathit{\boldsymbol{\mathcal{N}}}\left(p \right)\backslash q$在MRF中的解。当将BP模型应用到解式(4)的光流模型时，只需将1维的$f$符号变成2维的($u$，$v$)符号即可。

BP算法中信息传递的过程是：遍历每个像素点，每个像素点分别从上下左右4个方向将式(8)中计算的信息传递给相邻像素，经过一定次数的迭代后，信息传递收敛，不再进行信息传递。具体传递时为减少计算量，采用二分法每次循环交叉更新一半图像像素点的信息。二分法将整个像素网格分为两个一样大小的集合$\mathit{\boldsymbol{A}}$和$\mathit{\boldsymbol{B}}$。图 2显示了一个向右传播的示例。当循环次数为偶数时，对应左边的图，信息由$\mathit{\boldsymbol{B}}$(黑)向右传递给$\mathit{\boldsymbol{A}}$(白)。在传递之前，$\mathit{\boldsymbol{B}}$中的信息由其上下左3个邻域的$\mathit{\boldsymbol{A}}$的信息加上自身的信息更新。当循环次数为奇数时，对应右边的图，信息由$\mathit{\boldsymbol{A}}$向右传递给$\mathit{\boldsymbol{B}}$。在传递之前，$\mathit{\boldsymbol{A}}$中的信息由其上下左3个邻域的$\mathit{\boldsymbol{B}}$的信息加上自身的信息更新。

直接采用上述BP算法计算光流需要的迭代次数较多，计算时间长。Hossen等人^[15]在其提出的并行Lucas光流计算模型中采用了金字塔算法，大大减少迭代次数。金字塔将图像数据项依次进行二分降采样，从下到上构造一个多层的金字塔(图 3)，计算时从上到下进行计算。上层金字塔计算得到的信息为下一层金字塔信息计算的初始值，重复这一过程直到最底层，这样每层经过少数迭代即能得到较为准确的结果。

在阐述了MRF光流计算基本原理及如何采用BP算法最小化MRF光流能量方程之后，下面详细说明采用BP方法计算光流的步骤:

1) 计算数据项：载入连续两帧图像根据式(7)计算数据项，这样每个像素点需进行$L^2$次计算，得到$L^2$个数据项。假设一幅图像有$N$个像素点，本步计算复杂度为${\rm{O}}\left({N{L^2}} \right)$。

2) 计算并传递信息：因为光流计算是2维的，于是将式(8)的$p$点传递给$q$点的信息改为

$ \begin{array}{*{20}{c}} {m_{pq}^t\left( {{u_q},{v_q}} \right) = \mathop {\min }\limits_{{u_p},{v_p}} \left( {{{\left( {{u_q} - {u_p}} \right)}^2} + } \right.}\\ {\left. {{{\left( {{v_q} - {v_p}} \right)}^2} + h\left( {{u_p},{v_p}} \right)} \right)} \end{array} $

(10)

$ \begin{array}{*{20}{c}} {{u_p},{u_q},{v_q},{v_q} = - \left( {L/2} \right), \cdots ,\left( {L/2} \right) - 1}\\ {h\left( {{u_p},{v_p}} \right) = }\\ {\sum\limits_{s \in \mathit{\boldsymbol{\mathcal{N}}}\left( p \right)\backslash q} {\left( {m_{sp}^{t - 1}\left( {{u_p},{v_p}} \right) + \psi \left( {{D_p}\left( {{u_p},{v_p}} \right)} \right)} \right)} } \end{array} $

(11)

如采用常规方法求解式(10)，每个像素点的计算复杂度为${\rm{O}}\left({{L^4}} \right)$，最快方法是降维法^[17]，其计算复杂度只有${\rm{O}}\left({4{L^2}} \right)$。其计算原理为：降维法将式(10)的2维模型转化为下式所示两次1维模型来降维计算

$ \begin{array}{l} m_{pq}^t\left( {{u_q},{v_q}} \right) = \mathop {\min }\limits_{{u_p}} \left( {{{\left( {{u_q} - {u_p}} \right)}^2} + \mathop {\min }\limits_{{v_p}} \left( {{{\left( {{v_q} - {v_p}} \right)}^2} + } \right.} \right.\\ \left. {\left. {h\left( {{u_p},{v_p}} \right)} \right)} \right) = \mathop {\min }\limits_{{u_p}} \left( {{{\left( {{u_q} - {u_p}} \right)}^2} + m_{h\left| {{u_p}} \right.}^t\left( {{v_q}} \right)} \right) \end{array} $

(12)

$ m_{h\left| {{u_p}} \right.}^t\left( {{v_q}} \right) = \mathop {\min }\limits_{{v_p}} \left( {{{\left( {{v_q} - {v_p}} \right)}^2} + h\left( {{u_p},{v_p}} \right)} \right) $

(13)

如果还是采用常规计算，式(13)中对于单个$v_q$的计算复杂度为${\rm{O}}\left({{L^2}} \right)$，由于$v_q$取值有$L$种可能，故总体计算复杂度为${\rm{O}}\left({{L^3}} \right)$，则式(12)的计算复杂度为${\rm{O}}\left({2{L^3}} \right)$。降维法中采用了如图 4所示的抛物线下界法来进一步优化计算式(13)。

式(13)中的每一个平方项对应图 4中的一个抛物线，当$L$=5时，共有5个抛物线。可以看出抛物线1、3、5构成这些抛物线的下界，而这个下界在每个点上的长度(加粗线段)就是该点的$m_{h\left| {{u_p}} \right.}^t\left({{v_q}} \right), {v_q} = 1, \cdots, 5$。抛物线的下界可以由求抛物线的交点得到。假设$r$、$q$、$z$为3个由左向右的坐标点，坐标点对应的抛物线分别为${\left({x - r} \right)^2} + f\left(r \right)$，${\left({x - q} \right)^2} + f\left(q \right)$，${\left({x - z} \right)^2} + f\left(z \right)$则相邻抛物线的交点为

$ \begin{array}{l} {s_1} = \frac{{\left( {f\left( r \right) + {r^2}} \right) - \left( {f\left( q \right) + {q^2}} \right)}}{{2r - 2q}}\\ {s_2} = \frac{{\left( {f\left( q \right) + {q^2}} \right) - \left( {f\left( z \right) + {z^2}} \right)}}{{2q - 2z}} \end{array} $

(14)

若交点$s_2$在$s_1$的左边则抛物线$q$不构成抛物线的下界而被删除，这样采用一个堆栈结构从左到右扫描，将构成下界的抛物线放进堆栈，将不构成抛物线下界的抛物线弹出堆栈，即能快速计算式(13)，对于单个$v_q$的计算复杂度仅为${\rm{O}}\left({2{L}} \right)$，总体计算复杂度为${\rm{O}}\left({2{L^2}} \right)$，最终式(12)的计算复杂度只有${\rm{O}}\left({4{L^2}} \right)$。

3) 确定光流：信息传递结束后，对于每个像素点，找到使式(9)最小的那个$q$，即得到每个点的光流标记值(${u_q}, {v_q}$)。

假设金字塔层数为3，每层迭代次数为$T$，图像大小为$N = H \times W$，则1.5节中步骤1)的计算复杂度为${\rm{O}}\left({N \times {L^2}} \right)$，步骤2)采用二分法的计算复杂度为${\rm{O}}\left({0.5N \times 4{L^2}\left({1 + 1/4 + 1/16} \right) \times 4T} \right) = {\rm{O}}\left({10.5N \times {L^2} \times T} \right)$，其中$4{L^2}$是单个像素的信息计算复杂度，步骤3)的计算复杂度为：${\rm{O}}\left({N \times {L^2}} \right)$，则总的计算复杂度为${\rm{O}}\left({N \times {L^2} \times \left({10.5T + 2} \right)} \right)$。可见其计算主要负荷在于信息计算。

常规的并行计算方法用一个线程处理每一个像素点，可以将总体计算复杂度${\rm{O}}\left({N \times {L^2} \times \left({10.5T + 2} \right)} \right)$中的涉及像素点数目的项设为1，则其计算复杂度为${\rm{O}}\left({2{L^2} + 48{L^2} \times T} \right)$。这里存在两个问题：一是计算复杂度仍然较高，不能满足快速计算要求；二是单个线程计算负荷过重，系统难以分配足够资源。为解决这些问题，本文在CUDA平台上提出了优化的马尔可夫光流并行计算方法。

本文主要对最为复杂的信息计算及信息更新部分的并行计算做了创新性设计。由于需要对上下左右4个方向进行信息更新，而每个方向上的更新是独立的，为此本文同时在4个方向进行信息计算和更新，且对于最为复杂的信息计算步骤，采用了计算复杂度只有${\rm{O}}\left({4L} \right)$的并行降维法^[18]，这样3层金字塔总体计算复杂度为${\rm{O}}\left({2{L^2} + 12L \times T} \right)$，比常规并行方法少得多。

本并行计算方法计算过程如图 5所示，和串行方法一样共有3个步骤。CUDA中主要通过分配线程块和线程来实现并行计算，通常在处理图像时，将一幅图像分成由若干个线程块组成的2维网格，每个线程块中再分配若干个线程，每个线程对应一个图像像素点的操作。线程块和线程的划分对并行算法的效率有很大影响。下面分别介绍这3个步骤中并行方法是如何在CUDA中实现的。

1) 并行计算数据项，图像中每个点的数据项的计算都是独立的，只要在并行运算平台CUDA中为每个点分配一个线程，每个线程并行地根据式(7)即可计算完数据项。具体分配线程时采用最简单的2维线程块网格进行线程划分，假设$W$、$H$分别是图像的宽和高，每个线程块能够分配${n_{TX}} \times {n_{TY}}$个线程，则需要分配${W_b} \times {H_b}$个线程块，其中${W_b} = W/{n_{TX}}$，${H_b} = H/{n_{TY}}$，在CUDA编程时，进行线程块和线程分配，即

$ {\rm{Dim3}}\;{\rm{Block}}\left( {{W_b},{H_b}} \right) $

(15)

$ {\rm{Dim3}}\;{\rm{Threads}}\left( {{n_{TX}},{n_{TY}}} \right) $

(16)

进行线程与像素点映射，即

$ x = blockIdx.x \times blockDim.x + threadIdx.x $

(17)

$ y = blockIdx.y \times blockDim.y + threadIdx.y $

(18)

式中，$blockIdx.x$和$blockIdx.y$为某个线程块在网格中的$X$和$Y$方向的索引，$threadIdx.x$和$threadIdx.y$为某个线程块中的线程分别在$X$和$Y$方向的索引。$blockDim.x= W_b$，$blockDim.x= H_b$。这样每一个线程均对应一个像素点，同时计算每个像素点的数据项，整体计算复杂度为${\rm{O}}\left({{L^2}} \right)$。

2) 并行计算信息和二分法传递信息，本文主要对最为复杂的信息计算及信息传递部分的并行计算做了创新性设计。由于需要对上下左右4个方向进行信息更新，在分配线程块的时候采用了式(19)在CUDA中分配了如图 6所示的一个3维的线程块结构

$ {\rm{Dim3}}\;{\rm{Blocks}}\left( {{W_b}/2,{H_b},4} \right) $

(19)

该结构中共有4层，每一层负责计算一个方向的所有像素点的信息。由于在计算信息的时候每个线程的计算负荷很重，故本文分配多个线程采用并行的降维法计算一个像素点的信息，这样处理大大减少了每个线程的计算负荷，提高了计算效率。为此本文在线程块中分配线程时，在每一层中为每个线程块分配了一个3维线程结构

$ {\rm{Dim3}}\;{\rm{Threads}}\left( {{n_{TX}},{n_{TY}},L} \right) $

(20)

同样有${W_b} = W/{n_{TX}}, {H_b} = H/{n_{TY}}$，$L$对应$X$和$Y$方向上的采样点数，即图 6中最小方块所包含的线程数。每个大方块包含的小方块的数目由$n_{TX}$, $n_{TY}$来决定。为实现二分法信息传递，本文根据式(21)和式(22)利用线程块前2维索引($blockIdx.x$，$blockIdx.y$)和线程的前2维索引($threadIdx.x$，$threadIdx.y$)计算要操作的图像像素点的坐标$x$和$y$(即每个小方块的位置)，使得每个小方块对应着图像中的一个像素点。式(22)中的$t$为当前的循环次数，可以看到当$t+y$为奇数时对应图 6中的黑点，当$t+y$为偶数时对应图 6中的白点，这样每次循环时只计算一半的数据，实现了二分法的并行处理

$ y = blockIdx.y \times blockDim.y + threadIdx.y $

(21)

$ \begin{array}{*{20}{c}} {x = \left( {blockIdx.x \times blockDim.x + } \right.}\\ {\left. {threadIdx.x} \right) \times 2 + \bmod \left( {t + y,2} \right)} \end{array} $

(22)

$ {z_t} = threadIdx.z $

(23)

$ {z_b} = blockIdx.z $

(24)

图 6中每个小方块中又包含$L$个线程，这$L$个线程同时计算一个像素点的信息。每个小方块中的线程由第3维线程索引$z_t$决定，用这$L$个线程来并行计算2维降维法中的每一行或每一列的值。下面介绍并行降维法。

在降维法计算中，图像的$X$方向的每行和$Y$方向的每列计算是独立的，可以并行完成。并行降维法中包含两个步骤：

(1) 计算$m_{h\left| {{u_p}} \right.}^t\left({{v_q}} \right) = \mathop {\min }\limits_{{v_p}} \left({{{\left({{v_q} - {v_p}} \right)}^2} + h\left({{u_p}, {v_p}} \right)} \right)\left({p = 1, 2, \cdots, L; q = 1, 2, \cdots, L} \right)$。在图 6中的每个小方块中分配$L$个线程，每个线程对应一个纵向的$v_p$，即第$i$个线程对应的$v_p$=$i$，在横向每行依次从左向右进行扫描，即$u_p$从1变化到$L$，每个扫描点对应一个抛物线${\left({{v_q} - {v_p}} \right)^2} + h\left({{u_p}, {v_p}} \right)$，显然这些抛物线的下界构成$m_{h\left| {{u_p}} \right.}^t\left({{v_q}} \right)$。搜索抛物线下界的计算复杂度为${\rm{O}}\left(L \right)$，再根据抛物线的下界计算最小值，其计算复杂度也为${\rm{O}}\left(L \right)$，则该步骤的计算复杂度为${\rm{O}}\left({2L} \right)$。

(2) 计算$m_{pq}^t\left({{u_q}, {v_q}} \right)$，在图 6中的每个小方块中分配$L$个线程，每个线程对应一个横向的$u_p$，在纵向每列依次从上向下进行扫描，每个扫描点对应的抛物线为${\left({{u_q}, {v_p}} \right)^2} + m_{h\left| {{u_p}} \right.}^t\left({{v_q}} \right)$，同样找出这些抛物线的下界即可得到$m_{pq}^t\left({{u_q}, {v_q}} \right)$，则该步骤的计算复杂度也为${\rm{O}}\left({2L} \right)$。综合两步骤，总的计算复杂度为${\rm{O}}\left({4L} \right)$。

在信息传递时，根据式(24)，$z_b$对应信息更新的每一个方向，其中0表示向上，1表示向下，2表示向右，3表示向左。这样相当于为每个像素点分配了4$L$个线程，计算和更新信息时，每个方向的计算和降维法的每行(列)的计算都是并行的，大大提高了计算速度。

3) 并行输出各点光流值($u$，$v$)，图像中每个点的光流值$u$和$v$的输出计算都是独立的，只要在并行运算平台CUDA中为每个点分配一个线程，每个线程并行地根据式(9)即可计算得到输出光流值$u$和$v$，其计算复杂度为${\rm{O}}\left({{L^2}} \right)$。

在并行计算中，不同线程之间的同步至关重要，如果需要同步的地方没有同步，可能会出现随机结果。BP算法中采用的二分法信息传递不仅减少了计算量，还很好地解决了线程同步问题。由图 2和图 6可知一个像素点的信息计算和更新需要周围3个邻域点的信息，采用二分法能够保证信息的读取和更新是严格分离的。举个例子，在某次迭代中，更新白点的信息，则首先要读取邻域3个黑点的信息，而这3个黑点的信息由上次迭代计算得到，在本次迭代中不会被更新，随后在更新白点的信息时，白点的信息不会被任意一个邻域点读取，这样均保证了数据的一致性。于是本并行算法只要在信息更新或读取后插入同步指令即可实现数据同步，而不会出现随机读取及读写冲突现象。

本文的实验平台为联想台式计算机，Windows7 64位操作系统，intel i7-4790 CPU，主频3.6 GHz，8 GB内存，Nvidia GTX950显卡，算法均在CUDA3.2平台C语言实现。本文实验中选用的性能评估指标由Barron等人^[19]提出，分别是平均角速度误差(AAE)和平均端点误差(AEE)。AAE是计算光流矢量($\left({{u^\text{e}}, {v^\text{e}}} \right)$)与标准光流矢量($\left({{u^\text{c}}, {v^\text{c}}} \right)$)的平均角度，其表达式为

$ A = \frac{1}{N}\sum\limits_{i = 1}^N {{\phi _{\rm{e}}}\left( i \right)} $

(25)

$ \begin{array}{*{20}{c}} {{\phi _{\rm{e}}}\left( i \right) = }\\ {\arccos \left( {\frac{{u_i^{\rm{c}}u_i^{\rm{e}} + v_i^{\rm{c}}v_i^{\rm{e}} + 1}}{{\sqrt {{{\left( {u_i^{\rm{c}}} \right)}^2} + {{\left( {v_i^{\rm{c}}} \right)}^2} + 1} \sqrt {{{\left( {u_i^{\rm{e}}} \right)}^2} + {{\left( {v_i^{\rm{e}}} \right)}^2} + 1} }}} \right)} \end{array} $

(26)

AEE表示($\left({{u^e}, {v^e}} \right)$)和($\left({{u^c}, {v^c}} \right)$)之间的平均绝对误差，其表达式为

$ E = \frac{1}{N}\sum\limits_{i = 1}^N {\sqrt {{{\left( {u_i^{\rm{c}} - u_i^{\rm{e}}} \right)}^2} + {{\left( {v_i^{\rm{c}} - v_i^{\rm{e}}} \right)}^2}} } $

(27)

共采用3种图像序列与Horn算法^[1]、Weickert算法^[3]、Hossen^[15]提出的Lucas并行光流算法和Grauer-Gray等人^[16]提出的马尔可夫并行光流算法进行对比实验。

本实验采用合成旋转小球图像序列的第19帧和20帧进行实验，以验证变分光流算法对单一的旋转运动物体的处理能力。图像大小为200×200像素，其图像、真实光流及计算光流如图 7所示，计算误差见表 1。

图 7(c)为Horn算法计算得到的光流场结果，在球的内部光流矢量在出现了较明显的偏移, 其原因在于这些区域并不满足灰度守恒假设；图 7(d)为Weickert算法计算得到的光流场结果，在球的外部出现了数值非零的光流矢量，导致计算误差较大，其原因在于此区域较为平滑，梯度较小，梯度守恒假设难以满足；图 7(e)为Hossen算法计算得到的光流场结果，其计算时间最少，有些点出现较大误差，其原因在于Lucas算法在利用矩阵求逆计算这些区域的光流时，特征矩阵的行列式过小，导致计算结果错误，说明Lucas算法稳定性不佳。图 7(f)为Grauer-Gray算法计算得到的光流场结果，其计算时间小于1 s，结果中有较大的误差，其原因在于数据插值引起的误差；图 7(g)为本文算法计算得到的光流场结果，虽然由于MRF的平滑作用，在球的顶部有少许小光流矢量丢失，但是仍得到了最低的计算误差，由于对并行计算进行了优化，计算时间是Grauer-Gray算法的一半。表 1中还列出了串行和并行计算时间，可以看出采用并行算法，计算效率提高了314倍。

本实验采用Yosemite山谷真实图像序列的第10帧和11帧进行实验，以验证本文算法对有光照变化的运动图像的光流场计算精度。图像大小为316×252像素，其图像、真实光流及计算光流如图 8所示，计算误差见表 2。

图 8(c)为Horn算法计算得到的光流场结果，Horn光流场算法误差较大，虽然对物体的大致运动有着正确的判断，但在光照变化区域出现了比较大的误差，这说明采用灰度守恒假设的Horn算法对光照变化敏感；图 8(d)为Weickert算法计算得到的光流场结果，由于Weickert算法中采用了梯度守恒假设，对光照变化不敏感，故能正确反映图像中的物体运动情况，得到的误差较小。图 8(e)为Hossen算法计算得到的光流场结果，计算时间最少，同样由于Lucas算法的稳定性问题，在光照变化处有较大误差；图 8(f)为Grauer-Gray算法计算得到的光流场结果，由于仅采用了灰度信息，故在光照变化处有较大角度误差。图 8(g)为本文算法计算得到的光流场结果，由于既采用了灰度信息又采用了梯度信息，虽然在光照变化区域出现了一定的误差，但是计算误差比Horn算法要小很多，得到最小的AEE, AAE也仅略大于Weickert和Hossen算法。

本实验采用来自Hydrangea图像序列(http://vision.middlebury.edu/flow/data/)的第10帧和第11帧进行实验，以验证变分光流算法处理多物体复杂运动图像序列的计算能力。图像大小为584×388像素，其图像、真实光流及计算光流如图 9所示，计算误差见表 3。

一般情况下，不同物体之间的遮挡处可能既不满足灰度守恒也不满足梯度守恒，故在处理多运动图像时，图 9(c)中Horn算法计算得到的光流场结果和图 9(d)中Weickert算法计算得到的光流场结果误差较大。图 9(e)为Hossen算法计算得到的光流场结果，其计算时间最少，由于稳定性问题，其中有些点出现较大误差。图 9(f)为Grauer-Gray算法计算得到的光流场结果，在遮挡处存在较大误差。图 9(g)为本文算法计算得到的光流场结果，从光流图上可以看出，在遮挡处由于采用了非平方惩罚函数通过平滑运动矢量，本方法具有更好的效果，得到最小的AEE，AAE也仅略大于Hossen算法。

本文提出的光流计算模型中最重要的参数是式(1)中的$\lambda $与式(7)中的$\gamma $和$\alpha $。其中$\lambda $表示平滑项和数据项之间的比例关系，$\lambda $越小，平滑权重越高，得到的光流越光滑。$\gamma $和$\alpha $和光流位移量有关。当位移较小时，插值计算容易引起较大误差，应选用较大的$\gamma $，反之应选用较小的$\gamma $。当位移较小时，应选用较小的$\alpha $，以提高计算精度，当位移较大时，应选用较大的$\alpha $，以减少计算量。

图 10中列出了6组不同参数下对Yosemite山谷真实图像序列的第10帧和11帧进行实验的结果。可以看到当$\lambda $较小时(图 10(a))，得到的光流较平滑，在云彩光照变化大区域得到的误差较小；当$\lambda $较大时(图 10(b))，平滑作用小，在云彩光照变化大区域得到的误差较大。当$\alpha $较小时(图 10(c))，在光流位移大的区域(左下角)误差较大；当$\alpha $较大时(图 10(d))，在光流位移小的区域(中间偏右)得到的误差较大；当$\gamma $较小时(图 10(e))，在光流位移较小的区域(中间偏右)误差较大；当$\gamma $较大时(图 10(f))，在光流位移较大的区域(左下角)误差较大。本实验结果与理论分析基本吻合。

本文提出了一种基于马尔可夫随机场的并行变分光流计算方法，从实验结果来看，其计算时间接近Horn算法，远远小于Weickert算法，快于常规并行马尔可夫光流计算方法(Grauer-Gray算法)，仅比并行Lucas光流计算方法(Hossen算法)慢。可见本文实现了MRF并行变分光流的高效计算。从3个实验结果来看，Horn算法和Grauer-Gray算法仅采用了灰度信息进行光流计算，对光照变化敏感，对多运动目标之间的遮挡处理能力弱；Weickert算法采用了梯度信息进行光流计算，对光照变化不敏感，但是在运动边界常出现大的误差；Lucas算法在3个实验中的某些区域均会出现较大误差，存在稳定性问题。本文算法无论是处理旋转运动、多运动目标和光照变化情况，均取得最小的AEE误差。可见提出的马尔可夫并行变分光流计算模型具备一定的鲁棒性。虽然本文方法取得了很高的计算效率和较好的精度，但是本文方法对内存需求巨大，在处理高分辨率图像时，限制了其采样点数，难以计算大位移。故需对MRF光流计算方法做进一步研究，解决高分辨率图像大位移计算问题。