CN103020964B - 一种基于自适应分水岭的图割的双目立体匹配方法 - Google Patents

Abstract

本发明提出一种基于自适应分水岭的图割的双目立体匹配方法，本发明提出了新的自适应局部阈值方法，并将其应用于分水岭结合Prim方法的区域融合中。使用自适应分水岭对图像进行处理，使图中像素以一定关系将图像分割成不同的区域并分配标号来建立能量方程，并提出新的大领域无参数的平滑约束模型。最后通过优化的α-扩展法，利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则不搜索。大邻域能够使本发明更灵活地获得更好的目标边界，在不连续的边界和高折叠纹理区域都得到很好的效果。在不影响合成视图质量的同时解决了静态图像和动态视频序列的矛盾性，深度估计的结果不会呈现不连续的跳跃现象。

Description

一种基于自适应分水岭的图割的双目立体匹配方法

技术领域：

本发明涉及一种双目数字图像的处理方法，特别涉及一种基于自适应分水岭的图割的双目立体匹配方法。

背景技术：

随着多媒体技术的不断发展，图像和视频技术也由二维向三维发展，交互性将成为未来多媒体技术的一个主要特征。交互式三维视频系统的关键技术之一就是虚拟视点绘制合成。虚拟视点绘制合成是所有立体显示系统终端不可缺少的模块，在远程视频会议、自由视点立体电视等高端多媒体领域中也具有至关重要的作用。为了使用户可以在场景中漫游，实现“连续的环视”，在多视点视频采集的过程中，摄像机的数量应尽可能多，但由于放置无限个摄像机以实现视点无缝切换的不现实性，为了显示任意视点的视图，必须在客户端进行虚拟视点的合成，通过对已有视点的分析，合成用户所要观察的视点。因此，虚拟视点合成技术是多媒体领域一项非常重要的新兴技术。首先需要有精确的立体匹配方法获得视差图或深度图才能完成虚拟视点合成技术。

与灰度图像相比，深度图像具有物体三维特征信息，即深度信息。由于深度图像不受光源照射方向及物体表面的发射特性的影响，而且不存在影响，可以得到三维物体更可靠的几何信息，所以更准确地表现物体目标表面的三维深度信息。深度图在视频编码中的作用很大，可以有效的提高多视图图像传输和多视图视频传输的编码效率。正因为如此，深度图像分析越来越受到计算机视觉、图像分析等研究领域的重视，在工业领域特别是在机器人视觉、自动导航、工业零件的自动检测和自动装配等领域，得到了越来越广泛的应用。

对于深度图的获取方面，图割作为一种基于图论的组合优化技术，在用来最小化计算机视觉中的能量函数问题上被众多研究者所使用，目前已有许多新技术应用于该领域。RichardSzeliski等人把当前常用的几种能量函数最小化方法在解的质量和运行时间等方面进行了比较，发现与传统的8种方法（模拟退火方法、M-估计法等）相比，图割方法不仅总体精度高，而且在不连续区域和低纹理区域的精度也比其它方法都高。（参见DScharstein,RSzeliski.Ataxonomyandevaluationofdensetwo-framestereocorrespondencealgorithms[J].InternationalJournalofComputerVision,2002,47(1):7-42.）。不仅如此，即使有些方法（如模拟退化方法）的精度和图割方法的精度接近，但图割方法在优化过程中收敛更快（参见RSzelisk,IRZabih.Anexperimentalcomparisonofstereoalgorithms[A].ProceedingsoftheInternationalWorkshoponVisionAlgorithms:TheoryandPractice[C].Springer-VerlagLondon,UK.LectureNotesinComputerScience,2000,1883:1-19）。利用图割方法最小化能量方程函数可以将图像映射为网络图，图像的特征就可以用图论的方法进行处理，在优化能量函数时，运用图割的方法在二值标号问题中可以得到能量函数的全局最小；在多标号问题中可以得到带有很强特征的局部最小，图割方法还可以保证能量函数的解收敛到全局最小，且实际效率较高，得到的数值解有很强的鲁棒性。

在Middlebury网站中几乎所有的这些好的方法在立体匹配时都使用图像分割。基于分割的立体匹配方法能够很好实现平坦区域的重建。这些方法的出现有一个潜在的问题就是在应用在动态视频序列中。图像分割在视频帧间是相互矛盾的，深度估计的结果常常呈现不连续的跳跃现象。例如，HaiTao,HarpreetS,SawhneyI,RakeshKumar提出了将3D场景通过图像分割将颜色或分为不同平面，他们假设同一颜色的图像区域与三维表面一致。这个想法鼓舞了很多现有的有关立体匹配的研究。这种模型使用一个递增方程来进行估计。这个方程能够优化与空间颜色相关的一致性和平滑项的能量方程（参见H.Tao,H.S.Sawhney,R.Kumar.Dynamicdepthrecoveryfrommultiplesynchronizedvideostreams.CVPR,2001.）。O.J.Woodford,P.H.S.Torr,I.D.Reid,A.W.Fitzgibbon提出基于“QPBO”方法的扩展法来有效的优化能量方程，但是这会带来三倍于二阶约束项的计算量。然而这种方法对于平坦区域能够有很好的效果，对于处理纹理表面比如有折叠在不同方向上有不同纹理区域效果并不好（参见O.J.Woodford,P.H.S.Torr,I.D.Reid,A.W.Fitzgibbon.Globalstereoreconstructionundersecondordersmoothnesspriors.CVPR,2008.）。又如Tsin的方法核心是使用3D点作为无参数平滑约束。但是，没有提出全局方法来最小化它们的能量方程。相反，他们使用每个像素“winner-take-all”的估计策略，其对初始深度估计敏感（参见Y.Tsin.KernelCorrelationasanaffinitymeasureinpoint-sampledvisionproblems.PhDthesis,RoboticsInstitute,CarnegieMellonUniversity,September2003.）。王年，范益政，鲍文霞等提出一种基于图割的匹配方法。相对于以往的基于图割的方法该方法不需要摄像机的位置信息以及运动信息，并把标号从一维推广到二维，从而使本方法适用于更一般情形下的匹配问题，此外，为了更利于网络的构造和方法的实现，此文献中通过用像素的梯度向量的距离来代替平滑项（参见王年,范益政,鲍文霞等.基于图割的图像匹配方法.电子学报,2006,34(2):232-235.）。张令涛，曲道奎，徐方提出了一种基于图割的改进立体匹配方法，方法通过区域匹配方法得到每个像素的初始视差值，然后只保留完整网格图的部分可能的视差值，去除其余大部分的节点和边缘，建立简化的网格图，该方法大大缩减了网格图的容量，缩短匹配所用的时间，并且能够选用更大的视差范围（参见张令涛,曲道奎,徐方.一种基于图割的改进立体匹配方法.机器人,2010,32(1):104-108.）。朱程辉，任冉冉提出一种快速立体匹配方法，把图像分割成颜色单一的不同区域；计算初始视差图，利用可靠点求取各分割区域的平面模板参数，对模板参数相同的相邻区域进行融合；构造全局能量方程，采用图割方法求取全局能量最小的视差最优分配。该方法对低纹理区域和遮挡区域有较好的匹配结果（参见朱程辉,任冉冉.一种基于图割理论的快速立体匹配方法.微型机与应用,2010,10:35-38.)。本发明对Middlebury平台中，YangQ,YangR,DavisJ,NisterD的方法DoubleBP（参见YangQ,YangR,DavisJ,NisterD.Stereomatchingwithcolor-weightedcorrelation,hierarchicalbeliefpropagationandocclusionhandling[J].PatternAnalysisandMachineIntelligence,2009,31(3):492-504.）；KlausA,SormannM,KarnerK的方法AdaptBP（参见KlausA,SormannM,KarnerK.Segmentbasedstereomatchingusingbeliefpropagationandaselfadaptingdissimilaritymeasure[A].Proceedingsofthe18thInternationalConferenceonPatternRecognition(ICPR2006)[C],2006:15-18.）；伍春洪，付国亮的方法K-均值分割（参见伍春洪,付国亮.一种基于图像分割及邻域限制与放松的立体匹配方[J].计算机学报.2011,34(4):755-760.）；SangHwaLee,SiddharthSharma的SAD方法（SangHwaLee,SiddharthSharma.Real-timedisparityestimationalgorithmforstereocamerasystems[J].IEEETransactionsonConsumerElectronics,2011,57(3):1018-1026.）的相关数据来进行比较，详见表1。

虚拟视点绘制中立体匹配技术还有诸多环节和关键方法上存在较大的改进空间。基于分割的立体匹配方法能够很好地实现平坦区域的重建。这些方法应用在动态视频序列中就不能解决视频帧间相互矛盾的问题，深度估计的结果常常呈现不连续的跳跃现象。

发明内容：

本发明要解决的技术问题是：为克服现有技术的不足，本发明提供一种基于自适应分水岭的图割的双目立体匹配方法，实现一种在稀疏图下的像素匹配，使不连续的边界保留得很好而不需要图像分割作为预处理项。本发明避免了预处理中使用分割，并且在动态视频序列中能够恢复实时稳定的深度值，在每一帧单独处理时也能得到很好的效果。本发明的约束项模型是大领域无参数的。无参数模型将图像特征用深度值表示时不需要使用明确的固定阶的约束项，这样计算量就大大降低。大邻域能够使本发明更灵活地获得更好的目标边界，在不连续的边界和高折叠纹理区域都得到很好的效果。在不影响合成视图质量的同时解决了静态图像和动态视频序列的矛盾性，深度估计的结果不会呈现不连续的跳跃现象。实验结果表明：通过Middlebury平台对本发明方法定量评估得出在所有区域误匹配、非遮挡区域以及深度不连续区域的误匹配率都控制在8.5%以内，在Middlebury平台135组数据中排名第19位。

本发明具有稳定可靠、精度高、抗干扰能力强的优点。另外，通过对表1对比，可以得出：本发明属于全局优化方法，与局部优化方法K-均值分割和SAD方法相比较，从表1不同类方法错误百分比的数据比较中可以看出，本发明的各项错误率均明显低于局部优化方法。而同类方法中，本发明的错误百分比略高于平台上最优方法AdaptBP，但低于平台上最优方法DoubleBP。

本发明解决其技术问题所采用的技术方案包括：一种基于自适应分水岭的图割的双目立体匹配方法，其特征在于包括以下步骤：

（1）输入拍摄自同一场景，同一时刻的两幅图像，这两幅图像在拍摄视角上和平行度上都可以存在差异；

（2）判断两幅输入图像与摄像机个数是否匹配，若不同，提示错误并跳出；若相同，读取相关数据，加载图像并执行步骤（3）；

（3）一种自适应分水岭方法，提出了新的自适应局部阈值方法，并将其应用于分水岭结合Prim算法的区域融合中。具体是这个方法包括两个主要步骤：首先，使用分水岭分割方法将图像分割成大量的区域；第二步是一个重复的过程，在此区域被融合且达到局部阈值停止融合。融合过程的顺序参照Prim最小生成树方法，在图表中找出最小值生长树的方法，在融合过程中我们追踪每个区域的变化并将变化的特征保存下来；

（4）根据步骤（3）中求出的融合区域作为标号，建立能量方程，使用大领域无参数深度平滑模型来建立图割的能量方程的立体匹配方法；

（5）为一种能量函数最小化的方法即立体匹配的方法，具体是对于步骤（4）中的能量方程进行最小化过程。使用优化的α-扩展法，寻找匹配点时不需要对整幅图进行搜索，而是利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则不搜索。

上述所述步骤（3）中的基于自适应分水岭建立标号，采用以下步骤实现：

（ⅰ）对加载的图像应用Canny边缘检测算子来得到梯度度量的图像；

（ⅱ）经过分水岭分割，并使其投影到彩色图像，输出即是将I_waterleft，I_waterright分别分割成n个不重叠的过分割的区域图像；

（ⅲ）使用Prim方法来生成最小生成树，令G=(V，E)为RAG结构，表示对图像I_waterleft的初始分割，其中e(i,j)的权值为函数的值。产生MST的过程即融合区域的过程；

（ⅳ）一般的融合方法都是只设置单个阈值，达到这个值就停止融合，但是这样容易造成不必要的误差，本发明提出一种自动计算局部阈值（即自适应局部阈值），这样经过Prim算法后得到的每个区域由于阈值不同所以大小是不同的，并且可以根据所需精度进行调整。由于融合不同区域时，区域的同一性质就会产生比较大的变化，本发明利用这个变化来确定局部阈值，简单来说就是融合过程中，一旦融合的区域不同就停止融合；

（ⅴ）通过（ⅰ）-（ⅳ）步骤获得一个T树。继续使用Prim方法和自适应局部阈值方法得到最小生成树来完成估计最初密集图。从原始图排除树边界，剩余图仍然很密集。为了更好的估计，在剩余图的第二个树仍然使用Prim方法。此外，反复寻找T树，并且合并所有这些树去合成稀疏图以估计原始密集图，这样一个稀疏图G_s至多有T(L-1)个边界。

上述所述步骤（4）中，基于步骤（3）建立能量方程，方法如下：

（ⅰ）图的一致项其中d_p=D₁(p)是图Il中像素p的视差，q=p+D_l是图I_r中p的相对应的像素，d_q=D_r(q)是I_r中q的视差。

（ⅱ）选用的平滑项如下：

$E_{\mathrm{smooth}}\left(D\right)=\underset{p\∈I}{\mathrm{\Σ}}{E}_{\mathrm{smooth}}(d_p;{\left\{d_q\right\}}_{q\∈N_p})$

其中E_smooth模型连接像素p的视差d_p和像素q在p的邻域N_p的视差d_q。

（ⅲ）使用5D矢量f＝[x,c]作为特征矢量，包括位置x和颜色c。定义像素p以及其的邻域N_p，基于邻域N_p使用非参数归正预测视差d_q。建立连接在N_p视差d和特征f的分类：

$P(d,f|N_p)=\frac{1}{\left|N_p\right|}\underset{q\∈N_p}{\mathrm{\Σ}}{g}_d\left(\frac{d-d_q}{{\mathrm{\σ}}_d}\right)g_x\left(\frac{x-x_q}{{\mathrm{\σ}}_x}\right)g_c\left(\frac{c-c_q}{{\mathrm{\σ}}_c}\right)$

其中g_d，g_x和g_c分别是视差d，像素位置x和像素颜色c的核心方程。

上述所述步骤（5）中在能量方程最小化的方法即优化后的α-扩展方法实现为：

（ⅰ）初始化能量函数的值为0，根据视差范围设置缓冲区的个数，并初始化迭次数为0；

（ⅱ）产生标号的随机排列，依次选择标号进行α-扩展操作，利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则匹配点不在该区域内，不进行搜索；

（ⅲ）标号集中的标号循环一次后输出一个能量值，迭代次数增加1，重复（ⅱ）的操作；

（ⅳ）直到迭代次数达到最大值或者缓冲区个数减为0时，方法结束；

（ⅴ）根据步骤（ⅰ）-（ⅳ）完成能量函数最小化的过程之后，根据得到的视差分布最终得到深度图I_LD，I_RD。

本发明与现有技术相比所具有的优点：

（1）提出新的自适应局部阈值方法，通过识别融合过程中每个区域是从何处开始变得不同类的设置阈值，达到阈值停止融合，并将其应用于分水岭结合Prim方法的区域融合中。基于对现有图割方法的研究，进一步设计了基于自适应分水岭方法并且使用非参数深度平滑模型来建立图割的能量方程的立体匹配方法。所得的最小生成树组就为像素组层，该像素组层作为能量方程的标号；

（2）使用大领域的密集型估计建立一个非参数深度平滑模型将图像特征与深度值连接。这个邻域明确叙述了全局立体匹配，并且使用图割在这样两个或多个像素组层网络来匹配像素。立体匹配公式使用不同的平滑规则来恢复表面，比如高纹理细节和不连续区域。经实验证明，本发明的方法对高纹理遮挡区域的细节和不连续区域的边界都能得到很好的效果；

（3）能量函数最小化使用优化的α-扩展法，寻找匹配点时不需要对整个图进行搜索，而是利用最小生成树的区域（即像素组层）中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则不搜索，这样大大减小了搜索的数据量。

附图说明：

图1是本发明基于自适应分水岭的图割的双目立体匹配方法流程图；

图2是本发明中Prim方法流程图；

图3是本发明中能量方程最小化流程图；

图4是原始视图组1中的左视图；

图5是原始视图组1中的右视图；

图6是本发明结果中组1的左视图的深度图；

图7是本发明结果中组1的右视图的深度图；

图8是原始视图组2中的左视图；

图9是原始视图组2中的右视图；

图10是本发明结果中组2的左视图的深度图；

图11是本发明结果中组2的右视图的深度图。

具体实施方式：

下面将结合附图对本发明方法作进一步的详细说明，使用自适应分水岭方法将输入图像进行初始分割，每个区域作为一个点云，使用Prim方法将点云融合，通过识别融合过程中每个区域是从何处开始变得不同类的设置阈值，达到阈值停止融合，形成一棵最小生成树T。剩余点云继续使用Prim方法和自适应局部阈值方法得到其余的生成最小生成树，直到所有点云全部融合形成所需的像素组层。依据像素组层建立能量方程，使用优化后的α-扩展法，寻找匹配点，最终得到所需的深度图。

本实施例中所有“视图”、“图像”均指数字位图，横坐标为从左至右，纵坐标为从上至下，均从0开始计数，使用5D矢量f＝[x，c]作为特征矢量，包括位置x和颜色c。图1所示为本发明的基于自适应分水岭和图割的立体匹配方法流程图；图2所示为本发明中Prim方法的流程；图3所示为能量函数最小化的流程图，具体步骤如下：

（1）输入拍摄自同一场景，同一时刻的两幅图像图4和图5（或图8和图9），这两幅图像在拍摄视角上和平行度上都可以存在差异。对一组输入图像图4和图5（或图8和图9）满足双目原理，分别相当于用“左右眼”对场景进行观察和记录。而本发明所要完成的工作即是匹配并输出左右两幅图像相应的深度图；

（2）判断图4和图5（或图8和图9）两幅输入图像与摄像机个数是否匹配，若不同，提示错误并跳出；若相同，读取相关数据，加载图像并执行步骤（3）；

（3）一种自适应分水岭方法，提出了新的自适应局部阈值方法，并将其应用于分水岭结合Prim算法的区域融合中。具体是这个方法包括两个主要步骤：首先，使用分水岭分割方法将图像分割成大量的区域；第二步是一个重复的过程，在此区域被融合且达到局部阈值停止融合。融合过程的顺序参照Prim最小生成树方法，在图表中找出最小值生长树的方法，在融合过程中我们追踪每个区域的变化并将变化的特征保存下来；

（4）根据步骤（3）中求出的融合区域作为标号，建立能量方程，其中数据项搜索是对两幅图单独进行，平滑项使用大领域无参数模型，使用步骤（3）所得的区域进行处理；

（5）为一种能量函数最小化的方法即立体匹配的方法，具体是对于步骤（4）中的能量方程进行最小化过程。使用优化的α-扩展法，寻找匹配点时不需要对整个图进行搜索，而是利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则不搜索。

其中，上述所述步骤（3）中的基于自适应分水岭建立标号方法，采用以下步骤实现：

（ⅰ）对加载的图像I_l，I_r进行灰度处理并应用Canny边缘检测算子来得到梯度度量的图像I_lc，I_rc；

（ⅱ）经过分水岭分割，并使其投影到彩色图像，输出即是将I_waterleft，I_waterright分别分割成n个不重叠的过分割的彩色图像，为了得到所需要的图像，使用区域融合。将这n这区域表示为其中i∈[1,n]，r_i∈[1,R_i]。n为区域数量，R_i为融合过程中融合的数目。在这里假设(i∈[1,n])是原始区域的序列，是分水岭方法在重复融合之前的输出。为了确定融合顺序，设任意两个不同区域和间相异点函数为函数f包含色彩和边缘两个分量。对于每个区域它的色调分量平均值表示为而边缘分量使用梯度大小来表示。区域和的平均梯度为由两个区域共有的像素梯度计算出来的。图像I_waterleft中的像素值即为梯度值，M_ij表示区域和边界像素的集合。为：

${\mathrm{\μ}}_{\mathrm{waterleft}}(N_i^{r_i},N_j^{r_j})=\frac{{\mathrm{\Σ}}_{(x,y)\∈M_{\mathrm{ij}}}{I}_{\mathrm{waterleft}}(x,y)}{\left|M_{\mathrm{ij}}\right|}$

其中|M_ij|表示M_ij中的像素数量。

$f\left(N_i^{r_i}{N}_j^{r_j}\right)={\mathrm{\ω}}_1\×d({\mathrm{\μ}}_h\left(N_i^{r_i}\right),{\mathrm{\μ}}_h\left(N_j^{r_j}\right))+{\mathrm{\ω}}_2\×{\mathrm{\μ}}_{\mathrm{waterleft}}(N_i^{r_i},N_j^{r_j})$

其中是和的平均值之差，即：

$d({\mathrm{\μ}}_h\left(N_i^{r_i}\right),{\mathrm{\μ}}_h\left(N_j^{r_j}\right)=\min \{|{\mathrm{\μ}}_h\left(N_i^{r_i}\right)-{\mathrm{\μ}}_h\left(N_j^{r_j}\right)|,(360-|{\mathrm{\μ}}_h\left(N_i^{r_i}\right)-{\mathrm{\μ}}_h\left(N_j^{r_j}\right)\left|\right)\}$

ω₁和ω₂是预定义的恒定系数；

（ⅲ）使用Prim方法来生成最小生成树，令G=(V，E)为RAG结构，表示对图像I_water0的初始分割，其中e(i,j)的权值为函数的值。G的生长树为一个贯穿它所有节点的非周期连接子图。其中当|V|=n时，G的每个生成树有n-1个边。每个边都有1个权值，G的最小生成树就是总权值之和最小的生成树。在本发明中设边的集合为A。

一般的融合方法都是只设置单个阈值，达到这个值就停止融合，但是这样容易造成不必要的误差，本发明提出一种自动计算局部阈值（即自适应局部阈值），这样经过Prim算法后得到的每个区域由于阈值不同所以大小是不同的，并且可以根据所需精度进行调整。由于融合不同区域时，区域的同一性质就会产生比较大的变化，本发明利用这个变化来确定局部阈值，简单来说就是融合过程中，一旦融合的区域不同就停止融合。

令表示区域的亮度分量平均值，且令V(x,y)表示坐标(x,y)处的亮度分量值。任何一个区域i＝1,2,3…n，在它的第r_i次融合时产生的变化定义为：

$\mathrm{\σ}\left(N_i^{r_i}\right)=\frac{1}{\left|N_i^{r_i}\right|}\sqrt{\underset{(x,y)\∈N_i^{r_i}}{\mathrm{\Σ}}{(V(x,y)-{\mathrm{\μ}}_V\left(N_i^{r_i}\right))}^2}$

其中表示的大小。

区域在第r_i次融合时的变化为 $\mathrm{\Δ\σ}\left(N_i^{r_i}\right)|\mathrm{\σ}\left(N_i^{r_i}\right)-\mathrm{\σ}(N_i^{r_i}-1)|.$ 令J表示r_i=1,…,R_i局部最大值集合，局部阈值由下式决定：

其中 $\mathrm{\δ}=\frac{1}{K}\underset{r_i\∈J_i}{\mathrm{\Σ}}\mathrm{\Δ\σ}\left(N_i^{r_i}\right)$ 是的平均值。

产生MST的过程即融合区域的过程，记述如下，设TE是G上最小生成树中边的集合：

1、从U＝{u₀}(u₀∈V)，开始；

2、在所有u∈U，V∈V-U的边(u,v)∈E中找到一条权值的边(u₀,v₀)归入TE中，同时v₀并入U中；

3、重复执行步骤2，一旦符合局部阈值决定公式就停止融合形成一棵最小生成树T_i；

4、重复步骤2、3，直到U＝V为止。

（ⅳ）通过（ⅰ）-（ⅲ）步骤获得一个T树。继续使用Prim方法和自适应局部阈值方法得到最大生成树来完成估计最初密集图。从原始图排除树边界，剩余图仍然很密集。为了更好的估计，在剩余图的第二个树仍然使用Prim方法。此外，反复寻找T树，并且合并所有这些树去合成稀疏图以估计原始密集图，这样一个稀疏图G_s至多有T(L-1)个边界。

上述所述步骤（4）中，根据基于步骤（3）建立能量方程方法具体如下：

（ⅰ）图的一致项其中d_p=D₁(p)是图I_l中像素p的视差，q=p+D_l是图I_r中p的相对应的像素，d_q=D_r(q)是I_r中q的视差。e_data定义为：

$e_{\mathrm{data}}(d_p,d_q)=\left\{\begin{array}{ccc}0& \mathrm{if}& d_p<d_q\\ {\mathrm{\&rho;}}_{\mathrm{data}}\left({\left|\right|c_p-c_q\left|\right|}^2\right)& \mathrm{if}& d_p{=d}_q\\ 0& \mathrm{if}& d_p>d_q\end{array}\right.$

其中c_p=I_l(p)，c_q=I_r(q)，ρ_data是图一致性的一个常数。令ρ_data(x)=min(0,|x|-τ_data)，对于多目立体视频，图的一致性是建立在一组被选择的图像对中。

（ⅱ）选用的平滑项如下：

$E_{\mathrm{smooth}}\left(D\right)=\underset{p=I}{\mathrm{\&Sigma;}}{E}_{\mathrm{smooth}}(d_p;{\left\{d_q\right\}}_{q\&Element;N_p})$

其中E_smooth模型连接像素p的视差d_p和像素q在p的邻域N_p的视差d_q。

（ⅲ）使用5D矢量f＝[x,c]作为特征矢量，包括位置x和颜色c。定义像素p以及其的邻域N_p，基于邻域N_p使用非参数归正预测视差d_q。建立连接在N_p视差d和特征f的分类：

$P(d,f|N_p)=\frac{1}{\left|N_p\right|}\underset{q\&Element;N_p}{\mathrm{\&Sigma;}}{g}_d\left(\frac{d-d_q}{{\mathrm{\&sigma;}}_d}\right)g_x\left(\frac{x-x_q}{{\mathrm{\&sigma;}}_x}\right)g_c\left(\frac{c-c_q}{{\mathrm{\&sigma;}}_c}\right)$

其中g_d，g_x和g_c分别是视差d，像素位置x和像素颜色c的核心方程。对像素位置x和颜色c使用高斯核心方程：σ_d，σ_x和σ_c分别是相关的频带宽度。为预测d_p，根据给定的f_p计算出f_p所有可能值：

$P\left(d\right|f_p,N_p)=\underset{q\&Element;N_p}{\mathrm{\&Sigma;}}{w}_{p,q}{g}_d\left(\frac{d-d_q}{{\mathrm{\&sigma;}}_d}\right)$

其中

$w_{p,q}=\frac{g_x\left(\frac{x-x_q}{{\mathrm{\&sigma;}}_x}\right)g_c\left(\frac{c-c_q}{{\mathrm{\&sigma;}}_c}\right)}{{\mathrm{\&Sigma;}}_{q\&Element;N_p}{g}_x\left(\frac{x-x_q}{{\mathrm{\&sigma;}}_x}\right)g_c\left(\frac{c-c_q}{{\mathrm{\&sigma;}}_d}\right)}$

这样平滑项E_smooth就成为：

E_smooth(d_q;N_p)=-λlog(P(d_p|f_p,N_p))

其中λ是归一化系数。不需要使用整幅图像去评估P(d_p|f_p,N_p)，只要在核心方程的帮助中来评估分配加权w_p,q。

因为log(·)是一个凸函数。因此P(d_p|f_p)的表达式，有以下E_smooth的不等式：

$E_{\mathrm{smooth}}(d_p;N_p)<-\mathrm{\&lambda;}\underset{q\&Element;N_p}{\mathrm{\&Sigma;}}{w}_{p,q}\log \left(g_d\left(\frac{d-d_q}{{\mathrm{\&sigma;}}_d}\right)\right)$

让 ${\mathrm{\&rho;}}_{\mathrm{sm}}(d-d_q)\underset{\&OverBar;}{\underset{\&OverBar;}{\mathrm{def}}}-\log \left(g_d\left(\frac{d-d_g}{{\mathrm{\&sigma;}}_d}\right)\right).$ 得到E_smooth上限值只包含两个视差变量。使用E_smooth上限值来代替总能量方程中的立体模型，对E有上限值E^u：

E^u(D_l,D_r)=E_smooth(D_l,D_r)+E^u _smooth(D_l)+E^u _smooth(D_r)

原则上，图割能够最优化E^u；因为包括大量的对偶项，所以执行起来计算量过大。

为使最优化方程更有效，注意到很多w_p，q的权值很小，比如那些包括两个位置距离很远的像素值或两个颜色不相似的像素。选择忽略这些项，并且去寻找一组优势项来接近

通过步骤（3）得到的稀疏图来规整每个像素，在这个图中每个像素是一个结点和在p和q之间组成了一个权值边界w_p,q。这个图没有直接获得而是通过用一个权值w_p,q+w_q,p来连接边界p-q和q-p。这个图是密集的，我们希望通过一个总边界权值最大化的稀疏图G_s来估计密集图。我们最终目标是在这个稀疏图中边界相关项来估计用于定义估计出得上限值

$E_{\mathrm{smooth}}^u\&ap;E_{\mathrm{smooth}}^{\mathrm{au}}\left(D\right)=\mathrm{\&lambda;}\underset{(p,q)\&Element;G^s}{\mathrm{\&Sigma;}}{w}_{p,q}{\mathrm{\&rho;}}_{\mathrm{sm}}(d_p-d_q)$

其中上标“au”代表“估计的上限值”。用来取代方程E^u(D_l,D_r)=E_smooth(D_l,D_r)+E^u _smooth(D_l)+E^u _smooth(D_r)中的对于能量总方程中原始立体模型中的Φ有估计上限值E^au：

$E^{\mathrm{au}}(D_l,D_r)=E_{\mathrm{smooth}}^{\mathrm{au}}(D_l,D_r)+E_{\mathrm{smooth}}^{\mathrm{au}}\left(D_l\right)+E_{\mathrm{smooth}}^{\mathrm{au}}\left(D_r\right)$

上述所述步骤（5）中在能量方程最小化的方法即优化后的α-扩展方法实现为：

（ⅰ）初始化能量函数的值为0，根据视差范围设置缓冲区的个数，并初始化迭次数为0；

（ⅱ）产生标号的随机排列，依次选择标号进行α-扩展操作，利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点，否则匹配点不在该区域内，不进行搜索；

（ⅲ）标号集中的标号循环一次后输出一个能量值，迭代次数增加1，重复（ⅱ）的操作；

（ⅳ）直到迭代次数达到最大值或者缓冲区个数减为0时，方法结束；

（ⅴ）根据步骤（ⅰ）-（ⅳ）完成能量函数最小化的过程之后，根据得到的视差分布最终得到深度图I_LD，I_RD即图6和图7（或图10和图11）。

本发明方法的实验平台为：Core^TM2DuoProcessorCPU，主频1.86GHz，3.25GDDR2内存。本发明所提出的立体匹配系统的编程环境为MicrosoftVisualStudio2008。使用Middlebury平台^[19]提供的四组标准测试图像。实验结果表明：通过Middlebury平台对算法定量评估得出在所有区域误匹配、非遮挡区域以及深度不连续区域的误匹配率都控制在8.5%以内，在Middlebury平台135组排行数据中排名第19位。

如表1所示为基于Middlebury平台本发明方法与其它方法误匹配率(％)结果的比较，测试图像为Middlebury平台上所提供的标准测试图像Teddy和Cones。

表1

注：all表示所有的区域误匹配率，Nonocc表示非遮挡区域误匹配率，disc.表示深度不连续区域误匹配率，匹配结果相差2像素以上认为匹配错误。

通过表1中对比，可以得出：本文方法属于全局优化方法，与局部优化方法均值分割和SAD方法相比较，从表1中不同类方法错误百分比的数据比较中可以看出，本文方法的各项错误率均明显低于局部优化方法。而在同类方法中，本文方法的错误百分比略高于Middlebury平台上的最优方法AdaptBP（排名第4），但部分错误百分比低于Middlebury平台上的最优方法DoubleBP（排名第6）。

Claims (4)

1.一种基于自适应分水岭的图割的双目立体匹配方法，其特征在于包括以下步骤：

(1)拍摄同一场景和同一时刻的两幅图像，输入这两幅图像；

(2)判断两幅输入图像与摄像机个数是否匹配，若不同，提示错误并跳出；若相同，读取相关数据，加载图像并执行步骤(3)；

(3)一种自适应分水岭方法，提出了新的自适应局部阈值方法，并将其应用于分水岭结合Prim方法的区域融合中；具体这个方法包括两个主要步骤：首先，使用分水岭分割方法将图像分割成大量的区域；第二步按照Prim方法，即在图表中利用Prim方法对使用分水岭分割方法将图像分割成的大量区域进行融合，在融合过程中追踪每个区域的变化并将变化的特征保存下来；重复以上区域融合过程，直到达到局部阈值则停止融合；

(4)根据步骤(3)中求出的融合区域作为标号，建立能量方程，使用大邻域无参数深度平滑模型来建立图割的能量方程的立体匹配方法；

(5)一种能量函数最小化的方法即立体匹配的方法，具体是对于步骤(4)中的能量方程进行最小化；使用优化的α-扩展法，寻找匹配点时不需要对整幅图进行搜索，而是利用Prim方法求得的最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点。

2.根据权利要求1所述的一种基于自适应分水岭的图割的双目立体匹配方法，其特征在于：所述步骤(3)采用以下步骤实现：

(ⅰ)对加载的图像应用Canny边缘检测算子来得到梯度度量的图像；

(ⅱ)经过分水岭分割，并使其投影到彩色图像，输出即是将I_waterleft，I_waterrigh分别分割成n个不重叠的过分割的区域图像；

(ⅲ)经典最小生成树方法有Kruskal方法和Prim方法；通过比较，两者都属于贪心方法，而运行时间上，Prim方法优于Kruskal方法；本发明使用Prim方法来生成最小生成树，令G＝(V,E)为相邻区域邻接图结构，G＝(V,E)表示对图像I_waterleft的初始分割；为了确定融合顺序，设任意两个不同区域和间相异点函数为其中e(i,j)权值为函数的值；产生MST的过程即融合区域的过程；

(ⅳ)本发明提出一种自动计算局部阈值，这样经过Prim方法后得到的每个区域由于阈值不同所以大小是不同的；由于融合不同区域时，区域的同一性质就会产生比较大的变化，本发明利用这个变化来确定局部阈值，简单来说就是融合过程中，一旦融合的区域的同一性质不同就停止融合；

(ⅴ)通过(ⅰ)-(ⅳ)步骤获得一个T树；继续使用Prim方法和自适应局部阈值方法得到最小生成树来完成估计原始密集图；从原始图排除树边界，剩余图仍然很密集；为了更好的估计，在剩余图的第二个树仍然使用Prim方法；此外，反复寻找T树，并且合并所有这些树去合成稀疏图以估计原始密集图。

3.根据权利要求1所述的一种基于自适应分水岭的图割的双目立体匹配方法，其特征在于：所述步骤(4)中基于步骤(3)建立能量方程，方法如下：

(ⅰ)图的一致项d_p＝D_l(p)是图I_l中像素p的视差，q＝p+D_l是图I_r中p的相对应的像素，d_q＝D_r(q)是I_r中q的视差；e_data(d_p,d_q)定义为：

$e_{data}(d_p,d_q)=\left\{\begin{array}{ccc}0& if& d_p<d_q\\ {\mathrm{\&rho;}}_{data}\left(\right||c_p-c_q|{|}^2)& if& d_p=d_q\\ 0& if& d_p>d_q\end{array}\right.$

其中c_p＝I_l(p)，c_q＝I_r(q)，ρ_data是图一致性的一个常数；令ρ_data(x)＝min(0,|x|-τ_data)，对于多目立体视频，图的一致性是建立在一组被选择的图像对中；

(ⅱ)选用的平滑项如下：

$E_{smooth}\left(D\right)=\underset{p\&Element;I}{\&Sigma;}{E}_{smooth}(d_p;{\left\{d_q\right\}}_{q\&Element;N_p})$

其中E_smooth模型连接像素p的视差d_p和像素q在p的邻域N_p的视差d_q；

(ⅲ)使用特征矢量f＝[x,c]作为特征矢量，包括位置x和颜色c；定义像素p以及其邻域N_p，基于邻域N_p使用非参数归正预测视差d_q；建立连接在Np视差d和特征矢量f＝[x,c]的分类：

$P(d,f|N_p)=\frac{1}{\left|N_p\right|}\underset{q\&Element;N_p}{\&Sigma;}{g}_d\left(\frac{d-d_q}{{\mathrm{\&sigma;}}_d}\right)g_x\left(\frac{x-x_q}{{\mathrm{\&sigma;}}_x}\right)g_c\left(\frac{c-c_q}{{\mathrm{\&sigma;}}_c}\right)$

其中，g_d，g_x和g_c分别是视差d，像素位置x和像素颜色c的核心方程；σ_d、σ_x和σ_c分别是相关的频带宽度。

4.根据权利要求1所述的一种基于自适应分水岭的图割的双目立体匹配方法，其特征在于：所述步骤(5)中在能量函数最小化的方法即优化后的α-扩展方法实现为：

(ⅰ)初始化能量函数的值为0，根据视差范围设置缓冲区的个数，并初始化迭代次数为0；

(ⅱ)产生标号的随机排列，依次选择标号进行α-扩展操作，利用最小生成树的区域中像素范围来搜索，在范围内搜索区域像素，寻找匹配点；

(ⅲ)标号集中的标号循环一次后输出一个能量值，迭代次数增加1，重复(ⅱ)的操作；

(ⅳ)直到迭代次数达到最大值或者缓冲区个数减为0时，方法结束；

(ⅴ)根据步骤(ⅰ)-(ⅳ)完成能量函数最小化的过程之后，根据得到的视差分布最终得到深度图I_LD，I_RD。