泛化置信度传播的双目立体视觉匹配方法 【技术领域】
本发明涉及图像处理、 计算机视觉、 计算方法、 数学、 数值方法领域， 尤其是计算机 视觉的双目立体视觉匹配方法。背景技术
目前， 立体视觉匹配问题的研究已经取得了很大的进展。特别是基于全局优化的 匹配算法， 已经成为了解决匹配问题的主要方法， 得到了广泛的运用。 其能得到如此关注的 原因是因为匹配问题可以很好地被建模为一个马尔可夫随机场 (MRF) 或条件随机场 (CRF) 的优化问题。这类问题在许多学科中都有设计， 由此产生的很多算法都可以运用到匹配问 题的解决当中来。
其中， 基于置信度传播 (BeliefPropagation) 的算法是一种目前受到广泛关注 的方法。它的主要思想是通过节点之间置信度的传播实现整个模型的逐步收敛。其最初 是在文献 (Pearl J..Probabilistic reasoningin intelligent systems ： networks of plausible inference[M]( 智能系统中的概率推演 ： 可信度推演网络 )， San Francisco ： Morgan KaufmannPublishers Inc.， 1988.) 中 提 出 的， 然 后 在 文 献 (Sun Jian， Zheng Nan-Ning， et al..Stereo matching using belief propagation[J]( 使用置信度传播的 匹配算法 ).IEEE transactions on Pattern Analysis and Machine Intelligence.2003， 25(7) ： 787-800.) 中第一次引入到求解匹配问题中。
置信度传播算法的一个最大的问题是只有在无环的图结构中才被证明可以收 敛， 而像匹配问题中使用到的网格状的马尔可夫随机场模 型具有大量的环状结构， 这样 就会造成结果的不稳定。对于这个问题主要有三种解决方法， 第一种是将迭代的次数设 定在一个比较小的次数上， 不用等到算法收敛就直接停止迭代。这也是最简单， 最常用的 一种解决方法 ； 第二种是将邻近两次迭代过程中置信度变化不大或没有变化的节点移除 出迭代过程， 使得当迭代过程的不断进行， 节点数量逐渐减少， 当没有节点参与迭代时， 迭 代自动终止 ； 第三种方法是在每次迭代中都先构造一棵生成树， 迭代在该生成树中进行， 如文献 (Wainwright M.J.， Jaakkola T.S.， et al..MAP estimation via agreementon trees ： message-passing and linear programming[J]( 使 用 树 的 统 一 性 的 最 大 后 验 概率估计 ： 信 息 传 递 和 线 性 规 划 ).IEEE Transactions onInformation Theory.2005， 51(11) ： 3697-3717.) 和 (Kolmogorov V.Convergent tree-reweighted message passing for energy minimization[J]( 使 用 收 敛 的 树 型 信 息 传 递 的 能 量 最 小 化 方 法 ).IEEE Transactions onPattern Analysis and Machine Intelligence.2006.10， 28(10) ： 1568-1583.) 中所示。 由于生成树中没有环状结构， 所以不会产生不收敛的结果。 而且文献 (Szeliski R.， Zabih R.， et al..A comparative study of energyminimization methods for markov random fields with smoothness-basedpriors[J]( 对基于光滑先验概率的 马尔可夫随机场进行能量最小化的算法的比较研究 ).IEEE Transactions on Pattern Analysis and MachineIntelligence.2008， 30(6) ： 1068-1080.) 中的实验表明， 该方法在匹配问题中可以得到比传统的置信度传播算法更稳定的匹配结果。
另外一种置信度传播算法的扩展是泛化置信度传播算法。该算法是传统置信度 传播算法的一种扩展。它最早是由 Yedidia 及其他研究者提出的 (Yedidia J.S.， Freeman W.T.， et al..Generalized belief propagation[J]( 泛 化 置 信 度 传 播 算 法 ).Neural Information ProcessingSystems.2000， 13 ： 689-695.)。泛化置信度传播算法的总体思 想是将图结构中的各节点进行聚类， 在各聚类间进行信息传播。它并没有制定具体的聚类 或者说分块策略， 只是一个框架性的算法。文献 (YedidiaJ.S.， Freeman W.T.， et al.. Constructing free-energy approximations andgeneralized belief propagation algorithms[J]( 构造自由能估计算法和泛化置信度传播算法 ).IEEE transactions on information theory， 2005， 51(7) ： 2282-2312.) 对泛化置信度传播算法的聚类方式进行了 讨论， 提出了用一种区域图的图结构来表示各种不同的聚类结果， 同时还在区域图上讨论 了泛化置信度传播算法的收敛性质。由于该算法计算量较大， 所以在匹配问题中还没有出 现应用实例。 发明内容
为了克服现有的针对泛化置信度传播算法的双目图像匹配方法的较高复杂度， 较 大计算量的不足， 本发明提出了基于最小和的缓存加速策略， 使得泛化置信度传播算法的 性能得到了很大加强， 从而成功将其应用于双目视觉中匹配问题的求解过程中， 提供了一 种有效降低复杂度、 减少计算量的泛化置信度传播的双目立体视觉匹配方法。
本发明解决其技术问题所采用的技术方案是 ：
一种泛化置信度传播的双目立体视觉匹配方法， 所述双目立体视觉匹配方法包括 以下步骤 ：
1) 采集双目的左右两幅图像， 将左图中的每一个像素点都作为一个变量， 然后保 持这些变量在图像坐标中的相对位置不变， 进行 4 邻域连接， 得到马尔可夫随机场的拓扑 结构， 然后按公式 (1) 和公式 (2) 分别计算马尔可夫随机场中的各变量的各状态代价值， 以及各连接边的代价值 ：
V(fp， fq) ＝ min(|fp-fq|， K)(2)其中， λ 表示代价权重， 它影响到点代价在整个能量函数中所占有的比重 ； fp 和 fq 分别表示变量 p 和 q 的状态序号 ； T 表示截断值 ； 颜色向量距离采用欧氏距离来表征， K表 L R 示截断值 ； Ic (p) 和 Ic (p) 分别表示左图和右图中 p 点处的 c 通道的颜色值 ；
2) 产生多尺度马尔可夫随机场， 第 k 层的大小是第 k+1 层大小的四分之一 ；
3) 设在多尺度马尔可夫随机场中共有 n 层， 按从 1 到 n 的顺序分别对 n 个马尔可 夫随机场进行求解 ； 首先， 使用泛化置信度传播算法对每一层的马尔可夫随机场分别求解， 原始的泛化置信度传播算法使用公式 (3) 和公式 (4) 来进行信息传递 ：
其中， φs ＝ φs(xs) ＝ D(xs)， ms → u ＝ ms → u(xu) 表示当 变量 u 选定状态 xu 时， 变量 s 向变量 u 传递的点信息， mst → uv ＝ mst → uv(xu， xv) 表示当变量 u 和变量 v 选定状态 xu 和 xu 时， 变量 s 和变量 t 之间的边向变量 u 和变量 v 之间的边传递 的边信息 ；
对公式 (3) 和公式 (4) 进行负对数操作， 并对其中的独立计算进行缓存， 得到两条 新的公式， 即公式 (5) 和公式 (6) ：
其中， 上标表示当前的迭代序号。均是缓存变量 ；将公式 (6) 的计算过程描述成一个在二维网格中搜索最小值的问题， 首先选取 ：
sini ＝ min{Qsu′ (xs， xu)} (7)
tini ＝ min{Qtv′ (xt， xv)} (8)
(7)、 (8) 作为搜索的起始点， 然后使点 P 在水平方向上移动， 搜索在该方向上的最 小值， 找到最小值后， 将点 P 移动到该位置 ； 接着将点 P 沿垂直方向进行搜索， 同样找到该方 向上的最小值， 然后移到该位置， 反复循环迭代以上步骤， 直到 P 点不再移动位置时， 其所 处位置的值即被认定为该二维网格中的最小值 ；
然后分点传递和边传递两部分分别进行信息传递， 在点传递过程中， 首先选择所 有非相邻的变量同步进行上下左右四个方向的传递， 再选择在前一步中没有进行传递的变 量进行同样方式的传递 ； 在边传递过程中， 分为水平边传递和垂直边传递两部分， 其中， 水 平边传递分过程中， 首先选择非相邻的水平边进行同步上下两个方向的传递， 再选择在前 一步中没有进行传递的水平边按同样方式进行传递 ； 然后 进行垂直边传递， 其传递方式和 水平边传递相同 ； 每一层次设定迭代次数， 进行迭代 ；
迭代完成后， 采用常值递减方法设置各层次的最少状态数量， 其定义如下所示 ：
其中， Γ(i) 表示第 i 层中最少的状态数量， d 表示在最稀疏的层次上保有的状态 数量， L 表示多尺度空间中的尺度数量， η 表示常值递减速率 ；当第 i 层中各变量在信息传递了一定次数收敛后， 先将最终计算得到的每个变量 的置信度向量中的各元素排序， 挑出最好的 Γ(i) 个状态， 作为下一层次的候选状态 ； 然 后， 将由点代价函数计算得到的每个变量的各状态的代价值进行排序， 选出最好的 Γ(i) 个状态， 作为下一层次的候选状态 ； 最后， 将两次计算得到的候选状态向量进行 “或” 操作， 得到最终传入下一层次的候选状态向量 ；
将第 i 层的计算结果传递继承到第 i+1 层 ；
4) 在最底层马尔可夫随机场求解完成后， 按下式计算每个变量的代价值 ：
然后取代价值最小的那个状态作为该变量的最终状态， 即为该变量所对应图像中 点的视差值。
进一步， 所述步骤 1) 中， 点代价计算是在 CIELAB 颜色空间中进行。
再进一步， 所述步骤 3) 中， 第 i 层的计算结果传递继承到第 i+ 1 层的继承过程分 为点信息继承和边信息继承两部分， 在点信息继承中， 第 i 层中的任意一个变量都有 4 条点 信息需要继承， 每一条点信息由第 i+1 层中相对应的 4 个变量继承 ； 在边信息继承中， 分为 水平边信息继承和垂直边信息继承两部分， 水平边信息继承中， 第 i 层中的任意一条水平 边都有 2 条边信息需要继承， 每一条边信息由第 i+1 层中相应的 2 条边继承， 垂直边的继承 过程和水平边继承过程相同。
本发明的技术构思为 ： 使用二路筛选策略在尺度空间中对状态空间进行缩减， 从 而加快整个算法的计算速度。本发明的有益效果主要表现在 ： 大大加速了泛化置信度传播 算法的计算速度。
附图说明 图 1 是多尺度马尔可夫随机场的建立过程的示意图。
图 2 和图 3 分别是在泛化置信度传播算法的一次迭代过程中两种信息的传递顺序 的示意图。
图 4、 图 5、 以及图 6 是多尺度马尔可夫随机场层次之间的信息继承过程的示意图。
图 7 是使用方向集算法进行复杂度降低的流程图。
图 8 是二路状态空间筛选策略的筛选过程的示意图。
具体实施方式
下面结合附图对本发明做进一步描述。
参照图 1 ～图 8， 一种泛化置信度传播的双目图像匹配方法， 所述计算机双目立体 视觉匹配方法包括以下步骤 ：
1) 使用左右两幅图像， 计算得到相对应的马尔可夫随机场。
2) 产生多尺度马尔可夫随机场， 第 k 层的大小是第 k+1 层大小的四分之一。
3) 设在多尺度马尔可夫随机场中共有 n 层， 按从 1 到 n 的顺序分别对 n 个马尔可 夫随机场进行求解。在计算过程中， 将第 i 层的计算结果传递到第 i+1 层。
4) 在最底层马尔可夫随机场求解完成后， 计算每个点最终状态值， 即每个像素点的视差值。
在步骤 1) 中， 将左图中的每一个像素点都作为一个变量， 然后保持这些变量在图 像坐标中的相对位置不变， 进行 4 邻域连接， 得到的就是马尔可夫随机场的拓扑结构。然后 按公式 (1) 和公式 (2) 分别计算马尔可夫随机场中的各变量的各状态代价值， 以及各连接 边的代价值。
V(fp， fq) ＝ min(|fp-fq|， K) (2)
其中， λ 表示代价权重， 它影响到点代价在整个能量函数中所占有的比重 ； fp 和 fq 分别表示变量 p 和 q 的状态序号 ； T 表示截断值 ； 颜色向量距离采用欧氏距离来表征， L R K 表示截断值 ； Ic (p) 和 Ic (p) 分别表示左图和右图中 p 点处的 c 通道的颜色值 ； 需要提 及的是， 为了使代价计算更精确， 这里的点代价计算是在 CIELAB(the L*a*b*standardof Commission Internationale de L′ Eclairage) 颜色空间中进行的。颜色向量距离采用 欧氏距离来表征。K 表示截断值。
在步骤 2) 中， 将马尔可夫随机场按图 1 产生多尺度结构。将每个 4 点变量集合中 的左上变量作为上一层次马尔可夫随机场中的变 量。
在步骤 3) 中， 首先， 使用泛化置信度传播算法对每一层的马尔可夫随机场分别求 解。原始的泛化置信度传播算法使用公式 (3) 和公式 (4) 来进行信息传递。
其中， φs ＝ φs(xs) ＝ D(xs)， ms → u ＝ ms → u(xu) 表示当 变量 u 选定状态 xu 时， 变量 s 向变量 u 传递的点信息， mst → uv ＝ mst → uv(xu， xv) 表示当变量 u 和变量 v 选定状态 xu 和 xu 时， 变量 s 和变量 t 之间的边向变量 u 和变量 v 之间的边传递 的边信息 ；
为了避免使用计算量较大的乘法和除法操作， 本发明对以上两条公式进行负对数 操作， 并对其中的独立计算进行缓存， 得到两条新的公式， 即公式 (5) 和公式 (6)。
8其中，CN 101901483 A
说明书6/8 页都是缓存变量。
通过分析公式 (6)， 发现其计算复杂度为 O(n4)， 是泛化置信度传 播算法的主要复 杂度。 本发明通过将一个在二维空间中搜索最小值的问题分解到多个一维空间中搜索最小 值的问题， 将其复杂度从 O(n4) 降低到 O(n3)。使用的方法为方向集方法， 也被称为 Powell 方法， 在文献 (Press W.H.， Teukolsky S.A.， et al..Numerical Recipes in C[M](C 语言 数值计算摘要 )， second edition.Cambridge ： Cambridge UniversityPress， 1992.) 中有 详细描述。该方法的思想是将一个 N 维的搜索问题转化到多个一维搜索问题， 实现问题规 模的线性化。将公式 (6) 的计算过程描述成一个在二维网格中搜索最小值的问题， 首先选 取：
sini ＝ min{Qsu′ (xs， xu)} (7)
tini ＝ min{Qtv′ (xt， xv)} (8)
(7)、 (8) 作为搜索的起始点， 然后使点 P 在水平方向上移动， 搜索在该方向上的最 小值， 找到最小值后， 将点 P 移动到该位置 ； 接着将点 P 沿垂直方向进行搜索， 同样找到该方 向上的最小值， 然后移到该位置， 反复循环迭代以上步骤， 直到 P 点不再移动位置时， 其所 处位置的值即被认定为该二维网格中的最小值 ； 图 7 显示了以上所述过程。
然后分点传递和边传递两部分分别进行信息传递， 在点传递过程中， 首先选择所 有非相邻的变量同步进行上下左右四个方向的传递， 再选择在前一步中没有进行传递的变 量进行同样方式的传递 ； 在边传递过程中， 分为水平边传递和垂直边传递两部分， 其中， 水 平边传递分过程中， 首先选择非相邻的水平边进行同步上下两个方向的传递， 再选择在前 一步中没有进行传递的水平边按同样方式进行传递 ； 然后进行垂直边传递， 其传递方式和 水平边传递相同， 即按图 2 和图 3 方 式进行顺序计算， 其中， 箭头方向表明消息的传递起始 点。每一层次设定一定的迭代次数， 进行迭代。
迭代完成后， 本发明使用尺度空间的状态空间缩减策略来对状态空间进行缩减。 在使用马尔可夫随机场建模的匹配问题中， 每个变量的状态数量是一个无法避免的参数， 而且随着两个摄像机之间距离的增加， 以及图像分辨率的增加， 该参数将会急剧增加， 从而 导致算法复杂度的指数级增长。 在分析算法的过程中可以发现， 在没有优化的情况下， 很多 的计算集中在一些对最终结果影响很小的状态值上。从而算法的复杂度越大， 像泛化置信 度传播算法， 这样的冗余计算量就越大。本发明提出的基于尺度空间的状态空间缩减策略 在尺度空间中逐步去除无关的状态， 使得算法的复杂度增长始终在较低的水平， 从而起到 加速作用。首先， 给定每一个层次应该保持的最少状态数量。这样的设定是为了在每一个
层次上都保证一定数量的状态， 从而使得信息能够得到充分传播， 另外一方面也是为了保 证状态空间不会过快缩减， 避免算法过快收敛， 保证能得到一个较好的收敛值。 具体的设定 值需要考虑两方面的因素， 一个是在每一层次中要保证一个足够大的状态空间， 使得信息 得到充分传播， 另一个是要保证计算量保持在一个较低的水平上。
由于这种设置不受每个变量的实际状态数影响， 而且在计算量和保证信息传递的 状态空间大小的平衡方面很难做出一个最佳决策， 更进一步说， 该设定可能和具体的应用 相关。本节介绍的方法将使用一个简单的统一设定， 即常值递减方法。其定义如下所示 ：
其中， Γ(i) 表示第 i 层中最少的状态数量， d 表示在最稀疏的层次上保有的状态 数量， L 表示多尺度空间中的尺度数量， η 表示常值递减速率。当给定最稀疏层次上保有的 状态数量后， 其他层次的状态保有数量都可以通过上式计算出来。为了适应不同的最大视 差范围， 即当最大视差范围大于 d 时， 根据点代价函数得到的最佳 d 个状态将作为最稀疏层 次上的 d 个候选状态。这里还要提及的是， 不光是点代价函数， 任意可以产生匹配准确性排 序的方法都能够采用， 本节直接采用点代价函数是出于简便考虑， 也可采用更复杂但更稳 定的方法。
在设置好了各层次的最少状态数量后， 本文采用一种二路筛选策略来选择一定数 量的状态进入下一层次的传播。首先， 当第 i 层中各变量在信息传递了一定次数收敛后， 先 将最终计算得到的每个变量的置信度向量中的各元素排序， 挑出最好的 Γ(i) 个状态， 作 为下一层次的候选状态。然后， 将由点代价函数计算得到的每个变量的各状态的代价值进 行排序， 选出最好的 Γ(i) 个状态， 作为下一层次的候选状态。最后， 将两次计算得到的候 选状态向量进行 “或” 操作， 得到最终传入下一层次的候选状态向量。以上描述的过程如图 8 所示。 其中， 浅色方格表示被选定的状态， 深色方格表示被删除的状态， 越大的点表示该状 态的惩罚值越大， 在图示的例子中， 该变量在该层次的候选状态数量为 5， 下一层中的最少 候选变量数， 即 Γ(i+1) 为 3。
在经过状态空间缩减后， 那些和保留的状态相关的信息将由下一层次继承 ；
第 i 层的计算结果传递继承到第 i+1 层的继承过程分为点信息继承和边信息继承 两部分， 在点信息继承中， 第 i 层中的任意一个变量 都有 4 条点信息需要继承， 每一条点信 息由第 i+1 层中相对应的 4 个变量继承 ； 在边信息继承中， 分为水平边信息继承和垂直边信 息继承两部分， 水平边信息继承中， 第 i 层中的任意一条水平边都有 2 条边信息需要继承， 每一条边信息由第 i+1 层中相应的 2 条边继承， 垂直边的继承过程和水平边继承过程相同， 继承过程如图 4， 图 5， 图 6 所示。
再在步骤 4) 中， 按下式计算每个变量的代价值
然后取代价值最小的那个状态作为该变量的最终状态， 即为该变量所对应图像中 点的视差值。
表 1 给出了使用本发明进行加速和不使用本发明时的泛化置信度传播算法的运 行效率的比较 ：
表1
本实验的实验环境为具有 1.6GHz 主频 CPU 和 1G 内存的个人电脑。图像对为由 Middlebury 提供 (D.Scharstein and R.Szeliski. “Ataxonomy and evaluation of dense two-frame stereo correspondencealgorithms( 一种对双目密集点匹配算法的分类和评 价 )” ， InternationalJournal of Computer Vision， 2002， 47(1)， pp.7-42.) 的 “Tsukuba” 实验 图像对， 大小为 384×288 像素。参数设置为 T ＝ 30.0， λ ＝ 0.87， K ＝ 10.0， d ＝ 16， η ＝ 3， 以及多尺度空间的层数为 5， 每个层次的迭代次数为 4 次。