立体匹配是立体视觉的关键技术之一.立体匹配过程中，多采用前向平行表面模型的假设情况来实现^[1]，而真实场景中多数情况不符合该假设，如视差不连续区域、倾斜平面和非前向平行平面，这种情况下误匹配率会大幅提高.

为了提高非前向平面情况下的匹配精度，有学者采用基于图像分割的立体匹配算法^{[2, 3]}，通过将具有不同视差的像素分割为不同区域，估计分割区域的匹配视差平面进行匹配.还有学者采用Patchmatch随机一致性算法^{[4, 5]}，快速地计算两幅图像中块之间的近似稠密的最近邻估计，通过不断迭代逼近最优三维表面得到深度图.文献[6]提出了一种将Patchmatch与粒子置信度传播算法结合的方法，解决了Patchmatch算法在大面积低纹理区域误匹配率高的问题，但该算法在初始化随机匹配区域时，不能充分体现区域的边缘特性，且粒子置信度传播过程中的粒子重采样必须选取合适的建议分布，增加了算法难度，直接影响匹配精度.本文提出了一种基于改进Patchmatch及切片采样粒子置信度的立体匹配算法，在Patchmatch算法中采用具有边缘特性的相似性函数，应用切片采样^[7]马尔可夫链蒙特卡罗方法解决传播过程中粒子重采样更新问题.

基于Patchmatch和粒子置信度传播的视差估计算法^[6]利用最邻近域算法随机初始化三维平面，采用粒子置信度传播实现迭代传播，最终得到全局最优化匹配结果生成稠密视差图.

将图像中任意像素定义为三维空间中的一个三维平面上的一点，若其邻近的像素与它具有很大的相似性，则认为它们所在的平面也是近似相同的，这样更逼近真实场景中的倾斜平面及曲面.在Patchmatch算法随机初始化过程中，对于图像中每一个像素p存在一个三维随机平面，用3参数向量表示为f_p=(x_p,y_p,z_p)，则像素p的视差值用平面参数来表示为

令三维随机平面的单位法向量为n=(n_x,n_y,n_z)，则式(1)中参数可以表示为

粒子置信度传播算法^[8]是将粒子滤波引入置信度传播算法中，初始化定义一系列粒子P来表示置信度中的后验边缘估计，通过抽样或马尔可夫链蒙特卡罗(Markov chain Monte Carlo，MCMC)采样选择合适的建议分布从后验分布中采样，以构建一条收敛于后验分布的马氏链.

粒子BP算法应用MRF模型，采用最大积粒子BP算法^[9]，对每个节点s进行采样得到含有q个粒子的粒子集P_s={x_s⁽¹⁾…x_s^(q)}，任意节点处n次迭代后，节点s最大积置信度估计为b_sⁿ(x_s⁽ⁱ⁾)，与其对应的对数置信度为B_sⁿ(x_s⁽ⁱ⁾)=-lg(b_sⁿ·(x_s⁽ⁱ⁾)).节点s置信度定义为

从节点t到节点s的消息M_t→sⁿ(x_s)为

对于立体匹配算法，式(3)中ψ_s(u_s)为MRF模型中的数据项，用来计算局部势能，节点s表示为u_s=[d_s]=[a_sx_t+b_sy_t+c_s]，数据项ψ_s(u_s)定义为块内待匹配像素s的聚合匹配代价，表示为

ψ_st(u_s,u_t)为平滑项，代表{s,t}节点间的势能函数，求解在(x,y,disparity)视差空间中的两个局部平面差异.令节点s的平面法向量为n_s=orth([x_s,y_s,z_s]^T)，点u_s为这个平面上的一点，定义平滑项计算为

在迭代传播过程中，利用Metropolis-Hastings MCMC(MH MCMC)采样法^[9]对节点邻域进行重采样，在当前粒子集P_s及已传播的邻域候选粒子集C_s中选取能够使节点处数据项ψ_s最小的K个“最佳”粒子对当前粒子集进行更新.经n次消息和置信度的迭代计算，选择具有最小置信度的粒子集P_s作为立体匹配中的视差估计：

应用Patchmatch数据项计算匹配代价时通常采用尺寸较大的像素块，这样会丢失像素块的边缘特性.本文在Patchmatch相似度度量中，将像素色度空间和几何空间的邻近性引入权值计算，采用双边权重方法计算像素块之间的匹配代价，以此来达到还原输入视图细节，保持像素块的边缘特性，减少视差不连续区域的误匹配.像素匹配代价即数据项可以表述为

为了避免对建议分布的依赖，采用切片采样(slice sampling)粒子置信度传播来完成对节点邻域的重采样，选择合适的粒子集P_s来更新节点s的消息和置信度，使匹配能量最小.

将每个节点处置信度作为目标密度函数，节点s具有粒子集P_s={x_s⁽ⁱ⁾}_i=1,…,k，置信度传播的第n次迭代时，计算节点s的最大积置信度估计为b_sⁿ(x_s⁽ⁱ⁾)，对应对数置信度为B_sⁿ(x_s⁽ⁱ⁾).对节点s每个粒子x_s⁽ⁱ⁾的重采样从分布q(x_s⁽ⁱ⁾)中采样得到，构成MCMC采样链{x_s^(i)〈m〉}_m=1,…,M.其中节点s处的置信度q(x)可表示为

设辅助变量h∈R，则目标分布扩展为

通过对辅助变量h进行均匀采样形成切片，在范围h内每间隔A进行一次新的采样.

这里U_I为间隔I内的均匀分布，且间隔A={x;q(x)≥h^〈m〉}，图 1为切片采样示例.

图1(Figure 1)

图1 切片采样 Fig. 1 Slice sampling

令f_l(x)=m_t(l)→s(x_s⁽ⁱ⁾)，将q(x)分解，并通过L个辅助变量h₁,…,h_L得到采样，采样间隔为

即对于节点s的第i个粒子来说，在第m次MCMC迭代中以采样间隔A进行均匀采样，得到粒子x_s^(i)〈m〉(在下面论述中上标(i)〈m〉省略).经对数空间变换，节点s的对数置信度表示为

其中:F₀(x_s)=ψ_s(x_s);F_l(x_s)=M_t(l)→s(x_s);t^(l)表示节点s的第l个邻域t;|N_s|为最大邻域个数，在此采用节点四邻域关联则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为非负用户自定义常数.n_s和u_s表示为n_s=[nx_s,ny_s,nz_s],u_s=[x_s,y_s,z_s]，令′_l=(_l-B_t(x_t)+M_s→t (x_t))，经矩阵运算得到:

在采样间隔计算时，分别对x_s在(x,y,disparity)视差空间各个维度进行求解.因此，令Δx=|x_s－x_t|,Δy=|y_s－y_t|,Δz=|z_s－z_t|，因为在重采样过程中，要求重采样得到的粒子x_t应同样是处于像素块内，并且其视差应小于最大视差搜索范围，所以这里0≤Δx≤patchsize,0≤Δy≤patchsize,0≤Δz≤dip max，推导得到x_s在x维度采样间隔为

已知Δy和Δz取值范围，可求得X的取值范围，进而得到x维度采样间隔.同理可得x_s在y,z维度的采样间隔.

为了保证采样空间足够大，这里在计算X,Y,Z时，其中的Δx,Δy,Δz均取最小值.因此可以确定在不同维度下x_s的采样间隔.在进行采样时，则可在各自维度的采样间隔内依次在x维度、y维度、z维度抽取样本来构造马氏链.

将改进Patchmatch算法和切片采样粒子置信度算法相结合，视差估计算法分为初始化、传播、局部重采样、视差图后处理4个步骤.其中前3个步骤将视差估计转化为视图中节点的最佳粒子集及置信度的最优求解过程，算法具体描述：

输入：左右视图，P_s为视图中节点s的粒子集，K为粒子数，N为置信度传播迭代次数，M为切片MCMC重采样迭代次数.

输出：粒子及置信度最优解(P_s,B_s).

1) 对于视图中所有节点s进行初始化，每个节点s 进行K次随机采样，得到初始粒子集P_s={x_s⁽ⁱ⁾}_i=1,…,K；初始化消息及置信度m_t→s⁰(x_s)=0,B_s⁰(x_s⁽ⁱ⁾)=0；

2) for n=1 to N

3) for 所有节点s的i=1,…,K个粒子;

4) 初始化采样链x_s^(i)〈0〉←x_s⁽ⁱ⁾;

5) 按照传播策略函数对T(s)所有节点s的邻域N(s)进行均匀采样，得到候选粒子集 C_s=∪ {P_t|t∈N(s)∩T(s)};

6) C_s=P_s;

7) for m=1 to M

8) 对所有C_s中的粒子x_s(x_s∈C_s)在l=1,…,|N_s|邻域进行均匀采样h_l~U_[0,1](h)，采样得到_l=F_l(x_s^{(i)〈m－1〉})－lg(h_l);

9) 计算采样间隔A_(i)^〈m〉;

10) 采样得到x_s^(i)〈m〉~U_A(i)〈m〉(x);

11) 计算置信度.

12) 满足F_l(_s^(i)〈m〉) < _l，则接受_s^(i)〈m〉，即x_s^(i)〈m〉←_s^(i)〈m〉;

13) 更新粒子集P_s;

14) 对视图中所有节点s，返回粒子集P_s和置信度B_s，满足.

初始化及传播过程中在对节点s进行粒子随机采样时，可以对每个像素完成空间传播、视角传播及平面优化，对三维随机平面参数及视差进行随机选取，以得到近似平面.并通过视差图后处理步骤完成对遮挡区域的处理.

本文算法在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,σ)，则在对粒子集P_s进行粒子重采样时得到粒子集为

选取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 算法 Tsukuba Venus Teddy Cones 误匹配率 无遮挡区 整幅图片 视差不连续区 无遮挡区 整幅图片 视差不连续区 无遮挡区 整幅图片 视差不连续区 无遮挡区 整幅图片 视差不连续区 % Patchmatch 2.09 2.33 9.31 0.21 0.39 2.62 2.99 8.16 9.62 2.47 7.80 7.11 4.59 PM-Huber 3.49 4.09 9.12 0.22 0.43 2.50 3.38 5.56 10.7 2.15 6.69 6.40 4.56 PMBP 1.96 2.21 9.22 0.30 0.49 3.57 2.88 8.57 8.99 2.22 6.64 6.48 4.46 PMF 1.74 2.04 8.07 0.33 0.49 4.16 2.52 5.87 8.30 2.13 6.80 6.32 4.06 本文算法 1.85 2.31 8.03 0.35 0.37 2.43 3.02 5.52 9.53 2.10 6.75 6.19 4.42

表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个维度对粒子的切片采样范围进行计算，算法的执行时间有所增加，其执行效率有待提升.

1) 针对立体匹配时视差不连续区、倾斜平面及非前向平行平面误匹配较高问题，应用具有边缘特性的Patchmatch相似性函数，本文建立基于Patchmatch的非前向平行平面视差平面估计模型.

2) 采用切片采样马尔可夫链蒙特卡罗方法解决粒子置信度传播过程中粒子重采样更新问题，降低了算法复杂性，提高了匹配精度.

3) Middlebury图像数据集测试表明，该算法能够降低视差不连续区域的误匹配，有效提高了倾斜平面及非前向平行平面图像的匹配精度.