基于半全局和全局算法的立体匹配研究

基于半全局和全局算法的立体匹配研究摘要：传统的基于像素点的匹配算法常常是算出初始匹配代价后直接采用贪心策略求取视差，虽然速度较快，但往往是局部最优的，以至精确度很低。针对这一问题，目前策略主要有：（1）半全局优化算法：扫描线算法和动态规划算法；（2）全局优化算法：置信度算法和图割算法。本文旨在通过详细讨论这四种算法原理本质，算法步骤与算法运行，从而深刻分析各自的优点与缺点，为进一步改进其不足，进而研究新的算法打下基础。关键词：半全局优化，全局优化，扫描线，动态规划，置信度，图割立体匹配介绍图像的立体匹配即给定同一场景的两幅图像,寻找同一场景点投影到图像中的像素之间的对应关系。根据考虑的是基于像素点的还是基于区域块，可以分为基于像素点的匹配与基于区域的匹配。立体匹配算法通常是通过构建能量函数试图获得图像的某些全局性质,即全局能量最小化，但通常很难获得能量函数的全局最小化，鉴于此,很多学者更倾向于寻找局部小的求解.然而在一般情形下,局部小不能带来任何的全局性，所以匹配效果较差，准确率较低，基于像素点的匹配就是一种局部小的解，所以若想提高精度，研究的多是一种半全局或全局优化策略的区域匹配算法。立体匹配的通常包括以下四步：1)图像预处理(Preprocessing)—由于拍摄照片的时候难免会有传感器的噪声（sensornoise）和光度的扭曲（photometricdistortions）而这都会对视差的计算带来严重影响，常用的解决方法有，高斯拉普拉斯滤波(LaplacianofGaussian(LoG)filtering)[1]直方图均衡化(HistogramEqualization/Matching)，中值滤波(Subtractionofmeanvaluescomputedintheneighboursofeachpixel)[2]双边滤波(Bilateralfiltering)[3]。2)匹配代价计算（Matchingcostcomputation）—对匹配代价的计算通常有四种方法AD（1-1）、SAD（1-2）、SD（1-3）与SSD（1-4），计算公式，从而能得到元素的不同视差匹配代价所组成的初始视差空间。(1-1)(1-2)(1-3)(1-4)3)视差的计算(Disparitycomputation)—真实的像素视差是指这两个像素点具有高的相似性，传统的WTA（WinnerTakesAll）算法就是每个像素点选取最小的代价来求取视差，是仅仅考虑一个像素的基于像素点匹配算法。如图(1)所示图1基于像素点的视差计算Figure1disparitycomputation图1基于像素点的视差计算Figure1disparitycomputationbasedonthepixeltopixelstereo图2视差优化Figure2disparityrefinement4）视差的优化(disparityrefinement.）—大多数立体匹配算法计算出来的视差是离散的，常常视差值都是整数，然而世界实际上是连续的，若想将立体匹配算法用在较高精度的场合，如机器人视觉，精密三维重建，这种离散的视差值不进行后续处理就无法达到令人满意的效果。针对这一问题，在获取初始视差后可以采用一些措施对视差进行细化，非整数视差，如匹配代价的曲线拟合如图2所示，或者直接采用亚像素精度法(sub-pixeldisparityEstimate)，即将原图像进行水平拉伸，再对行像素点进行模糊。本文将详细研究的就是立体匹配的视差计算阶段，就是在初步计算出来的视差空间中进行半全局或全局优化，从而求出更好的视差值。与基于像素点的匹配方法不同，基于半全局或全局的算法通常是将匹配问题转换为一个能量方程，然后通过求解该能量方程的最小值来求取视差值。能量方程通常具有以下的形式(1-5)其中Cdata(dx)是数据项用来约束像素点在偏移前后的变化尽量小,V(dx,dx-1)是光滑项，约束像素点在偏移前后与周围像素点的关系变化尽量小。不同的优化算法的数据项和光滑项的定义不同，本文采取的数据项计算方法为AD。而光滑项则根据不同的算法用不同的模型方法。二．匹配算法与实现 2.1扫面线优化方法SO(ScanlineOptimization)方法属于一种半全局能量最小化优化算法，比基于像素点匹配的算法有更好的适应性，通过最小化自身的能量方程来求得像素对应的视差，但是由于其只考虑一行的最优，并未考虑所有的像素点，所以本文称之为半全局能量最小化算法。2.1.1算法原理：SO的能量方程的光滑项通常为：(2-1)(2-2)其中opt_smooth表示图片中像素点x的光滑值，opt_grad_penalt用来对同纹理区域(灰度梯度值小于给定的阈值opt_grad_penalty)的区域进行惩罚的系数，从而不会在同纹理区域由于匹配代价相同而导致视差计算错误。SO是逐行考虑的算法，本质就是计算从每一行第一个像素点到当前像素点x在视差为dx时候的匹配代价，用sumcost(x,dx)表示，(2-3)其中(2-4)而Cost(dx)表示从前一个像素点视差为dx-1到当前像素点在视差为dx时候的匹配代价。(2-5)当保存到本行最后一列元素x=W时候，此时sumcost保存了x=W处的所有视差对应的最小匹配代价，选取最小的sumcost，进行回溯，此时即得到本行最优视差路径。然后同理计算下一行元素。从而能得出所有像素点的视差值，而SO的最小化能量方程转化为：(2-6)2.1.2算法步骤：初始化第一列元素(x=1):令sumcost(x;1,2,3,…,dmax)=Cdata(dx=1,2,3,…,dmax)；从下一列(x=x+1)开始,取不同的视差d=1，2,3…dmax,计算像素x-1从不同的视差dx-1=1,2,…,dmax到当前具有视差d的像素x所对应的最小代价bestcost，从而得出该元素的sumcost(x,dx),并且保存x各个视差所对应的前一个像素的视差值best_d。判断若x<=W则转步骤（2）否则继续；此时的sumcost(w;1,2,…dmax)保存了从第一列到当前列的最小能量，则此时从min(sumcost(x;1,2,3,…,dmax))所对应的d开始回溯,根据之前保存的各个像素对应的best_d来得到本行的视差集；考虑下一行元素，同样采取这样的方法。2.1.3算法运行详解：为简单起见设一行仅含有5个像素点，视差的范围为d=1,2,3,4。边的权值表示从像素点x-1的各个视差到像素点x所对应的视差d的最小能量代价bestcost(x,d)。如图3每个点记录了从x=1到当前像素点在视差d时候的最佳匹配代价sumcost(x,dx)，例如在x=2,d=2处一共有四条路径(路径1，路径2，路径3，路径4)从中选取bestcost(2)所对应的路径(路径4)。最后在x=5处求出最小的匹配代价即min(sumcost(5,d=1,2,3,4))，此时进行回溯则找到图中的红色路径即为该行的最佳视差即d=(4,2,2,4,1)。图3SO算法运行图Figure3SOAlgorithmdiagram图3SO算法运行图Figure3SOAlgorithmdiagram比起基于像素点的立体匹配，它能保证在一行是最优的，速度快，而且与下面的动态规划相比，允许相邻匹配点拥有更大的视差变化范围。不足是，首先它仅仅对水平扫描线进行优化寻径，造成扫描线间的视差不连续；其次，它使用固定的视差搜索范围，导致一定的冗余计算。目前的改进方法是添加更多的约束条件，来减少冗余的视差扫描，同时点使扫描线间的视差也有一定的约束作用。2.2动态规划优化方法在实际应用中通常立体视觉匹配算法要具有较高的准确性和实时性。为了提高立体匹配的精度，一般采用全局算法。它使用全局约束来解决匹配问题，通常利是用全局优化算法来最小化能量函数，但全局算法的时间消耗异常的大，而动态规划(DynamicProgramming)是属于半全局算法中的一种优化方法，它避免了NP问题，而且具有计算效率高、匹配效果较好的特点，因而这种算法在实际场合中是最常被使用。动态规划立体匹配算法分为两个阶段：1）计算左右两幅图像的原始匹配代价，建立视差空间图像；2）在视差空间图像上运用动态规划算法进行优化，基于动态规划的立体匹配算法的实质就是要在视差空间图像上，考虑各种视差变化情况下，求出一条最优路径，使能量函数E取值最小，从而求取视差。为进一步限制路径搜索的范围，降低算法的复杂度，文献[4]引入了唯一性约束和顺序性约束，从而限定了从前一阶段到当前阶段的转换只有7种可能的形式，如图4所示。分别用m、l、r表示前一阶段(pre_state)的匹配点、左遮挡点（即只能在参考图中被看到）、右遮挡点（即只能在匹配图中看到），用M、L、R表示当前阶段(current_state)的匹配点、左遮挡点、右遮挡点。我们称此动态规划为三状态动态规划。图47种形式转换图Figure47formsconversionchart2.2.1算法原理对于DP算法的核心部分就是从每一行第一个像素点到当前像素点x在视差为d，状态为current_state时候的最小匹配代价，用sumcost(d,x,current_state)表示。(2-7）其中sumcost(pre_d,pre_x,pre_state)表示像素点x的前一个位置(此时的视差为pre_d,状态为pre_state,位置为pre_x)最小匹配代价。Cost表示从前一个位置到当前位置的转移代价，转移代价参考图5，则可得计算方法为(2-8)CoccL与CoccR为左遮挡和右遮挡的惩罚代价，光滑约束项V(dx,dx-1)的定义与SO算法相同。图5代价转换Figure5Costoftransition直到本行最后一列元素x=W时候，此时sumcost(d,W,current_state)共保存了dmax个最小匹配代价即sumcost(d=1,2,…,dmax,W,current_state)，与SO算法相同，选取最小的匹配代价集min(sumcost(d=1,2,3,…,dmax,W,current_state))，进行回溯，此时即得到本行最优视差路径,然后同理计算下一行元素。与上面的SO算法不同的是，DP在计算匹配代价的时候考虑了遮挡问题,从而添加了遮挡惩罚代价，能较好的处理遮挡问题。对常规能量方程式使用动态规划求的最佳路径可转化为：(2-9)2.2.2算法步骤：(1)初始化的第一个点x=1的sumcost，设置其为匹配状态；(2)转到x的下一个位置,讨论从前一个像素点位置pre_x到当前像位置x的变换方式t与从第一个点到当前点的匹配代价sumcost。方法为：根据前一个像素点位置pre_x，状态pre_state与视差pre_d讨论像素点pre_x到像素点x的变换t，计算出转换代价Cost(t)，从而得到第一个像素点到当前像素点的最小匹配代价sumcost(d,x,current_state);并且保存x各个视差所对应的pre_x处的最佳转换方式best_t。若x<W则继续转(2)否则继续执行；(3)回溯寻找最优的路径,先求出最后一列元素即x=W时的最小匹配代价（注意默认其为无遮挡匹配像素）sumcost(d,W,current_state=1)根据当前的d,计算出从前一个像素点到当前像素点的转换方式best_t;从而根据转换度，得出前一个位置的视差pre_d注意，若当前的状态为匹配则视差值设置为d，否则标志其为空洞，即：(4)填充空洞，根据回溯时候的视差dispMap，若当前像素点为遮挡点occLable，则使其视差值为最近一处非遮挡点的视差值。2.2.3算法运行详解：如图所示，假设有六个像素点，视差d=1,2,3,在最末处x=6处寻找最小的sumcost(d=1,2,3,6,current=1)，然后进行回溯，根据保存的best_t寻找x=5处的视差，以此类推，最后得到一条从x=1到x=6的路径，即得到的视差为d={1,1,1,2,3,2,1,1},而x=3,x=4处为左遮挡区域，x=5处为右遮挡区域其余为匹配点。图6DP算法运行图Figure6DPAlgorithmdiagramDP不仅能保证在一行是最优的，速度快，而且能考虑到遮挡区域。但是由于每一条扫描线被单独地处理，缺少了扫描线之间的约束限制，所以该算法和SO算法一样存在条纹瑕疵改进方法常常是利用扫描线间的相关信息来约束动态规划过程，使得处理每一条扫描线时都能充分利用以前扫描线的匹配信息。例如Stephen[5]等引入了控制点修正技术，改善了视差图质量，大大减少了条纹瑕疵.；GongandYang[6]提出了一个基于可靠性的二通道的DP算法，它沿着水平和垂直路径利用动态规划，充分利用了像素间信息；Coxetal[7]提出了利用扫描线间的一致性约束，实现了二维的内聚性。虽然这些改进对DP都起到了一定的优化作用，但在优化的同时常常会引进新的缺点，所以都不是很理想。2.3置信度优化方法2.3.1算法原理：很多计算机视觉问题都具有不适定性，处理这类问题要引入合适的先验约束，把它转化为适定的问题。其中马尔科夫随机场是描述计算机视觉问题的常用方法。2003年Sun[8]等人将置信传播算法(BeliefPropagation)用在立体匹配中，取得了异常好的效果，此后基于置信度传播的立体匹配算法成为人们研究的热点。BP算法任然是将视差表示成标号形式,为匹配问题构建能量函数,在适当的约束下,通过邻接及转移关系建立马尔可夫网络,巧妙的使能量与网络中的信任度关联。最后通过求取全局能量函数的最小值来求取视差。由于基于置信传播算法的立体匹配不是单单一行，所以其能量模型应变为(2-10)：（2-10）N表示相邻像素对的集合，P表示所有的像素点集合对式子两边进行指数变换可得（2-11）令则能量方程可变为（2-12）从而可以将立体匹配全局能量最小化问题转化为最大积马尔科夫随机场问题，从而利用置信传播算法求解立体匹配问题的全局近似最优解。BP算法的主要思想是消息的迭代传输，假设xs在xt的左边，用msti+1(xs,xt)表示第i次迭代中节点xs传递给xt的左消息，记作msti+1(xleftmessage）,传递信息计算参考图7,bs(xs)记为结点xs的置信度，要不断的更新，更新方式如图8。其计算方法分别为：（2-13）即得到xt点处的左信息矩阵xleftmessage（2-14）其中mprawmessage，mrightmessage，mleftmessage，mupmessage，mdownmessage分别表示xs处，节点的信息，右信息，左信息，上信息和下信息最后xs点的置信度为（2-15）图7x信息传递给y

图7x信息传递给y

图8图8信任度传播2.3.2算法步骤初始化四个存储信息传递矩阵，初始值均设置为1即Leftmessage,当信息从左向右传播时候，保存左信息；rightmessage,当信息从右向左传播时候，保存右信息；upmessage,当信息从上向下传播时候，保存上信息；downmessage当信息从下向上传播时候，保存下信息；prawmessage,源信息矩阵，即为初始匹配代价；根据公式（2-13）更新每个点的信息矩阵，即右信息矩阵mrightmessage，左信息矩阵mleftmessage，上信息矩阵mupmessage和下信息矩阵mdownmessage。根据公式（2-14）计算节点的最佳置信度xsopt同时记录下各元素的最大的置信度（2-15）所对应的视差d，此时便求得该处的视差值。最后判断若达到迭代次数则停止否则转(3)2.3.3算法运行详解：如图显示了在MarkovNetwork网络中消息传递的方式，在最大积算法中，消息传递从x6到x7便为x7的leftmessage，即同理求出x7处的rightmessage,upmessage,downmessage,则可以得到x7处的信任度b7为此时根据最小的b7所对应的视差来得出x7处的视差值best_d。图9BP算法运行图BP算法的优点是目标函数的凹凸性不受限制,计算复杂度高,获得的全局最优路径与真实路径有一定差距。目前针对BP的不足改进方法有：Felzenszwalb[9]等人采用的分层等技术对SBP加速，使得全局能量能够较快收敛；Yang[10]等人利用快速收敛策略也使算法能够较快的执行；BeliefProg[11]算法利用极值位置关系进行加速处理，能使运行效率大约提高了30％~60%。2.4图割优化方法GC(GraphicCut)[12]图割算法是另外一种全局优化算法，即用标号表示视差,建立能量函数,把匹配问题转化为能量函数最小化问题，再通过构造网络,使能量与网络的割的容量相联系，最后利用图的网络流理论给出能量函数的最小化,从而获得图像匹配的视差数据，根据图的构造方式不同，图割算法共有两个即Swap算法与算法2.4.1Swap算法原理Swap的思想是通过不断的选取两个不同的视差，通过其构造网格来与能量方程想关联，P为所有元素的集合，而Pαβ是视差为的元素集合，网格构成如图11所示，图的构成元素有：顶点集边集为顶点，Pαβ中的元素分别与相连接构成t-link边，而Pαβ中内部元素相连接构成n-link边,权值赋值的方式分别为(t-link)(n-link)则割集即其中则将最小割集与能量函数连接了起来，求出最小割集则为最小的能量，同时也能得出视差集。2.4.2swap算法步骤用简单的方法比如WTA方式获得初始视差集labels也可以默认都为1，设置标志success=0从labels中任意选出两种视差，将原来的视差为的像素点和视差值为的像素点形成一个标号集f，并用此标号集按照swap构造方法构造一个网格，计算出此标号集下f的能量方程值E，将能量方程对应到这个网格中，求出此网格的最小割，此时与源source相连接的像素点视差为，与源sink相连接的像素点的视差为，此时得到一个新的标号集f’，计算在此情况下的能量方程新的值E’若E’<E则表示在此情况下的标号集比原来的标号集能量值更低则得到一个新的标号集f=f’设置能量变小标志success=1判断若success=1则转步骤（2），否则表示能量已经不能再减小优化完毕，则停止循环2.4.3swap算法运行详解：假设已经由简单的匹配算法如WTA获得初始视差，每个像素点获得相应的标号与视差值标号集L={f1,f2,f3,…,f12}视差集P={1,2,3,4}图10带初始视差的像素点随机选取交换的视差=3与=2，从而构成新的标号集L’={f4,f5,f7,f`10,f12}而视差集即为P’={3,2}图11Swap算法运行图对应于此例如图所示,求出此时的最小割集C，即为求的最小的能量E’若此时的E’<E则开始修改视差，即f4,f5,f12与相连则这三个像素点的视差值均为=3，f7,f10与相连则该两个像素点的视差均为=2，其中f12的视差得到了修正2.4.4算法原理：的思想是通过对视差集中的每一个视差，通过其构造网格来与能量方程想关联，如图所示，图的构成有顶点集边集L中的元素分别与相连接构成t-link边，而L中元素相连接构成n-link边权值赋值的方式分别为t-linkn-link则割集即则将最小割集与能量函数连接了起来，求出最小割集则为最小的能量，同时得出视差集。2.4.5算法步骤（1）用简单的方法比如WTA方式获得初始视差集labels也可以默认都为1，设置标志success=0（3）对于P中每一个视差，将原来的视差为的像素点和其他视差值的像素点形成一个标号集f，并用此标号集构造一个网格，计算出此时的能量方程值E将能量方程对应到这个网格中，与swap不同的是以前只是考虑两种视差之间的变换，现在是视差为的像素点与其他所有像素点之间的变换，而且添加了辅助结点，求出此网格的最小割，此时与源source相连接的像素点视差为，与源sink相连接的像素点的视差为，此时得到一个新的标号集f’，计算在此情况下的能量方程新的值E’若E’<E则表示在此情况下的标号集比原来的标号集能量值更低则得到一个新的标号集f=f’设置能量变小标志success=1（4）判断若success=1则转步骤（2），否则表示能量已经不能再减小优化完毕，则停止循环2.4.6算法运行详解：仍然假设已经由简单的匹配算法如WTA获得初始视差，每个像素点获得相应的标号与视差值标号集L={f1,f2,f3,…,f12}视差集P={1,2,3,4}如上图所示对于视差集中每一个视差，图12算法运行图对应于此例如图所示,求出此时的最小割集C，即为求的最小的能量E’若此时的E’<E则开始修改视差，f1,f2,f3,f12与相连则这三个像素点的视差值均为f4,f11与相连则该两个像素点的视差保持不变f12的视差得到了修正由于GC是全局优化，计算量相当的大，所以在时间和空间的消耗上都比半全局优化大的多。针对效率的改进方法很多，例如国内的一种受限α-扩展(α-expansion)操作[13],该操作根据灰度连通性和特征点匹配的结果对每次网络构造的顶点进行控制,减少网络中顶点和边的数目,可有效提高计算效率。三实验分析3.1参数定义与设置：本文用文献[14]的方法来评估算法的性能，即分别对各自算法总的错误像素点百分比bad_pixels_all，无纹理区域T错误像素点百分比bad_pixels_textureless，不连续区域D错误像素点百分比bad_pixels_discont做了统计。(3-1)(3-2)(3-3)dC(x,y)为算法计算所得的视差，dT(x,y)为真实的视差；（eval_bad_thresh）为容许的视差错误，在本实验取值为1.0；无纹理区域T即为在指定大小的方形区域(宽度为eval_textureless_width)内视差平均值小于给定的阈值(eval_textureless_thresh)；本实验指定eval_textureless_width=3,eval_textureless_thresh=4不连续区域D即为在像素点的视差与其相邻(evaldiscontwidth）.的视差值大于给定的阈值(evaldispgap)；本实验指定eval_discont_width=9，eval_disp_gap=2;此外由于很多算法没有考虑边界问题，所以本文在比较算法的时候除去了边界像素点，边界定义的宽度为eval_ignore_border=18;本文图片来自/网站的数据集，选用了tsukuba图片与均为经过校正的图片，去除了干扰噪声，从而能将注意力放在算法的本身性能。右图左图右图左图3.2首先对tsukuba图片进行评估对于SO算法关键的参数就是对于数据项Cdata计算时候的截断值Truncation，与计算光滑项时候用到的opt_smooth，对这两个参数分别做实验DP算法根据DP算法的原理可以知道opt_smooth的值越高，视差图就越平滑，同时也会丢失一些细节信息；opt_smooth的值越低，视差图就更精细，但是会引入更多的异常值。CoccR/L的作用是控制视差图中允许的遮挡数量，过高会导致一些正确的遮挡丢失，相反，如果过低，那么视差图中会出现一些多余的遮挡。BP算法BP算法是通过迭代来进行优化的，由于BP算法每次迭代都非常耗时，所以希望迭代次数越少越好，实验证明tsukuba仅仅迭代3~4次已经能达到收敛BP算法的截断函数与前面的不同，采用了GC算法参考文献[1]T.Kanadeetal.,“DevelopmentofaVideo-RateMachine”,Pro.ofInt.RoboticandSystemsConf.(IROS’95),PittsburghPA.,Aug.,1995.[2]Faugerasetal.Real-Timecorrelation-basedstereo:Algorithm,ImplementationandApplications,INRIATR2013,1993.[3]Ansaretal.Enhancedreal-timestereousingbilateralfiltering.CVPR2004.[4]AaronF.Bobick,StephenS.Intille.Largeocclusionstereo[J].InternationalJournalofComputerVision,1999,33(3):181-200.[5]StephenS.Intille,AaronF.Bobick.Disparity-spaceimagesandlargeocclusionstereo[C]//Stockholm,Sweden:ThirdEuropeanConfComputerVision,1994[6]MinglunGong,Yee-HongYang.Nearreal-timereliablestereomatchingusingprogrammablegraphicshardware[J].ComputerVisionandPatternRecognition,2005,1(1):924-931[7]Ingemar