, 王爱侠1,马帅21. 东北大学 信息科学与工程学院，辽宁 沈阳 110819;
2. 辽宁大学 信息学院，辽宁 沈阳 110036
收稿日期: 2015-04-07
基金项目: 辽宁省教育厅科学研究项目(L2012003); 沈阳市科技局项目(F12277181).
作者简介: 李晶皎(1964-),女,辽宁沈阳人，东北大学教授，博士生导师.
摘要: 针对立体匹配时视差不连续区、倾斜平面及非前向平行平面误匹配较高的问题，提出了一种基于改进Patchmatch及切片采样粒子置信度的立体匹配算法.定义了具有边缘特性的Patchmatch相似性函数，并建立基于Patchmatch的非前向平行平面视差平面估计模型.利用粒子置信度传播代替原有的最近邻搜索，使用较少的粒子近似目标分布，并采用切片采样马尔可夫链蒙特卡罗方法解决传播过程中粒子重采样更新问题.Middlebury图像数据集测试表明，该算法能够降低视差不连续区域的误匹配，有效地提高了倾斜平面及非前向平行平面图像的匹配精度.
关键词: 立体匹配Patchmatch粒子置信度传播切片采样马尔可夫链蒙特卡罗
Stereo Matching Algorithm Based on Improved Patchmatch and Slice Sampling Particle Belief Propagation
LI Jing-jiao1, MA Li1,2
, WANG Ai-xia1,MA Shuai21.School of Information Science & Engineering, Northeastern University, Shenyang 110819, China;
2. School of Information, Liaoning University, Shenyang 110036, China.
Corresponding author: MA Li, E-mail: mali@lnu.edu.cn
Abstract: The high erroneous results of the stereo matching occur at the following three cases where there are depth discontinuity region, the slanted surface or the non-fronto-parallel surface. A stereo matching algorithm was proposed based on the improved Patchmatch and slice sampling particle belief propagation. An edge-preserving similarity function of the Patchmatch was defined. Then, a model of the depth estimation for the non-fronto-parallel surface was introduced. The nearest neighbor search was replaced with the particle belief propagation, and the target distribution was approximated with a finite set of particles. At the same time, the sampled particles from the belief distribution was typically done by using slice sampling Markov chain Monte Carlo method to solve the particle update problem. The experiments on the Middlebury indicate that the mismatching at the depth discontinuity region can be reduced, and the match accuracy for the slanted surface and the non-fronto-parallel surface can be improved.
Key words: stereo matchingPatchmatchparticle belief propagationslice samplingMarkov chain Monte Carlo 立体匹配是立体视觉的关键技术之一.立体匹配过程中，多采用前向平行表面模型的假设情况来实现[1]，而真实场景中多数情况不符合该假设，如视差不连续区域、倾斜平面和非前向平行平面，这种情况下误匹配率会大幅提高.
为了提高非前向平面情况下的匹配精度，有学者采用基于图像分割的立体匹配算法[2, 3]，通过将具有不同视差的像素分割为不同区域，估计分割区域的匹配视差平面进行匹配.还有学者采用Patchmatch随机一致性算法[4, 5]，快速地计算两幅图像中块之间的近似稠密的最近邻估计，通过不断迭代逼近最优三维表面得到深度图.文献[6]提出了一种将Patchmatch与粒子置信度传播算法结合的方法，解决了Patchmatch算法在大面积低纹理区域误匹配率高的问题，但该算法在初始化随机匹配区域时，不能充分体现区域的边缘特性，且粒子置信度传播过程中的粒子重采样必须选取合适的建议分布，增加了算法难度，直接影响匹配精度.本文提出了一种基于改进Patchmatch及切片采样粒子置信度的立体匹配算法，在Patchmatch算法中采用具有边缘特性的相似性函数，应用切片采样[7]马尔可夫链蒙特卡罗方法解决传播过程中粒子重采样更新问题.
1 基于Patchmatch与粒子置信度传播的立体匹配算法基于Patchmatch和粒子置信度传播的视差估计算法[6]利用最邻近域算法随机初始化三维平面，采用粒子置信度传播实现迭代传播，最终得到全局最优化匹配结果生成稠密视差图.
1.1 基于Patchmatch的非前向平行视差平面估计模型将图像中任意像素定义为三维空间中的一个三维平面上的一点，若其邻近的像素与它具有很大的相似性，则认为它们所在的平面也是近似相同的，这样更逼近真实场景中的倾斜平面及曲面.在Patchmatch算法随机初始化过程中，对于图像中每一个像素p存在一个三维随机平面，用3参数向量表示为fp=(xp,yp,zp)，则像素p的视差值用平面参数来表示为
令三维随机平面的单位法向量为n=(nx,ny,nz)，则式(1)中参数可以表示为
1.2 粒子置信度传播粒子置信度传播算法[8]是将粒子滤波引入置信度传播算法中，初始化定义一系列粒子P来表示置信度中的后验边缘估计，通过抽样或马尔可夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)采样选择合适的建议分布从后验分布中采样，以构建一条收敛于后验分布的马氏链.
粒子BP算法应用MRF模型，采用最大积粒子BP算法[9]，对每个节点s进行采样得到含有q个粒子的粒子集Ps={xs(1)…xs(q)}，任意节点处n次迭代后，节点s最大积置信度估计为bsn(xs(i))，与其对应的对数置信度为Bsn(xs(i))=-lg(bsn·(xs(i))).节点s置信度定义为
从节点t到节点s的消息Mt→sn(xs)为
对于立体匹配算法，式(3)中ψs(us)为MRF模型中的数据项，用来计算局部势能，节点s表示为us=[ds]=[asxt+bsyt+cs]，数据项ψs(us)定义为块内待匹配像素s的聚合匹配代价，表示为
其中:函数ρ(s,t)计算像素s和t的像素信息相似度;w(s,t)定义为像素差异权重.
ψst(us,ut)为平滑项，代表{s,t}节点间的势能函数，求解在(x,y,disparity)视差空间中的两个局部平面差异.令节点s的平面法向量为ns=orth([xs,ys,zs]T)，点us为这个平面上的一点，定义平滑项计算为
式中:ωst为平滑项的数据依赖项;β为用来约束平滑项的常数权重.
在迭代传播过程中，利用Metropolis-Hastings MCMC(MH MCMC)采样法[9]对节点邻域进行重采样，在当前粒子集Ps及已传播的邻域候选粒子集Cs中选取能够使节点处数据项ψs最小的K个“最佳”粒子对当前粒子集进行更新.经n次消息和置信度的迭代计算，选择具有最小置信度的粒子集Ps作为立体匹配中的视差估计：
2 算法改进2.1 Patchmatch算法改进应用Patchmatch数据项计算匹配代价时通常采用尺寸较大的像素块，这样会丢失像素块的边缘特性.本文在Patchmatch相似度度量中，将像素色度空间和几何空间的邻近性引入权值计算，采用双边权重方法计算像素块之间的匹配代价，以此来达到还原输入视图细节，保持像素块的边缘特性，减少视差不连续区域的误匹配.像素匹配代价即数据项可以表述为
其中:W为权重归一化参数;w′(s,t)为双边权重函数：
式中：‖s－t‖表示像素s和t的空间距离；‖Is－It‖为采用CIELab色彩空间中3个分量的欧式距离；σs和σr是用来控制空间和颜色差异的参数.
2.2 切片采样粒子置信度传播为了避免对建议分布的依赖，采用切片采样(slice sampling)粒子置信度传播来完成对节点邻域的重采样，选择合适的粒子集Ps来更新节点s的消息和置信度，使匹配能量最小.
将每个节点处置信度作为目标密度函数，节点s具有粒子集Ps={xs(i)}i=1,…,k，置信度传播的第n次迭代时，计算节点s的最大积置信度估计为bsn(xs(i))，对应对数置信度为Bsn(xs(i)).对节点s每个粒子xs(i)的重采样从分布q(xs(i))中采样得到，构成MCMC采样链{xs(i)〈m〉}m=1,…,M.其中节点s处的置信度q(x)可表示为
设辅助变量h∈R，则目标分布扩展为
通过对辅助变量h进行均匀采样形成切片，在范围h内每间隔A进行一次新的采样.
这里UI为间隔I内的均匀分布，且间隔A={x;q(x)≥h〈m〉}，图 1为切片采样示例.
图1(Figure 1)
 | 图1 切片采样 Fig. 1 Slice sampling |
令fl(x)=mt(l)→s(xs(i))，将q(x)分解，并通过L个辅助变量h1,…,hL得到采样，采样间隔为
即对于节点s的第i个粒子来说，在第m次MCMC迭代中以采样间隔A进行均匀采样，得到粒子xs(i)〈m〉(在下面论述中上标(i)〈m〉省略).经对数空间变换，节点s的对数置信度表示为
其中:F0(xs)=ψs(xs);Fl(xs)=Mt(l)→s(xs);t(l)表示节点s的第l个邻域t;|Ns|为最大邻域个数，在此采用节点四邻域关联则l={1，…，4}.因此可以将采样间隔分解为
为负对数空间的均匀采样值.通过推导而知，依据马尔可夫模型的势能函数ψs(xs)和ψst(xs,xt)，采样间隔A仅由数据项势能函数ψs和平滑项势能函数ψst定义的给定间隔Aψs(
)和Axtψst()决定：
由式(11)可知
本算法中，对于式(5)数据项难以进行不等式求解区间，则在确定采样间隔Aψs(
)时，xs像素选择范围定义为像素块内所有节点，而对于xs的视差选择范围定义为ds∈[0,disp max－1]. Axtψst()则通过式(6)及式(17)可以表示为
β和ωst为非负用户自定义常数.ns和us表示为ns=[nxs,nys,nzs],us=[xs,ys,zs]，令
′l=
(l-Bt(xt)+Ms→t (xt))，经矩阵运算得到:
在采样间隔计算时，分别对xs在(x,y,disparity)视差空间各个维度进行求解.因此，令Δx=|xs－xt|,Δy=|ys－yt|,Δz=|zs－zt|，因为在重采样过程中，要求重采样得到的粒子xt应同样是处于像素块内，并且其视差应小于最大视差搜索范围，所以这里0≤Δx≤patchsize,0≤Δy≤patchsize,0≤Δz≤dip max，推导得到xs在x维度采样间隔为
已知Δy和Δz取值范围，可求得X的取值范围，进而得到x维度采样间隔.同理可得xs在y,z维度的采样间隔.
为了保证采样空间足够大，这里在计算X,Y,Z时，其中的Δx,Δy,Δz均取最小值.因此可以确定在不同维度下xs的采样间隔.在进行采样时，则可在各自维度的采样间隔内依次在x维度、y维度、z维度抽取样本来构造马氏链.
3 算法描述将改进Patchmatch算法和切片采样粒子置信度算法相结合，视差估计算法分为初始化、传播、局部重采样、视差图后处理4个步骤.其中前3个步骤将视差估计转化为视图中节点的最佳粒子集及置信度的最优求解过程，算法具体描述：
输入：左右视图，Ps为视图中节点s的粒子集，K为粒子数，N为置信度传播迭代次数，M为切片MCMC重采样迭代次数.
输出：粒子及置信度最优解(Ps,Bs).
1) 对于视图中所有节点s进行初始化，每个节点s 进行K次随机采样，得到初始粒子集Ps={xs(i)}i=1,…,K；初始化消息及置信度mt→s0(xs)=0,Bs0(xs(i))=0；
2) for n=1 to N
3) for 所有节点s的i=1,…,K个粒子;
4) 初始化采样链xs(i)〈0〉←xs(i);
5) 按照传播策略函数对T(s)所有节点s的邻域N(s)进行均匀采样，得到候选粒子集 Cs=∪ {Pt|t∈N(s)∩T(s)};
6) Cs=Ps;
7) for m=1 to M
8) 对所有Cs中的粒子xs(xs∈Cs)在l=1,…,|Ns|邻域进行均匀采样hl~U[0,1](h)，采样得到
l=Fl(xs(i)〈m－1〉)－lg(hl);9) 计算采样间隔A(i)〈m〉;
10) 采样得到xs(i)〈m〉~UA(i)〈m〉(x);
11) 计算置信度.
12) 满足Fl(
s(i)〈m〉) < l，则接受s(i)〈m〉，即xs(i)〈m〉←s(i)〈m〉;13) 更新粒子集Ps;
14) 对视图中所有节点s，返回粒子集Ps和置信度Bs，满足
.初始化及传播过程中在对节点s进行粒子随机采样时，可以对每个像素完成空间传播、视角传播及平面优化，对三维随机平面参数及视差进行随机选取，以得到近似平面.并通过视差图后处理步骤完成对遮挡区域的处理.
4 实验结果与分析本文算法在Linux环境下采用C++编程语言实现，实验采用Middlebury大学立体视觉算法评估数据库提供的测试图像完成.本算法中像素块尺寸为20×20像素，视图中每个节点选取粒子数K=5，BP迭代次数为n=10，MCMC迭代次数为m=10，双边权重中空间和颜色差异参数σs=15和σr=15，平滑项势能权重β=7.5.
本文对采用MH MCMC重采样和切片重采样的粒子置信度传播视差估计算法进行了评估.这里对MH MCMC重采样粒子置信度传播的建议分布采用高斯分布N(0,σ)，则在对粒子集Ps进行粒子重采样时得到粒子集为
选取MCMC重采样不同迭代次数情况下，生成的稠密视差图与真实视差图的均方根误差(RMSE)来评估不同重采样方法对视差估计算法性能的影响.图 2为不同迭代次数情况下，对存在倾斜平面及非前向平行平面的测试图Flowerpot和Wood1左视图的视差图性能估计结果.
图2(Figure 2)
 | 图2 切片采样MCMC与MH MCMC视差估计性能比较 Fig. 2 Comparison of the disparities estimation produced by slice sampling MCMC and MH MCMC |
应用Middlebury数据库中经典测试图Tsukuba，Venus，Teddy和Cones对Patchmatch[4]，PM-Huber[10]，PMBP[6]，PMF[5]和本文算法进行性能分析，性能评估结果见表 1.由表 1可见，本文在具有非前向平行表面的测试图的视差不连续区域误匹配、总体匹配优于其他算法，具有较好的边缘特性.而在纹理稀疏的测试图中算法对纹理较为敏感，所以误匹配要高.图 3为本算法对4幅测试图最终视差图的评估结果.
表1(Table 1)
 表1 本文算法与Middlebury平台基于Patchmatch算法的性能评估Table 1 Evaluation results of proposed algorithm and Patchmatch algorithm of Middlebury datasets
| 表1 本文算法与Middlebury平台基于Patchmatch算法的性能评估 Table 1 Evaluation results of proposed algorithm and Patchmatch algorithm of Middlebury datasets | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
图3(Figure 3)
 | 图3 本算法对Middlebury数据集测试结果 Fig. 3 Qualitative results of the proposed algorithm on Middlebury datasets |
由于切片采样粒子置信度传播过程中需要在3个维度对粒子的切片采样范围进行计算，算法的执行时间有所增加，其执行效率有待提升.
5 结论1) 针对立体匹配时视差不连续区、倾斜平面及非前向平行平面误匹配较高问题，应用具有边缘特性的Patchmatch相似性函数，本文建立基于Patchmatch的非前向平行平面视差平面估计模型.
2) 采用切片采样马尔可夫链蒙特卡罗方法解决粒子置信度传播过程中粒子重采样更新问题，降低了算法复杂性，提高了匹配精度.
3) Middlebury图像数据集测试表明，该算法能够降低视差不连续区域的误匹配，有效提高了倾斜平面及非前向平行平面图像的匹配精度.
参考文献
| [1] | Scharstein D,Szeliski R.A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[J].International Journal of Computer Vision,2002,47 (1/2/3):7-42.( 1) |
| [2] | Wang Z F,Zheng Z G.A region based stereo matching algorithm using cooperative optimization [C]// IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Anchorage,2008:1-8.(1) |
| [3] | Damjanovic S,Heijden F,Spreeuwers L J.Local stereo matching using adaptive local segmentation[J].International Scholarly Research Notices Machine Vision,2012(2012):1-15.(1) |
| [4] | Bleyer M,Rhemann C,Rother C.Patchmatch stereo-stereo matching with slanted support windows[J].British Machine Vision Conference,2011,14:1-11.(2) |
| [5] | Lu J,Yang H,Min D,et al.Patch match filter:efficient edge-aware filtering meets randomized search for fast correspondence field estimation[J].IEEE Computer Society Conference on Computer Vision and Pattern Recognition,2013,9(4):1854-1861.(2) |
| [6] | Besse F,Rother C,Fitzgibbon A,et al.PMBP:patch match belief propagation for correspondence field estimation[J].International Journal of Computer Vision,2014,110(1):2-13.(3) |
| [7] | Muller O,Yang M Y,Rosenhahn B.Slice sampling particle belief propagation[C]// IEEE International Conference on Computer Vision.Sydney,2013:1129-1136.(1) |
| [8] | Trinh H,Mcallester D.Particle-based belief propagation for structure from motion and dense stereo vision with unknown camera constraints[M].Berlin:Springer,2008:16-28.(1) |
| [9] | Kothapa R.Max-product particle belief propagation[D].Rhode Island:Brown University,2011.(2) |
| [10] | Heise P,Klose S,Jensen B,et al.PM-Huber:patch match with huber regularization for stereo matching[C]// IEEE International Conference on Computer Vision.Sydney,2013:2360-2367.(1) |
本文献在全文中的定位：
... 技术之一.立体匹配过程中，多采用前向平行表面模型的假设情况来实现[1]，而真实场景中多数情况不符合该假设，如视差不连续区域、倾斜平面 ...
1
本文献在全文中的定位：
... 前向平面情况下的匹配精度，有学者采用基于图像分割的立体匹配算法[2, 3]，通过将具有不同视差的像素分割 ...
1
本文献在全文中的定位：
... 采用基于图像分割的立体匹配算法[2, 3]，通过将具有不同视差的像素分割为不同区域，估计分割区域的匹配视 ...
2
本文献在全文中的定位：
... 割区域的匹配视差平面进行匹配.还有学者采用Patchmatch随机一致性算法[4, 5]，快速地计算两幅图像中块之间的 ...
... 应用Middlebury数据库中经典测试图Tsukuba，Venus，Teddy和Cones对Patchmatch[4]，PM-Huber[10]，PMB ...
2
本文献在全文中的定位：
... 有学者采用Patchmatch随机一致性算法[4, 5]，快速地计算两幅图像中块之间的近似稠密的最近邻估计，通过不断迭 ...
... PMBP[6]，PMF[5]和本文算法进行性能分析，性能评估结果见 ...
3
本文献在全文中的定位：
... 似稠密的最近邻估计，通过不断迭代逼近最优三维表面得到深度图.文献[6]提出了一种将Patchmatch与粒子置信度传播算法结合的方法，解决了Patchmat ...
... 基于Patchmatch和粒子置信度传播的视差估计算法[6]利用最邻近域算法随机初始化三维平面，采用粒子置信度传播实现迭代 ...
... PM-Huber[10]，PMBP[6]，PMF[5]和本文算法进行性能 ...
1
本文献在全文中的定位：
... 算法，在Patchmatch算法中采用具有边缘特性的相似性函数，应用切片采样[7]马尔可夫链蒙特卡罗方法解决传播过程中粒子重采样更新问题. ...
1
本文献在全文中的定位：
... 粒子置信度传播算法[8]是将粒子滤波引入置信度传播算法中，初始化定义一系列粒子 ...
2
本文献在全文中的定位：
... 粒子BP算法应用MRF模型，采用最大积粒子BP算法[9]，对每个节点s进行采样得到含有q个粒子的粒子集 ...
... 在迭代传播过程中，利用Metropolis-Hastings MCMC(MH MCMC)采样法[9]对节点邻域进行重采样，在当前粒子集Ps及已传播的邻域 ...
1
本文献在全文中的定位：
... 和Cones对Patchmatch[4]，PM-Huber[10]，PMBP[6]，PMF ...
