一种基于误差分布式多层网格的快速光流场计算方法

摘要

本发明属于机器视觉领域，公开了一种基于误差分布式多层网格的快速光流场计算方法，包括：步骤一，输入图像，构建线性方程Ax=f；步骤二，建立多层图像金字塔；步骤三，进行前优化，消除高频分量；步骤四，进行残差传递，消除低频分量；步骤五，重复步骤三、四，直至残差传递至最粗层；步骤七，从最粗层开始进行误差回传；步骤八，进行细网格误差校正；步骤九，进行后优化迭代，提高解的稳定性；步骤十，重复步骤七～九，直至残差传递至最细层。本发明所述的方法是一种加速方程优化求解的有效方法，高频误差收敛很快，能够显著提高视觉光流场的计算速度，与变分法相比，该方法的收敛速度能够提高到3.5倍以上。

说明书

技术领域

本发明属于机器视觉领域，涉及一种基于误差分布式多层网格的快速光流场计算方法。

背景技术

在视觉运动分析理论中，当相机与场景目标间有相对运动时所观察到的图像亮度模式运 动称之为光流(optical flow)。由于序列视频图像之间的光流矢量场表达了图像的变化，能够定 量地描述帧间同名像点之间的匹配关系，并给出像面偏移量，因此为观察者提供了关于目标 的运动及结构信息。

变分法是求解光流的主流方法。其指导思想是在梯度平滑和数据保真的意义上构建能量 泛函，重点通过求解最优化能量函数的欧拉-拉格朗日方程(偏微分方程)得到致密的光流矢量 场。变分光流的优势在于，模型的通用性好，求解的精确度高，且光流场致密。在变分框架 下，Horn和Schunck提出的方法可以根据不同的假设条件得到不同的光流计算模型，其能量 泛函模型包括平滑项和数据项，可解析描述为：

$\underset{u}{\min }\{\underset{\Omega }{\int }\alpha ({|▿u|}^2+{|▿v|}^2)\mathrm{d\Omega }+\underset{\Omega }{\int }{(I_1(x+u\left(x\right))-I_0\left(x\right))}^2\mathrm{d\Omega }\}---\left(1\right)$

式中，I₀、I₁分别表示摄像机采集的相邻两帧图像，x=(x,y)^T表示图像上的某个像点， 用u(x)=(u(x),v(x))^T表示该点位置的光流矢量，其中u(x)和v(x)分别为该矢量的水平和垂直分 量。符号▽为梯度算子，Ω表示成像平面。前面的积分项为平滑项，后面的积分项为数据项， α为调节二者权重的常系数。

将上式灰度守恒项I₁(x+u(x))-I₀(x)进行泰勒展开，得到：

I₁(x+u(x))-I₀(x)＝I_xu+I_yv+I_t           (2)

式中，x、y分别表示图像I的水平与垂直分量，t表示时间。右下角下标表示相应的偏 导数，即$I_x=\frac{\partial I}{\partial x},$$I_y=\frac{\partial I}{\partial y},$$I_t=\frac{\partial I}{\partial t}·$

将式(2)带入式(1)，上述能量泛函转换为：

$\underset{u}{\min }\{\underset{\Omega }{\int }({|▿u|}^2+{|▿v|}^2)\mathrm{d\Omega }+\lambda \underset{\Omega }{\int }{(I_{x}u+I_{y}v+I_t)}^2\mathrm{d\Omega }\}$

利用变分方法，得到其相对应的欧拉-拉格朗日方程为：

$\left(\begin{array}{c}{I}_x^2·u+I_x{I}_y·v-\mathrm{\alpha \Delta u}=-I_x{I}_t\\ I_x{I}_y·u+I_y^2·v-\mathrm{\alpha \Delta v}=-I_y{I}_t\end{array}\right)---\left(3\right)$

式中，u和v分别为光流矢量的水平和垂直分量，符号Δ表示拉普拉斯算子。该方程组 离散化后，若用传统的数值求解法，如Gauss–Seidel，SOR等方法需要迭代数千次才能求得 较理想的结果，因此算法的实时性较差。

在传统数值迭代算法中，分析误差的收敛过程可以发现，高频误差分量在迭代过程中迅 速衰减，而处在低频的误差分量却衰减缓慢。因此，可用一系列分辨率不同的多层网格来快 速求解同一数值问题。这是因为，在细网格上利用较少的迭代次数即可消除高频误差；遗留 的低频误差相对于粗网格转变为高频误差，将映射至较粗的网格上，可继续在粗网格上利用 较少的迭代次数消除这部分低频误差，已达到加速收敛的目的。通常设定多层网格，将误差 层层传递，实现误差的分布式估计。在具体实现时，需将式(3)转换为线性方程组。将式(3) 稍加变形，可得到如下的形式：

β(x)=f                          (4)

其中：

$\beta =\left(\begin{array}{cc}{I_x}^2& I_x{I}_y\\ I_x{I}_y& {I_y}^2\end{array}\right)+\left(\begin{array}{cc}-\mathrm{\alpha \Delta }& 0\\ 0& -\mathrm{\alpha \Delta }\end{array}\right)$

$f=\left(\begin{array}{c}-I_x{I}_t\\ -I_y{I}_t\end{array}\right)$

$x=\left(\begin{array}{c}u\\ v\end{array}\right)$

可见，函数β中不仅包含常系数，还包含非线性的拉普拉斯算子由此导致 光流场的计算不能直接应用多层网格算法求解。

发明内容

根据光流计算方法中存在的上述问题，本发明提出了一种基于误差分布式多层网格的快 速光流场计算方法，可以加速光流迭代计算的收敛速度，提高计算机视觉系统的实时性。

多层网格是一种加速方程优化求解的有效方法。其中，构造方程Ax=f是首要的一步。 由于在光流方程(4)中，函数β中含有非线性的拉普拉斯算子，因此本发明将函数β转化为不 含算子的系数矩阵A，以便在多层网格算法框架下进行求解。

一种基于误差分布式多层网格的快速光流场计算方法，其特征在于包括如下步骤：

步骤一，输入图像，构建线性方程Ax=f，方法如下：

(1)计算时空梯度张量矩阵，将待解方程组离散化；

假设在某一帧图像平面上,（i,j）表示第i行第j列像素点的位置，u_i,j和v_i,j分别表示该 点在相邻帧间发生的横向x轴和纵向y轴上的位移场。计算图像的空间偏导数f_x,f_y以及时间 偏导数f_t,可得到图像中每个像素的时空梯度张量J矩阵：

$J=\left(\begin{array}{ccc}{f}_x{f}_x& f_x{f}_y& f_x{f}_t\\ f_y{f}_x& f_y{f}_y& f_y{f}_t\\ f_t{f}_x& f_t{f}_y& f_t{f}_t\end{array}\right)=\left(\begin{array}{ccc}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)---\left(5\right)$

若用[J]_i,j表示J矩阵在像素点(i,j)处的值，用h表示空间步长，则欧拉-拉格朗日方程(式 (3))可离散为：

${\left[J_{11}\right]}_{i,j}{u}_{i,j}+{\left[J_{12}\right]}_{i,j}{v}_{i,j}-\frac{\alpha }{h^2}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-{4u}_{i,j})=-{\left[J_{13}\right]}_{i,j}---\left(6\right)$

${\left[J_{12}\right]}_{i,j}{u}_{i,j}+{\left[J_{22}\right]}_{i,j}{v}_{i,j}-\frac{\alpha }{h^2}(v_{i+1,j}+v_{i-1,j}+v_{i,j+1}+v_{i,j-1}-{4v}_{i,j})=-{\left[J_{23}\right]}_{i,j}---\left(7\right)$

(2)将差分方程(6)转化成线性常系数方程组A₁x=f₁；

化简方程(6)，得到求解区域D内的差分方程组：

$(1+\frac{h^2}{4\alpha }{\left[J_{11}\right]}_{i,j})u_{i,j}-\frac{1}{4}(u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1})+\frac{h^2}{4\alpha }\left[J_{12}\right]v_{i,j}=-\frac{h^2}{4\alpha }{\left[J_{13}\right]}_{i,j}$

假设图像总共有N+1行，除去第零行和第N行，第一列各个像素的差分方程组有如下形 式：

$\left(\begin{array}{ccccccc}1+\frac{h^2}{4\alpha }{\left[J_{11}\right]}_{11}& -\frac{1}{4}& 0& ···& 0& 0& 0\\ -\frac{1}{4}& 1+\frac{h^2}{4\alpha }{\left[J_{11}\right]}_{21}& -\frac{1}{4}& ···& 0& 0& 0\\ ···& ···& ···& ···& ···& ···& ···\\ 0& 0& 0& ···& -\frac{1}{4}& 1+\frac{h^2}{4\alpha }{\left[J_{11}\right]}_{(N-2)1}& -\frac{1}{4}\\ 0& 0& 0& ···& 0& -\frac{1}{4}& 1+\frac{h^2}{4\alpha }{\left[J_{11}\right]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{u}_{11}\\ u_{21}\\ ···\\ u_{(N-1)1}\end{array}\right)$

$-\frac{1}{4}\left(\begin{array}{cccc}1& 0& ···& 0\\ 0& 1& ···& 0\\ ···& ···& ···& ···\\ 0& 0& ···& 1\end{array}\right)\left(\begin{array}{c}{u}_{12}\\ u_{22}\\ ···\\ u_{(N-1)2}\end{array}\right)+\left(\begin{array}{cccc}\frac{h^2}{4\alpha }{\left[J_{21}\right]}_{11}& 0& ···& 0\\ 0& \frac{h^2}{4\alpha }{\left[J_{21}\right]}_{21}& ···& 0\\ ···& ···& ···& ···\\ 0& 0& ···& \frac{h^2}{4\alpha }{{[J}_{21}]}_{(N-1)1}\end{array}\right)\left(\begin{array}{c}{v}_{11}\\ v_{21}\\ ···\\ v_{(N-1)1}\end{array}\right)$

$=\left(\begin{array}{c}-\frac{h^2}{4\alpha }{{[J}_{13}]}_{11}+\frac{1}{4}{g}_{01}+\frac{1}{4}{g}_{10}\\ -\frac{h^2}{4\alpha }{{[J}_{13}]}_{21}+\frac{1}{4}{g}_{20}\\ ···\\ \frac{h^2}{4\alpha }{{[J}_{13}]}_{(N-1)1}+\frac{1}{4}{g}_{N1}+\frac{1}{4}{g}_{(N-1)0}\end{array}\right)=b_1---\left(8\right)$

式中，g_ij表示求解区域D边界上的元素u_i，即g_ij＝u_ij。

若引入光流分量u_ij和v_ij上的列向量：

$u_i=\left(\begin{array}{c}{u}_{\mathrm{li}}\\ u_{2i}\\ ···\\ u_{(N-1)i}\end{array}\right),$$v_i=\left(\begin{array}{c}{v}_{\mathrm{li}}\\ v_{2i}\\ ···\\ v_{(N-1)i}\end{array}\right)$

且令：

${\left[J_{\mathrm{mn}}\right]}_i=\left(\begin{array}{cccc}{{[J}_{\mathrm{mn}}]}_{\mathrm{li}}& 0& ···& 0\\ 0& {{[J}_{\mathrm{mn}}]}_{2i}& ···& 0\\ ···& ···& ···& ···\\ 0& 0& ···& {{[J}_{\mathrm{mn}}]}_{(N-1)i}\end{array}\right)---\left(9\right)$

对于像素点(m,n)处的J矩阵元素的J_mn，式(9)中矩阵[J_mn]_i相当于N-1阶的单位阵与J_mn第i列这一列向量的乘积。因此，第一列各个像素的差分方程组(8)可化为：

$(G+\frac{h^2}{4\alpha }{{[J}_{11}]}_1)u_1-\frac{1}{4}{\mathrm{Iu}}_2+\frac{h^2}{4\alpha }{{[J}_{21}]}_1{v}_1=b_1$

式中，I为N-1阶的单位阵，b₁为式(6)中定义的列向量，G为下面的N-1阶矩阵：

$G=\left(\begin{array}{ccccccc}1& -\frac{1}{4}& 0& ···& 0& 0& 0\\ -\frac{1}{4}& 1& -\frac{1}{4}& ···& 0& 0& 0\\ ···& ···& ···& ···& ···& ···& ···\\ 0& 0& 0& ···& -\frac{1}{4}& 1& -\frac{1}{4}\\ 0& 0& 0& ···& 0& -\frac{1}{4}& 1\end{array}\right)$

同理，第二列各个像素的差分方程为：

$-\frac{1}{4}{\mathrm{Iu}}_1+(G+\frac{h^2}{4\alpha }{{[J}_{11}]}_2)u_2-\frac{1}{4}{\mathrm{Iu}}_3+\frac{h^2}{4\alpha }{{[J}_{21}]}_2{v}_2=b_2$

式中：

$b_2=\left(\begin{array}{c}-\frac{h^2}{4\alpha }{{[J}_{13}]}_{12}+\frac{1}{4}{g}_{02}\\ -\frac{h^2}{4\alpha }{{[J}_{13}]}_{22}\\ ······\\ -\frac{h^2}{4\alpha }{{[J}_{13}]}_{(N-1)2}+\frac{1}{4}{g}_{N2}\end{array}\right)$

将图像所有列均进行转化，则方程组可以写为A₁x=f₁的形式，式中：

$A_1=\left(\begin{array}{cccccccccccc}G+\frac{h^2}{4\alpha }{{[J}_{11}]}_1& \frac{1}{4}I& 0& ···& 0& 0& 0& \frac{h^2}{4\alpha }{{[J}_{21}]}_1& 0& ···& 0& 0\\ \frac{1}{4}I& G+\frac{h^2}{4\alpha }{{[J}_{11}]}_2& \frac{1}{4}I& ···& 0& 0& 0& 0& \frac{h^2}{4\alpha }{{[J}_{21}]}_2& ···& 0& 0\\ ···& ···& ···& ···& ···& ···& ···& ···& ···& ···& ···& ···\\ 0& 0& 0& ···& \frac{1}{4}I& G+\frac{h^2}{4\alpha }{{[J}_{11}]}_{N-2}& \frac{1}{4}I& 0& 0& ···& \frac{h^2}{4\alpha }{{[J}_{21}]}_{N-2}& 0\\ 0& 0& 0& ···& 0& \frac{1}{4}I& G+\frac{h^2}{4\alpha }{{[J}_{11}]}_{N-1}& 0& 0& ···& 0& \frac{h^2}{4\alpha }{{[J}_{21}]}_{N-1}\end{array}\right)$

$x=\left(\begin{array}{c}{u}_1\\ u_2\\ ···\\ u_{N-1}\\ v_1\\ v_2\\ ···\\ v_{N-1}\end{array}\right)$$f_1=\left(\begin{array}{c}{b}_1\\ b_2\\ ···\\ b_{N-1}\end{array}\right)$

(3)将差分方程(7)转化成线性常系数方程组A₂x=f₂；

(4)由A₁x=f₁与A₂x=f₂共同构成方程组：

Ax=f                    (10)

步骤二，建立多层图像金字塔；

多层图像金字塔的示意图如图2所示。最上面的一层图像为原分辨率图像，下面的各层 图像表示逐层缩小的降分辩率图像。由于像素本身即为均匀离散化的，因此可以将多层图像 看做多层网格。

通常，图像的尺寸（长和宽）介于300-800之间。按照设定的比例因子η（0.5<η<0.95） 将原分辨率图像逐层缩小，并将尺寸取整，得到多层图像金字塔。一般地，层数N=4～5层。 每层光流的初值设置为零。

步骤三，进行前优化，消除高频分量；

在当前层i(i=1，2，…，N-1)，以零矩阵为光流初值对细网格方程迭代m次，以消除高 频分量。设x_i为初值，得到近似估计值记为：

${\bar{x}}_i={\mathrm{Relax}}^m(x_i,A_i,r_i)$

式中，Relax^m表示Gauss-Seidel迭代过程，m为设定的迭代次数。

A_i由式(10)确定；当i=1时，r由式(10)确定，即r₁=f。否则，r等于上一层传递至本 层的残差。本层相应的残差更新为：

$r_i=A_i{\bar{x}}_i-f_i$

步骤四，进行残差传递，消除低频分量；

前优化过程旨在消除高频误差分量，而遗留误差主要为低频分量。因此，将残差传递至 粗网格上，以消除低频分量。限制到粗网格上的余量r_i+1为：

$r_{i+1}=I_i^{i+1}{r}_i$

式中，为细网格至粗网格上的映射算子。

步骤五，当i=1～N-1时，重复步骤三、四，直至残差传递至最粗层；

步骤六，求解粗网格上的方程组：

A_Ne_N＝r_N

式中，e_N为最粗层上精确求解的光流误差。

通过矩阵运算得到：

e_N＝(A_N)^-1r_N

步骤七，从最粗层(j=N)开始进行误差回传；

当前层j，将求得的误差校正量回传到上一层较细网格j-1上，即：

$e_{j-1}=I_j^{j-1}{e}_j$

式中，为粗网格至细网格上的映射算子。

步骤八，进行细网格误差校正；

将计算的初始值加上由粗网格传回的误差，得到细网格校正后的解即：

${\hat{x}}_{j-1}={\hat{x}}_{j-1}+e_{j-1}$

步骤九，进行后优化迭代，提高解的稳定性；

细网格校正之后，以校正后的解为初值，执行后优化迭代n次，得到更新后的校正值其表达式为：

$x_{j-1}^{\mathrm{new}}={\mathrm{Relax}}^n({\hat{x}}_{j-1},A_{j-1},f_{j-1})$

步骤十，当j=N，N-1，…，2时，重复步骤七～九，直至残差传递至最细层。

与现有技术相比，本发明具有以下有益效果：

本发明将光流求解方程按列展开，将其化为线性方程组的形式，使得光流的复杂计算过 程能够通过线性多层网格算法加速实现。其次，采用了不同分辨率的多层网格，将不同频率 的误差分量分布在不同的网格上，使得误差相对于本层网格成为高频误差。由于高频误差收 敛的很快，因而该方法能够显著提高视觉光流场的计算速度，与变分法相比，该方法的收敛 速度能够提高到3.5倍以上。

附图说明

图1为本发明所涉及方法的流程图；

图2为多层图像金字塔示意图；

图3为本发明应用实例的第一组测试图像RubberWhale；

图4为本发明应用实例的第一组测试图像。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步说明。

本发明的实施依托P3800M主频、512M内存及以上配置的微机一台，做为硬件设备平 台。运行环境为Windows XP操作系统和Matlab软件平台。

本实施方式的具体过程为，首先根据待求解的光流模型，构造线性方程组Ax=f。然后建 立4～5层图像金字塔，在每一层分别进行优化迭代，并将迭代残差逐层映射到下一层更粗的 网格上。在最粗层，求解最粗网格上的方程组，得到精确的误差估计。此后，沿相反方向回 传误差，在较细的网格上逐层进行误差校正。图1是本发明所涉及方法的流程图，具体包括 以下步骤：

步骤一，输入图像，构建线性方程Ax=f；

步骤二，建立多层图像金字塔；

步骤三，进行前优化，消除高频分量；

步骤四，进行残差传递，消除低频分量；

步骤五，当i=1～N-1，重复步骤三、四，直至残差传递至最粗层；

步骤六，求解粗网格上的方程组：A_Ne_N＝r_N

步骤七，从最粗层(j=N)开始进行误差回传；

步骤八，进行细网格误差校正；

步骤九，进行后优化迭代，提高解的稳定性；

步骤十，当j=N，N-1，…，2时，重复步骤七～九，直至残差传递至最细层。

下面给出本发明的一个应用实例。

选择两组测试图像对本发明提出的基于误差分布式多层网格的快速光流场算法进行验 证。一组是国际上广泛采用的Middlebury标准库中的一对测试图像RubberWhale，如图3所 示；另一组是Wedel等人在他们发表的学术论文“An improved algorithm for TV-L¹optical flow” 中采用的一对测试图像，如图4所示。

实验时，分别采用传统变分算法和本发明提出的基于误差分布式多层网格的快速光流场 算法计算每组图像对之间的光流场。通过软件统计两种方法的运行时间，实验结果如表1所 示。表1中，实验一采用图3所示图像，实验二采用图4所示图像。由表1可知，基于误差 分布式多层网格的快速光流场算法相比传统变分算法的运算时间大大缩短，两个实验的收敛 速度均提高到3.5倍以上，显著提高了计算效率。

表1传统变分算法与本发明所述方法运行时间的对比