CN104869387A - 基于光流法的双目图像最大视差获取方法 - Google Patents

Abstract

基于光流法的双目图像最大视差获取方法，通过对左右摄像机采集到的图像进行一系列处理，得到图像上各点的运动矢量，对运动矢量进行分析判断最后得到最大视差，以便于进行下一步处理。由于光流法对弱纹理及重复纹理，噪声干扰等的影响较不敏感，对光流算法中数据项的梯度算子以及窗口权重函数的改进。考虑到最大视差唯一性的特点，分简单和复杂两次计算光流，有效减小计算量，加快了计算速度。在CLG-TV金字塔光流法的基础上，对光流算法中平滑项引入全局扩散系数，对光流算法中数据项的梯度算子、窗口权重函数进行改进，使其更加适用于计算双目图像只有水平位移的情况，对金字塔计算方法，分简单和复杂两次计算光流，有效减小计算量，加快了计算速度。

Description

基于光流法的双目图像最大视差获取方法

技术领域

本发明涉及一种基于光流法的双目图像最大视差获取方法，属于计算机视觉领域。

背景技术

视觉是人类了解世界、认知世界的重要方式，人类视觉系统通过双眼成像来感知现实世界，具有从二维景物图像感知三维信息的能力。立体视觉系统类似地利用双目相机模拟人的双眼来获取外部信息，利用左右图像获取视差图，视差图可以用于三维重建，物体识别与跟踪，3D视频编码等领域，非常具有研究价值。在某些情况下，我们更关心的是视差中的最大值即最大视差。

现有的双目图像最大视差确定方法一般是由立体匹配算法得到，立体匹配算法通过寻找同一实际场景通过不同视角获得的图像投影点之间的一一对应关系，有局部立体匹配和全局立体匹配等方法。立体匹配方法最终得到该实际场景的视差图，然后在视差图中获取最大值。

普通匹配方法存在问题，得到的视差图中存在很多斑块，这是由于左右摄像机采集到的图像中存在大量弱纹理结构区域引起的。实际上由于光照变化、深度不连续、存在弱纹理及重复纹理和噪声干扰等的影响，同一景物在左右视点下的成像会有很大差异，会对获得精确的视差值造成困难，鲁棒性有待提高，同时立体匹配的计算量较大，使得最大视差值获取的实时性有待提高。

发明内容

本发明的主要目的是提供一种基于光流法的双目图像最大视差获取方法。

不同于传统的立体匹配方法，本发明的具体创新点在于采用了光流法来获取最大视差。通过对左右摄像机采集到的图像进行一系列处理，得到图像上各点的运动矢量，对运动矢量进行分析判断最后得到最大视差，以便于进行下一步处理。由于光流法对弱纹理及重复纹理，噪声干扰等的影响较不敏感，可以得到较准确的最大视差值。

本发明的技术方案具体来说，主要包括如下技术内容：

1、由粗到精的多尺度计算：按照一定的采样标准，把原始图像分为多层，形成金字塔形式的图像结构，原始图像位于金字塔底部，而金字塔顶部则是图像分层后最低分辨率的近似。这样的多分辨率思想在使用光流分析大位移运动时特别适用，同时使用较小的取样因子使计算量减少，加快了计算速度。

2、CLG-TV光流法求解光流场。使用CLG-TV光流法处理左右图像获取光流运动矢量。

3、对光流算法中平滑项引入全局扩散系数，有效地在不连续处进行精确计算。

4、对光流算法中数据项的梯度算子以及窗口权重函数的改进。使用一种更加适用于计算最大视差的梯度算子的计算方法，并通过改进的双边滤波器来计算窗口权重函数。

5、考虑到最大视差唯一性的特点，分简单和复杂两次计算光流，有效减小计算量，加快了计算速度。

本发明的流程图如图1所示，实施流程如下：

步骤1、使用左右摄像机采集到的双目图像序列输入计算机，首先对输入图像做预处理，如果源图像的分辨率较大可以先做下采样使图像缩小到VGA级别以减少计算量，由于本方法不使用图像的彩色信息，首先将彩色图转为灰度图，公式如下：

I＝0.299R+0.587G+0.114B

然后对图像进行高斯平滑操作来消除摄像头采集过程中产生的随机噪声。

步骤2、使用CLG-TV光流法求解光流场，得到初步估计。

光流用来描述物体在观测到的图像中的二维运动，光流法运动估计的计算结果为图像中像素的运动矢量。理想情况下，对双目图像计算光流场将得到图像中像素点的水平运动矢量，这些水平运动矢量就等于立体匹配算法得到的各像素点的视差，可见光流法计算的运动矢量包含有视差信息，可以应用于双目图像最大视差的获取。光流法的基本假设将亮度约束和平滑约束相结合，既保证了运动矢量计算的准确性，又使得各区域内运动矢量能够平滑过渡，是一种高精度的运动估计方法。实际图像存在复杂多变的情况，同一景物在左右视点下的成像会有很大差异，会对获得精确的视差值造成困难。利用光流法求解最大视差的优势就在于光流法的精度能够达到亚像素级，完全可以胜任于最大视差的计算；光流法对传统方法中比较难以处理的深度不连续、弱纹理及重复纹理，噪声干扰等问题的鲁棒性强。

CLG-TV(combined local-global with total variation regularization)光流法是结合全局与局部的光流法，基于全局的光流法对小运动精确度高，基于局部的光流法对大运动有更好的精度，将二者结合起来可以增强鲁棒性，CLG-TV光流法使用的误差计算式将LK方法和HS方法的误差式结合起来。另外，HS和LK的误差计算式都采用了L2型范数。L2型范数有一些缺点：对噪声敏感度高，过度扩张，选择性差，以及不保留运动边界等等，为克服L2模型的缺点提出了采用L1型范数的方法将误差式改写，其形式为：

$E_{\mathrm{CLG}-\mathrm{TV}}=\underset{\mathrm{\Ω}}{\∫}[\left(\underset{\mathrm{region}}{\mathrm{\Σ}}W\left(x\right){(I_{x}u+I_{y}v+I_t)}^2\right)+\mathrm{\λ}\left(\right||\▿u||+||\▿v|\left|\right)]---\left(1\right)$

设I(x,y,t)是图像点(x,y)在时刻t的亮度，u(x,y)和v(x,y)是该点光流的x和y分量，假定点在t+δt时运动到(x+δx,y+δy)时，亮度保持不变，即

I(x+uδt,y+vδt,t+δt)＝I(x,y,t)

假设亮度随x、y、t光滑变化，将上式的左边用泰勒级数一阶展开：

$I(x,y,t)+\mathrm{\δx}\frac{\∂I}{\∂x}+\mathrm{\δy}\frac{\∂I}{\∂y}+\mathrm{\δt}\frac{\∂I}{\∂t}+e=I(x,y,t)$

其中e是关于δx、δy、δt的二阶和二阶以上的项。上式两边的I(x,y,t)相互抵消，两边除以δt，并取极限δt→0，得到这就是亮度不变性约束。

(1)式中左式是LK方法的数据项，其本质是一个区域上的所有像素点满足亮度不变性约束I_xu+I_yv+I_t＝0。其中亮度对空间的梯度I_x和I_y的物理意义是亮度在水平、竖直方向的变化率，亮度对时间的梯度I_t的物理意义是亮度随时间的变化率，u和v代表像素点分别在水平，竖直方向的变化速度。

(1)式中的右式是L1型的平滑项，由于实际场景中的物体都是连续，其运动矢量理应是大面积相同或相近的，因此设计平滑项这一项来满足这一特点。表示水平运动速度的空间梯度，反映了水平运动速度u在空间上的连续程度，L1型范数是指模的大小，其值越小连续程度越高；竖直方向同理。整个平滑项的物理意义是指求得的整个运动场应满足所有像素点的之和最小，即整个运动场的大部分点是足够连续，足够平滑的。这样就克服了深度不连续、弱纹理及重复纹理和噪声干扰造成的影响。

对于亮度对空间的梯度I_x和I_y，采用了Barron算子与horn算子结合的方式。Barron算子是一维的，其形式为1/12(180-8-1)；horn算子是三维的，它的计算方法如下式所示：

$I_x=\frac{1}{4}(I_{i+1,j,k}+I_{i+1,j+1,k}+I_{i+1,j,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i,j+1,k}+I_{i,j,k+1}+I_{i,j+1,k+1})$

$I_y=\frac{1}{4}(I_{i,j+1,k}+I_{i+1,j+1,k}+I_{i,j+1,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i+1,j,k}+I_{i,j,k+1}+I_{i+1,j,k+1})$

$I_t=\frac{1}{4}(I_{i,j,k+1}+I_{i+1,j,k+1}+I_{i,j+1,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i+1,j,k}+I_{i,j+1,k}+I_{i+1,j+1,k})$

Barron算子适合于计算空间梯度，horn算子适合于计算时间梯度。对于亮度对空间的梯度I_x和I_y，将Barron算子与horn算子计算出的梯度进行加权作为最后的梯度，对于亮度对时间的梯度I_t，将horn算子算出的梯度作为最后的梯度。

通过改进的双边滤波器来计算窗口权重函数，增强数据项中对边缘的保留，在双边滤波器的基础上增加阻塞度一项，使得权重正比于空间距离，亮度差别，阻塞度差别：

$\mathrm{BFW}\left(x\right)=\exp \{-{\mathrm{\λ}}_s\frac{(x-x_0)}{{2\mathrm{\σ}}_1^2}-{\mathrm{\λ}}_i\frac{(I_x-I_{x_0})}{{2\mathrm{\σ}}_2^2}\}{\mathrm{\λ}}_o\frac{o\left(x\right)}{o\left(x_0\right)}---\left(2\right)$

其中，阻塞度o(x)定义为

$o\left(x\right)=\exp \{-\frac{d\left(x\right)}{{2\mathrm{\σ}}_d^2}-\frac{I_t\left(x\right)}{{2\mathrm{\σ}}_e^2}\}---\left(3\right)$

其中d(x)为速度向量的单边散度

$d\left(x\right)=\left\{\begin{array}{c}\mathrm{div}(x_u,x_v),\mathrm{div}(x_u,x_v)<0\\ 0,\mathrm{otherwise}\end{array}\right.---\left(4\right)$

双边滤波器将高斯空间滤波与距离滤波相结合，可以起到既去除噪声又保护边缘的作用，将双边滤波结果作为数据项权重表示当一个区域内的像素点离中心点空间距离越远，那么权重越小；当一个区域内的像素点灰度值与中心点灰度值相差越大，那么权重越小。在此基础上，引入阻塞度这一概念，阻塞度表示了当前像素点的运动矢量的杂乱程度，当一个区域内的像素点阻塞度与中心点阻塞度比值越小，那么权重越小。总之加强了对图像细节的精确计算。

对光流算法中平滑项引入全局扩散系数，增强平滑项中对边缘的保留。全局扩散系数定义为：

$D\left(\right||\&dtri;I|\left|\right)=e^{-\mathrm{\&alpha;}{\left|\right|\&dtri;I\left|\right|}^{\mathrm{\&beta;}}}---\left(5\right)$

可以利用图像中的亮度信息来获取一扩散系数，扩散系数表明了图像中物体的连续程度，用它来对平滑项作加权处理使得图像中无纹理的连续区域平滑程度高，而边角这种非连续区域平滑程度低。总之进一步加强了弱纹理处的平滑程度。

这样，对两项系数改进后得到的误差式为：

$E_{\mathrm{CLG}-\mathrm{TV}}=\underset{\mathrm{\&Omega;}}{\&Integral;}[\left(\underset{\mathrm{region}}{\mathrm{\&Sigma;}}\mathrm{BFW}\left(x\right){(I_{x}u+I_{y}v+I_t)}^2\right)+\mathrm{\&lambda;}\left(\right||D\&CenterDot;\&dtri;u||+||D\&CenterDot;\&dtri;v|\left|\right)]---\left(6\right)$

该误差式的平滑项是TV-L1模型，其求解属于全变分正则化(total variationregularization)问题，直接求解是有困难的，现将其分为三部分：

$E_{\mathrm{CLG}-\mathrm{TV}-1}=\underset{\mathrm{\&Omega;}}{\&Integral;}[\left(\underset{\mathrm{region}}{\mathrm{\&Sigma;}}\mathrm{BFW}\left(x\right){(I_x\hat{u}+I_y\hat{v}+I_t)}^2\right)+\frac{1}{2\&CenterDot;\mathrm{\&theta;}}{(u-\hat{u})}^2+\frac{1}{2\&CenterDot;\mathrm{\&theta;}}{(v-\hat{v})}^2]---\left(7\right)$

$E_{\mathrm{TV}-u}=\underset{\mathrm{\&Omega;}}{\&Integral;}[\frac{1}{2\&CenterDot;\mathrm{\&theta;}}{(u-\hat{u})}^2+||D\&CenterDot;\&dtri;u|\left|\right]---\left(8\right)$

$E_{\mathrm{TV}-v}=\underset{\mathrm{\&Omega;}}{\&Integral;}[\frac{1}{2\&CenterDot;\mathrm{\&theta;}}{(v-\hat{v})}^2+||D\&CenterDot;\&dtri;v|\left|\right]---\left(9\right)$

使用交替求解的方法来最小化上面三个式子：

(6)式的求解是凸函数最优化问题，可以直接求其极值点作为最优解，作为未知量，将第一个式子分别对求偏导，并令其为零，整理化简可以得到方程组：

$\left(\begin{array}{cc}1+2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_x^2& 2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_x{I}_y\\ 2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_x{I}_y& 1+2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_y^2\end{array}\right)\&CenterDot;\left(\begin{array}{c}\hat{u}\\ \hat{v}\end{array}\right)=\left(\begin{array}{c}u-2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_x{I}_t\\ v-2\mathrm{\&theta;\&Sigma;BFW}\&CenterDot;I_y{I}_t\end{array}\right)---\left(10\right)$

可以通过求解以上线性方程组得到第一个式子的解

对于TV-u和TV-v式，对于这种形式的积分式应该构造欧拉--拉格朗日方程然后使用原始对偶算法来求解，TV-u的欧拉--拉格朗日方程为：

$-\mathrm{div}(D\&CenterDot;\frac{\&dtri;u}{\left|\right|\&dtri;u\left|\right|})+\frac{1}{\mathrm{\&lambda;}}\&CenterDot;(u-\hat{u})=0---\left(11\right)$

使用原始对偶算法，定义对偶向量为

$\stackrel{\&RightArrow;}{p_u}=\frac{\&dtri;u}{\left|\right|\&dtri;u\left|\right|}\&DoubleRightArrow;\stackrel{\&RightArrow;}{p_u}\&CenterDot;\left|\right|\&dtri;u\left|\right|-\&dtri;u=0,\left|\right|\stackrel{\&RightArrow;}{p_u}\left|\right|\&le;1---\left(12\right)$

将欧拉—拉格朗日方程改写为

$-\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\frac{1}{\mathrm{\&lambda;}}\&CenterDot;(u-\hat{u})=0\&DoubleRightArrow;u=\mathrm{\&lambda;}\&CenterDot;\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\hat{u}---\left(13\right)$

再将(13)式代入(12)式可得：

$\begin{array}{c}\stackrel{\&RightArrow;}{p_u}\&CenterDot;\left|\right|\&dtri;(\mathrm{\&lambda;}\&CenterDot;\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\hat{u})\left|\right|-\&dtri;(\mathrm{\&lambda;}\&CenterDot;\mathrm{div}(D\&CenterDot;{\stackrel{\&RightArrow;}{p}}_u)+\hat{u})=0\\ \&DoubleRightArrow;\stackrel{\&RightArrow;}{p_u}\&CenterDot;\left|\right|\&dtri;(\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\frac{\hat{u}}{\mathrm{\&lambda;}})\left|\right|-\&dtri;(\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\frac{\hat{u}}{\mathrm{\&lambda;}})=0\end{array}---\left(14\right)$

为求上式中结合方程形式可以使用固定点迭代的方法求解：

$\stackrel{\&RightArrow;}{{p_u}^{n+1}}=\frac{\stackrel{\&RightArrow;}{{p_u}^n}+\mathrm{\&tau;}\&CenterDot;\&dtri;(\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\frac{\hat{u}}{\mathrm{\&lambda;}})}{1+\mathrm{\&tau;}\&CenterDot;\left|\right|\&dtri;(\mathrm{\&lambda;}\&CenterDot;\mathrm{div}(D\&CenterDot;\stackrel{\&RightArrow;}{p_u})+\hat{u})\left|\right|}---\left(15\right)$

此外，求解光流场还需要用到两个策略：

1)由粗到精的多尺度计算，降低图像的分辨率，得到金字塔图像。由于当像素点位移大于滤波窗口的尺寸时，使用线性化的数据项直接计算光流会显著增加光流的误差，为了适用于计算较大的光流，一种解决方法是使用由粗到精的多尺度计算，较大的光流在低分辨率尺度变得很小，足以满足光流约束方程。本方法使用自动决定层数的CLG-TV金字塔光流法，原始图像是最底层图像，高一层的图像由比它低一层的图像经过特定的高斯滤波然后进行下采样得到，此处选取每层的图像的长和宽都是比它低一层的图像的四分之一。由粗到精的多尺度计算很好的适应了了计算较大的光流的要求，另外，将金字塔底层图像的分辨率从VGA级别降低至较低的分辨率，减少了计算量。

2)翘曲迭代(warp)方法

实际计算中，运动矢量的求解是经过多次翘曲得到的，将每次的计算结果作为右眼图像的自变量数组进行双线性差值翘曲，然后将左眼图像和翘曲后的右眼图像的运算结果更新上次的结果。

整个CLG-TV金字塔光流法的计算流程如图2所示。最后，将程序的步骤以伪代码的形式列在下面：

1)建立图像金字塔I1，I2，并计算各层梯度。

2)对各层I1计算全局扩散系数D。

3)每层进行如下步骤：

a)如果在最高层，将u,v,初始化为0。

b)如果不在最高层，将上层的u,v,上采样，作为本层的初始值。

c)for i＝1to warps(采用翘曲迭代方法)

i)对u,v中值滤波，减少异常值。

ii)用双线性插值对I₁,I₂,进行warp操作

iii)for k＝1to iterations

A.使用(10)式计算

B.使用(13)式更新u,v。

C.使用(15)式更新

本步骤在CLG-TV光流法的基础上，对光流算法中平滑项进行了改进，使用了各向异性的平滑权重，改善了平滑效果；对光流算法中数据项的梯度算子、窗口权重函数进行改进，加强了对图像细节的精确计算；使用由粗到精的多尺度计算，翘曲迭代方法，加强了对大位移运动的鲁棒性。求解光流场仅计算到第3层为止，目的是对最大视差的像素所在位置作一个初步的估计，通过很小的计算量准确地得到最大视差在整幅图像中的大致位置，以待之后进行精确计算。

步骤3、对步骤2计算得到的光流场进行处理，找出其中的最大值u_max以及最大值所在的位置(i_umax,j_umax)，在I1中的该位置截取原图像的一部分，并找到I2中经最大值移动后的该位置(i_umax+u_max·scale,j_umax)，在此处截取原图像的一部分。光流法的精度保证了截取的小图的计算准确度，同时，截取一小幅图片大大减少了数据量使得计算速度提高。

步骤4、使用CLG-TV光流法精确求解光流场，本步骤与步骤2的不同之处在于：本步骤使用的输入图像I1，I2是步骤3输出的截取图像，而不是原图像；为保证精度，本步骤使用的金字塔采样系数为0.75，并计算到第1层为止；本步骤的LK窗口宽度为11，而步骤2中为5。本层的输出为截取图像的光流。

步骤5、将步骤4的输出换算为第0层的值后取其水平运动矢量，作出直方图，从中计算最大值时，对直方图进行一阶差分，然后对一阶差分序列从最小值开始考察，当一阶差分大于一阈值时，将低于该阈值的运动矢量取出，计算其方差，并根据方差做出判断：

$u_{\max }=\left\{\begin{array}{c}\min \left(u_{\mathrm{\&psi;}}\right),D\left(u_{\mathrm{\&psi;}}\right)<\mathrm{\&alpha;}\\ \max \left(u_{\mathrm{\&psi;}}\right),\mathrm{otherwise}\end{array}\right.$

将此时得到的最大值与步骤2得到的最大值作算术加，得到的结果为整个系统的输出即最大视差。

与现有技术相比，本发明具有如下优点：

不同于传统的立体匹配方法，本发明的具体创新点在于采用了光流法来获取最大视差。由于光流法对弱纹理及重复纹理，噪声干扰等的影响较不敏感，可以得到较准确的最大视差值。在CLG-TV金字塔光流法的基础上，对光流算法中平滑项引入全局扩散系数，对光流算法中数据项的梯度算子、窗口权重函数进行改进，使其更加适用于计算双目图像只有水平位移的情况，对金字塔计算方法，分简单和复杂两次计算光流，有效减小计算量，加快了计算速度。经过多组测试序列实验，本发明可以获取的最大视差正确率可达90％。

附图说明

图1为本发明具体流程图。

图2为CLG-TV金字塔光流法具体流程图。

具体实施方式

本发明的最大视差获取方法基本流程如图1所示，具体包括以下步骤：

1)将双目图像序列输入计算机，首先对输入图像做预处理，如果源图像的分辨率较大可以先做下采样使图像缩小到VGA级别，将彩色图转为灰度图，然后对图像进行高斯平滑操作来减轻摄像头采集过程中产生的随机噪声。

2)使用CLG-TV光流法求解光流场，得到初步估计。

CLG-TV金字塔光流法的计算流程如图2所示，程序的伪代码如下：

I)建立图像金字塔I1，I2，并计算各层梯度。

II)对各层I1计算全局扩散系数D。

III)每层进行如下步骤：

A)如果在最高层，将u,v,初始化为0。

B)如果不在最高层，将上层的u,v上采样，作为本层的初始值。

C)for i＝1to warps(采用翘曲迭代方法)

i)对u,v中值滤波，减少异常值。

ii)用双线性插值对I₁,I₂,进行warp操作

iii)for k＝1to eq_iterations

a.使用(10)式计算

b.使用(13)式更新u,v。

c.使用(15)式更新

对于亮度对空间的梯度I_x和I_y，采用了Barron算子与Horn算子结合的方式。

Barron算子是一维的，其形式为1/12(180-8-1)；horn算子是三维的，它的计算方法如下式所示：

$I_x=\frac{1}{4}(I_{i+1,j,k}+I_{i+1,j+1,k}+I_{i+1,j,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i,j+1,k}+I_{i,j,k+1}+I_{i,j+1,k+1})$

$I_y=\frac{1}{4}(I_{i,j+1,k}+I_{i+1,j+1,k}+I_{i,j+1,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i+1,j,k}+I_{i,j,k+1}+I_{i+1,j,k+1})$

$I_t=\frac{1}{4}(I_{i,j,k+1}+I_{i+1,j,k+1}+I_{i,j+1,k+1}+I_{i+1,j+1,k+1})-\frac{1}{4}(I_{i,j,k}+I_{i+1,j,k}+I_{i,j+1,k}+I_{i+1,j+1,k})$

将Barron算子与Horn算子计算出的梯度分别赋予0.8和0.2的权重相加后得到最后的梯度，对于亮度对时间的梯度I_t，直接计算Horn算子的梯度作为它的值。

使用改进的双边滤波器来计算窗口权重函数，参数λ_s＝25,λ_i＝100,λ_o＝10,σ₁＝7,σ₂＝0.7,σ_d＝0.3,σ_e＝20。

对光流算法中平滑项引入全局扩散系数，参数α＝5，β＝0.5。

使用光流法求解光流场时金字塔采样系数为0.25，从最高层计算到第3层为止，LK窗口的宽度为5，warp次数为5，求解等式迭代次数为5。

3)对2)计算得到的光流场进行处理，找出其中的最大值u_max以及最大值所在的位置(i_umax,j_umax)，在I1中的该位置截取原图像的一部分，并找到I2中经最大值移动后的该位置(i_umax+u_max·scale,j_umax)，在两图的此处作为中点截取长宽各为原图像1/3一部分。

4)使用CLG-TV光流法精确求解光流场，本步骤中梯度算子、窗口权重函数以及全局扩散系数的参数与2)相同。金字塔采样系数为0.75，从最高层计算到第1层为止，LK窗口的宽度为11，warp次数为5，求解等式迭代次数为10。

5)将4)的输出换算为第0层的值后取其水平运动矢量，作出直方图，从中计算最大值时，对直方图进行一阶差分，然后对一阶差分序列从最小值开始考察，当一阶差分大于一阈值时，将低于该阈值的运动矢量取出，计算其方差，并根据方差做出判断：

$u_{\max }=\left\{\begin{array}{c}\min \left(u_{\mathrm{\&psi;}}\right),D\left(u_{\mathrm{\&psi;}}\right)<\mathrm{\&alpha;}\\ \max \left(u_{\mathrm{\&psi;}}\right),\mathrm{otherwise}\end{array}\right.$

将此时得到的最大值与步骤2得到的最大值作算术加，得到的结果为整个系统的输出即最大视差。

Claims (2)

1.基于光流法的双目图像最大视差获取方法，其特征在于：本方法采用了光流法来获取最大视差；通过对左右摄像机采集到的图像进行一系列处理，得到图像上各点的运动矢量，对运动矢量进行分析判断最后得到最大视差，以便于进行下一步处理；由于光流法对弱纹理及重复纹理，噪声干扰的影响较不敏感，得到较准确的最大视差值；

具体而言，该方法的具体实施过程如下，

1)由粗到精的多尺度计算：按照一定的采样标准，把原始图像分为多层，形成金字塔形式的图像结构，原始图像位于金字塔底部，而金字塔顶部则是图像分层后最低分辨率的近似；这样的多分辨率思想在使用光流分析大位移运动时特别适用，同时使用较小的取样因子使计算量减少，加快了计算速度；

2)CLG-TV光流法求解光流场；使用CLG-TV光流法处理左右图像获取光流运动矢量；

3)对光流算法中平滑项引入全局扩散系数，有效地在不连续处进行精确计算；

4)对光流算法中数据项的梯度算子以及窗口权重函数的改进；使用一种更加适用于计算最大视差的梯度算子的计算方法，并通过改进的双边滤波器来计算窗口权重函数；

5)考虑到最大视差唯一性的特点，分简单和复杂两次计算光流，有效减小计算量，加快了计算速度；

实施流程如下，

步骤1、使用左右摄像机采集到的双目图像序列输入计算机，首先对输入图像做预处理，如果源图像的分辨率较大先做下采样使图像缩小到VGA级别以减少计算量，由于本方法不使用图像的彩色信息，首先将彩色图转为灰度图，公式如下：

I＝0.299R+0.587G+0.114B 

然后对图像进行高斯平滑操作来消除摄像头采集过程中产生的随机噪声；

步骤2、使用CLG-TV光流法求解光流场，得到初步估计；

光流用来描述物体在观测到的图像中的二维运动，光流法运动估计的计算结果为图像中像素的运动矢量；理想情况下，对双目图像计算光流场将得到图 像中像素点的水平运动矢量，这些水平运动矢量就等于立体匹配算法得到的各像素点的视差，光流法计算的运动矢量包含有视差信息，应用于双目图像最大视差的获取；光流法的基本假设将亮度约束和平滑约束相结合，既保证了运动矢量计算的准确性，又使得各区域内运动矢量能够平滑过渡，是一种高精度的运动估计方法；实际图像存在复杂多变的情况，同一景物在左右视点下的成像会有很大差异，会对获得精确的视差值造成困难；利用光流法求解最大视差的优势就在于光流法的精度能够达到亚像素级，完全可以胜任于最大视差的计算；光流法对传统方法中比较难以处理的深度不连续、弱纹理及重复纹理，噪声干扰等问题的鲁棒性强；

CLG-TV光流法是结合全局与局部的光流法，基于全局的光流法对小运动精确度高，基于局部的光流法对大运动有更好的精度，将二者结合起来增强鲁棒性，CLG-TV光流法使用的误差计算式将LK方法和HS方法的误差式结合起来；另外，HS和LK的误差计算式都采用了L2型范数；L2型范数有一些缺点：对噪声敏感度高，过度扩张，选择性差，以及不保留运动边界，为克服L2模型的缺点提出了采用L1型范数的方法将误差式改写，其形式为：

设I(x,y,t)是图像点(x,y)在时刻t的亮度，u(x,y)和v(x,y)是该点光流的x和y分量，点在t+δt时运动到(x+δx,y+δy)时，亮度保持不变，即

I(x+uδt,y+vδt,t+δt)＝I(x,y,t)

亮度随x、y、t光滑变化，将上式的左边用泰勒级数一阶展开：

其中e是关于δx、δy、δt的二阶和二阶以上的项；上式两边的I(x,y,t)相互抵消，两边除以δt，并取极限δt→0，得到这就是亮度不变性约束；

(1)式中左式是LK方法的数据项，其本质是一个区域上的所有像素点满足亮度不变性约束I_xu+I_yv+I_t＝0；其中亮度对空间的梯度I_x和I_y的物理意义是亮 度在水平、竖直方向的变化率，亮度对时间的梯度I_t的物理意义是亮度随时间的变化率，u和v代表像素点分别在水平，竖直方向的变化速度；

(1)式中的右式是L1型的平滑项，由于实际场景中的物体都是连续，其运动矢量理应是大面积相同或相近的，因此设计平滑项这一项来满足这一特点； 表示水平运动速度的空间梯度，反映了水平运动速度u在空间上的连续程度，L1型范数是指模的大小，其值越小连续程度越高；竖直方向同理；整个平滑项的物理意义是指求得的整个运动场应满足所有像素点的 之和最小，即整个运动场的大部分点是足够连续，足够平滑的；这样就克服了深度不连续、弱纹理及重复纹理和噪声干扰造成的影响；

对于亮度对空间的梯度I_x和I_y，采用了Barron算子与horn算子结合的方式；Barron算子是一维的，其形式为1/12(1 8 0 -8 -1)；horn算子是三维的，其计算方法如下式所示：

Barron算子适合于计算空间梯度，horn算子适合于计算时间梯度；对于亮度对空间的梯度I_x和I_y，将Barron算子与horn算子计算出的梯度进行加权作为最后梯度，亮度对时间的梯度I_t，将horn算子算出的梯度作为最后的梯度；

通过改进的双边滤波器来计算窗口权重函数，增强数据项中对边缘的保留，在双边滤波器的基础上增加阻塞度一项，使得权重正比于空间距离，亮度差别，阻塞度差别：

其中，阻塞度o(x)定义为

其中d(x)为速度向量的单边散度

双边滤波器将高斯空间滤波与距离滤波相结合，起到既去除噪声又保护边缘的作用，将双边滤波结果作为数据项权重表示当一个区域内的像素点离中心点空间距离越远，那么权重越小；当一个区域内的像素点灰度值与中心点灰度值相差越大，那么权重越小；在此基础上，引入阻塞度这一概念，阻塞度表示了当前像素点的运动矢量的杂乱程度，当一个区域内的像素点阻塞度与中心点阻塞度比值越小，那么权重越小；总之加强了对图像细节的精确计算；

对光流算法中平滑项引入全局扩散系数，增强平滑项中对边缘的保留；全局扩散系数定义为：

利用图像中的亮度信息来获取一扩散系数，扩散系数表明了图像中物体的连续程度，用它来对平滑项作加权处理使得图像中无纹理的连续区域平滑程度高，而边角这种非连续区域平滑程度低；总之进一步加强了弱纹理处的平滑程度；

这样，对两项系数改进后得到的误差式为：

该误差式的平滑项是TV-L1模型，其求解属于全变分正则化问题，直接求解是有困难的，现将其分为三部分：

使用交替求解的方法来最小化上面三个式子：

(6)式的求解是凸函数最优化问题，直接求其极值点作为最优解，作为未知量，将第一个式子分别对求偏导，并令其为零，整理化简得到方程组：

通过求解以上线性方程组得到第一个式子的解

对于TV-u和TV-v式，对于这种形式的积分式应该构造欧拉--拉格朗日方程然后使用原始对偶算法来求解，TV-u的欧拉--拉格朗日方程为：

使用原始对偶算法，定义对偶向量为

将欧拉—拉格朗日方程改写为

再将(13)式代入(12)式可得：

为求上式中结合方程形式可以使用固定点迭代的方法求解：

此外，求解光流场还需要用到两个策略：

1)由粗到精的多尺度计算，降低图像的分辨率，得到金字塔图像；由于当像素点位移大于滤波窗口的尺寸时，使用线性化的数据项直接计算光流会显著增加光流的误差，为了适用于计算较大的光流，一种解决方法是使用由粗到精的多尺度计算，较大的光流在低分辨率尺度变得很小，足以满足光流约束方程；本方法使用自动决定层数的CLG-TV金字塔光流法，原始图像是最底层图像，高一层的图像由比它低一层的图像经过特定的高斯滤波然后进行下采样得到，此处选取每层的图像的长和宽都是比它低一层的图像的四分之一；由粗到精的多尺度计算很好的适应了了计算较大的光流的要求，另外，将金字塔底层图像的分辨率从VGA级别降低至较低的分辨率，减少了计算量；

2)翘曲迭代方法

实际计算中，运动矢量的求解是经过多次翘曲得到的，将每次的计算结果作为右眼图像的自变量数组进行双线性差值翘曲，然后将左眼图像和翘曲后的右眼图像的运算结果更新上次的结果；

整个CLG-TV金字塔光流法的计算流程；最后，将程序的步骤以伪代码的 形式列在下面：

1)建立图像金字塔I1，I2，并计算各层梯度；

2)对各层I1计算全局扩散系数D；

3)每层进行如下步骤：

a)如果在最高层，将u，初始化为0；

b)如果不在最高层，将上层的u，上采样，作为本层的初始值；

c)for i＝1 to warps(采用翘曲迭代方法)

i)对u,v中值滤波，减少异常值；

ii)用双线性插值对I₁，I₂，进行warp操作

iii)for k＝1 to iterations

A.使用(10)式计算

B.使用(13)式更新u,v；

C.使用(15)式更新

本步骤在CLG-TV光流法的基础上，对光流算法中平滑项进行了改进，使用了各向异性的平滑权重，改善了平滑效果；对光流算法中数据项的梯度算子、窗口权重函数进行改进，加强了对图像细节的精确计算；使用由粗到精的多尺度计算，翘曲迭代方法，加强了对大位移运动的鲁棒性；求解光流场仅计算到第3层为止，目的是对最大视差的像素所在位置作一个初步的估计，通过很小的计算量准确地得到最大视差在整幅图像中的大致位置，以待之后进行精确计算；

步骤3、对步骤2计算得到的光流场进行处理，找出其中的最大值u_max以及最大值所在的位置(i_umax,j_umax)，在I1中的该位置截取原图像的一部分，并找到I2中经最大值移动后的该位置(i_umax+u_max·scale,j_umax)，在此处截取原图像的一部分；光流法的精度保证了截取的小图的计算准确度，同时，截取一小幅图片大大减少了数据量使得计算速度提高；

步骤4、使用CLG-TV光流法精确求解光流场，本步骤与步骤2的不同之处在于：本步骤使用的输入图像I1，I2是步骤3输出的截取图像，而不是原图像；为保证精度，本步骤使用的金字塔采样系数为0.75，并计算到第1层为止；本步骤的LK窗口宽度为11，而步骤2中为5；本层的输出为截取图像的光流；

步骤5、将步骤4的输出换算为第0层的值后取其水平运动矢量，作出直方图，从中计算最大值时，对直方图进行一阶差分，然后对一阶差分序列从最小值开始考察，当一阶差分大于一阈值时，将低于该阈值的运动矢量取出，计算其方差，并根据方差做出判断：

将此时得到的最大值与步骤2得到的最大值作算术加，得到的结果为整个系统的输出即最大视差。

2.根据权利要求1所述的基于光流法的双目图像最大视差获取方法，其特征在于：

最大视差获取方法基本流程，具体包括以下步骤：

1)将双目图像序列输入计算机，首先对输入图像做预处理，如果源图像的分辨率较大先做下采样使图像缩小到VGA级别，将彩色图转为灰度图，然后对图像进行高斯平滑操作来减轻摄像头采集过程中产生的随机噪声；

2)使用CLG-TV光流法求解光流场，得到初步估计；

CLG-TV金字塔光流法的计算流程，程序的伪代码如下：

I)建立图像金字塔I1，I2，并计算各层梯度；

II)对各层I1计算全局扩散系数D；

III)每层进行如下步骤：

A)如果在最高层，将u,v,初始化为0；

B)如果不在最高层，将上层的u,v,上采样，作为本层的初始值；

C)for i＝1to warps(采用翘曲迭代方法)

i)对u,v中值滤波，减少异常值；

ii)用双线性插值对I₁,I₂,进行warp操作

iii)for k＝1to eq_iterations

a.使用(10)式计算

b.使用(13)式更新u,v；

c.使用(15)式更新

对于亮度对空间的梯度I_x和I_y，采用了Barron算子与Horn算子结合的方式；

Barron算子是一维的，其形式为1/12(180-8-1)；horn算子是三维的，其计算方法如下式所示：

将Barron算子与Horn算子计算出的梯度分别赋予0.8和0.2的权重相加后得到最后的梯度，对于亮度对时间的梯度I_t，直接计算Horn算子的梯度作为其值；

使用改进的双边滤波器来计算窗口权重函数，参数λ_s＝25,λ_i＝100,λ_o＝10,σ₁＝7,σ₂＝0.7,σ_d＝0.3,σ_e＝20；

对光流算法中平滑项引入全局扩散系数，参数α＝5，β＝0.5；

使用光流法求解光流场时金字塔采样系数为0.25，从最高层计算到第3层为止，LK窗口的宽度为5，warp次数为5，求解等式迭代次数为5；

3)对2)计算得到的光流场进行处理，找出其中的最大值u_max以及最大值所在的位置(i_umax,j_umax)，在I1中的该位置截取原图像的一部分，并找到I2中经最大值移动后的该位置(i_umax+u_max·scale,j_umax)，在两图的此处作为中点截取长宽各为原图像1/3一部分；

4)使用CLG-TV光流法精确求解光流场，本步骤中梯度算子、窗口权重函数以及全局扩散系数的参数与2)相同；金字塔采样系数为0.75，从最高层计算到第1层为止，LK窗口的宽度为11，warp次数为5，求解等式迭代次数为10；

5)将4)的输出换算为第0层的值后取其水平运动矢量，输出直方图，从中计算最大值时，对直方图进行一阶差分，然后对一阶差分序列从最小值开始考察，当一阶差分大于一阈值时，将低于该阈值的运动矢量取出，计算其方差，并根据方差做出判断：

将此时得到的最大值与步骤2得到的最大值作算术加，得到的结果为整个系统的输出即最大视差。