本文是立体视觉部分的第四篇，立体匹配。主要介绍了立体匹配的算法思路，详细介绍了SGM算法，并在最后给出了代码实例。关于立体匹配的其他内容，请移步本系列另外几篇博客。

立体匹配

获得两幅行对齐的图像后，就可以设计高效的立体匹配算法了。已知左图上的点(x0, y0)，右图与之匹配的点一定在(mindisparity, Maxdisparity)之间。其中, NumDisparities一般是预设的固定值。

已知左图某点(x0,y0)，为其在右图上寻找匹配点，最朴素的想法是，计算该点与右图上(x0-MaxDisparity, y0)到(x0-minDisparity, y0)之间所有像素的差异(代价，cost)。取差异最小的一个作为他的匹配点，这就是所谓的WTA(Winner take all)算法。

对左图每一个像素计算代价，然后找出最优，得到全图的视差图。

过程中，所有像素的dmin到dmax的代价组成了DSI空间(disparity space image). DSI中的元素C(x,y,d)代表了Ir(x, y)与It(x+d, y)的相关性(代价)。

诸多问题

如上文所描述，寻找匹配点似乎并不是一件很困难的事。然而实际情况并非如此。童话里的故事都是骗人的。

传感器感光的差异和噪声。

左右摄反射光线的差异。

透视收缩。

纹理缺失

重复纹理

遮挡和不连续性

以上这些都会影响匹配的结果。那么如何设计一个好的匹配算法？

论文A taxonomy and evaluation of dense two-frame stereo correspondence algorithms, D. Scharstein and R. Szeliski 给出了立体匹配算法的几个步骤，

计算匹配的代价(Matching cost)
代价聚合(Cost aggregation)
视差计算
视差优化

另外，从种类上来说立体匹配算法可以分为局部算法，全局算法以及半全局算法。这三种算法的主要区别在于第二步上。

代价计算

无论什么样的算法，第一步都需要定义像素和像素之间的相关性，也就是匹配代价(matching cost)。常见的代价有，

绝对值差异（Absolute differences）

$C(x,y,d) = \begin{vmatrix} I_{R}(x,y)-I_{T}(x+d,y) \end{vmatrix}$

平方差 (Squared differences)

$C(x,y,d)=(I_{R}(x,y)-I_{T}(x+d,y))^_{2}$

BT(Birchfied and Tomasi). 出自一篇经典论文，S. Birchfield and C. Tomasi. A pixel dissimilarity measure that is insensitive to image sampling. Opencv中SGBM采用了此方法。

Census。目前基于Census的变种使用非常广泛，效果也不错。原始论文 Non-parametric local transforms for computing visual correspondence

首先对左右图的每个像素做Census变换，

定义P为任一像素，D为其邻域像素集合， $\xi$ (P, P+[i,j])的意思是，如果P+[i,j]小于P，函数返回1，否则返回0。最后将所有的返回值组合成二进制。

C(x, y, d)是两个census的汉明距离。

例如，

左图点P和其邻域的灰度值为，

127 127 129

126 128 129 => census: 11010101 （中间像素与8邻域作比较，邻域像素小于中间像素，置1)

127 131 111

右图点P+d(其实应该是-d)及其邻域灰度值为，

121 180 110

130 131 111 => census: 10111111

127 129 110

两个census的汉明距离： 4

代价聚合

如果直接采用上述的代价做WTA(winner take all)寻找d，那么相当于只考虑了单像素的匹配程度。假若像素处于弱纹理或者重复纹理区域，这个代价就不能准确表达像素的差异。代价聚合建立了相邻像素之间的联系，建立了平滑性约束，即假设相邻的像素应有相近的视差。这个思路广泛用于局部算法和半全局算法。

局部匹配(Local Matching)

在为每一个像素计算代价后，局部匹配算法会采用基于窗口的代价聚合算法作平滑性处理。以SSD算法(Sum of Squared differences) 算法为例，他的算法流程如下，

通过计算灰度的平方差，计算像素所有d的代价（DSI）。
在一个固定窗口内，为该像素的每个d做代价聚合S。比如说，某一个特定的d，该点的聚合代价等于这个窗口内的所有像素视差为d时的代价和。 S为(x,y,d)所在的窗口像素的集合。 $S(x, y, d)= \sum_{(x',y')\in S}C(x',y',d)= \sum_{(x',y')\in S}(I_{R}(x',y')-I_{T}(x'+d, y'))^{2}$
WTA。对代价聚合的结果S，选聚合代价最小的d作为视差。

类似的算法还有SAD(Sum of Absolute difference)。

这种固定窗口的算法很明显有很多弊端。比如说，窗口内采用同一视差的代价和，其实是在假设窗口内的视差一致，但实际上是没法保证这一点的。

因此，有很多进一步优化的算法，如Box-filtering优化等，就不全面展开了。

全局匹配(Global)

典型的全局匹配算法是不做代价聚合的。他将视差问题转换成MRF问题，通过优化一个全局能量函数估计视差。能量函数包含数据项和平滑项，通过最小化这个函数，获得最优的视差。

$E(d)=\sum_{p\in R}C(p,d)+\sum_{q,p\in R}P(d_{q}-d_{p})$

p代表某一像素，q是其相邻像素，R是像素集合。 $\sum_{p\in R}C(p,d)$ 是数据项，定义了像素代价的匹配程度。 $\sum_{q,p\in R}P(d_{q}-d_{p})$ 是平滑项，约束了邻域之间关系。P是两个相邻像素p, q的惩罚函数，相邻像素之间，视差差别越大，代价越大，所以称之为平滑项。这种二维优化是一个NP-hard问题，常见的求近似解的方法有图割(Graph cut)，置信度传播(Belief Propagation)等。

Study of Energy Minimization Methods for Markov Random Fields with Smoothness-Based Priors 这篇论文对几种MRF的求解做了对比。

求解MRF的方法很多，无论是采用Graph cut还是BP, 亦或TRW-S，程序的性能都非常差。因此，出现了令外一种优化算法，半全局优化(SGM)。

半全局匹配(Semi-Global Matching)

半全局优化算法(SGM)将2维优化问题转换到1维，因此可以采用动态规划(DP)来优化。

SGM的原始论文，Stereo Processing by Semi-Global Matching and Mutual Information. Heiko Hirschmuller

这篇论文采用Mutual Information计算代价，再通过动态规划做代价聚合。之所以叫半全局，是因为他是在scanline上优化1D的能量函数，然后聚合相加。

抛开Mutual Information不谈，我们直接说他的代价聚合过程。

与全局匹配算法相同，SGM需要定义能量函数，

数据项(data term)与之前的定义的能量函数相同，使用代价表示匹配度。平滑项(Smooth term)引入了两个惩罚系数P1, P2。这两项惩罚相邻的像素，视差不同的情况，即相邻的像素，视差差异越大，惩罚越大。P1定义为小差异惩罚，即|Dp-Dq|为1的情况，P2定义视差差异大于1时的惩罚。

SGM对每一个像素在N个ScanLine上优化能量函数。N=4或者8，或者16。如下图，点P的scanline。

每一条Scan上，聚合函数定义如下，

这个公式一眼看上去有点啰嗦，举个例子就清楚了。

假设dmax=8，像素(x,y), d=4的某一个scanline上的聚合代价为，

$L_{r}(x,y,4)=C(x,y,4)+min\begin{Bmatrix} L_{r}(x-1,y,4)\\ {\color{Blue} L_{r}(x-1,y,5)+P_{1}}\\ {\color{Blue} L_{r}(x-1,y,3)+P_{1}}\\ {\color{Red} L_{r}(x-1,y,7)+P_{2}}\\ {\color{Red} L_{r}(x-1,y,6)+P_{2}}\\ {\color{Red} L_{r}(x-1,y,2)+P_{2}}\\ {\color{Red} L_{r}(x-1,y,1)+P_{2}}\\{\color{Red} L_{r}(x-1,y,0)+P_{2}}\\ \end{Bmatrix}- \min_{k=0...7} L_{r}(x-1,y,k)$

C(x, y, 4)是该点的代价，后面，取出前一个像素(x-1,y)的，

d=4时的聚合代价
与d=4相差1(3,5)的聚合代价加P1
与d=4相差大于1(7,6,2,1,0)的聚合代价加P2

取他们的最小值。再减掉L(x-1, y, k)， {k=0~7}的最小值。公式13和公式11是一致的。

将所有Scanline的聚合代价L相加，就得到像素最终的聚合代价。

对S，采用WTA策略，获得像素p的最终视差。

视差优化

亚像素差值(subpixel interpolation)

WTA算法挑选出来的视差是离散的整数值，为进一步提高精度，还需计算亚像素的视差。如下图，WTA给出的结果d=13, 在d=12到d=14之间，对cost做抛物线差值，得到亚像素视差d=12.8

对最终的视差图做滤波如中值滤波，双边滤波等
采用划分(segmentation)优化边缘
左右一致性检查排除错误视差(L-R consistency check)

通常T=1

等等等等。。。

SGM代码实例

https://github.com/linusyue/sgm.git

简单写了一份代码，放到了Github上。这份代码采用了SGM算法计算视差，并且集成了Kitti数据集的误差计算，根据数据集提供的GroundTruth计算了误差。原始SGM算法来自gishi523。

运行参数，

./sgm [dir]

dir为图像文件路径，需要按照Kitti数据集的结构放置源文件和Groundtruth。代码目录data下面放了几个例子。可以直接运行，

mkdir build
cd build && cmake ..
make./sgm ../data/training

程序流程

输入左右图像，首先对他们做census变换。代码提供了两种census的计算。第一种是CENSUS_9x7采用传统的计算census的方式。SYMMETRIC_CENSUS_9x7是在原始census算法上的简单改进。原始census比较的是中心像素与8邻域，而Symmetric census是让对称邻域两两比较。

    if (param_.censusType == CENSUS_9x7){census[0].create(h, w, CV_64F);census[1].create(h, w, CV_64F);census9x7<0>(I1, census[0]);census9x7<1>(I2, census[1]);int dir;OMP_PARALLEL_FORfor (dir = 0; dir < numDirections; dir++){L[dir].create(3, dims);scanCost<uint64_t>(census[0], census[1], L[dir], param_.P1, param_.P2, ru[dir], rv[dir]);}}else if (param_.censusType == SYMMETRIC_CENSUS_9x7){census[0].create(h, w, CV_32S);census[1].create(h, w, CV_32S);symmetricCensus9x7<0>(I1, census[0]);symmetricCensus9x7<1>(I2, census[1]);int dir;OMP_PARALLEL_FORfor (dir = 0; dir < numDirections; dir++){L[dir].create(3, dims);scanCost<uint32_t>(census[0], census[1], L[dir], param_.P1, param_.P2, ru[dir], rv[dir]);}}

一般来说，计算两幅图的census后，就可以计算相关像素的汉明距离来获取DSI空间了。这个程序考虑了执行效率的优化，把汉明计算这一步推迟到了聚合算法里去。scanCost计算某一路的聚合代价。所有路相加得到最终的聚合结果。

    template <typename T>
static void scanCost(const cv::Mat& C1, const cv::Mat& C2, cv::Mat1b& L, int P1, int P2, int ru, int rv)
{const int h = L.size[0];const int w = L.size[1];const int n = L.size[2];const bool forward = rv > 0 || (rv == 0 && ru > 0);int u0 = 0, u1 = w, du = 1, v0 = 0, v1 = h, dv = 1;if (!forward){u0 = w - 1; u1 = -1; du = -1;v0 = h - 1; v1 = -1; dv = -1;}for (int vc = v0; vc != v1; vc += dv){const T* _census1 = C1.template ptr<T>(vc);const T* _census2 = C2.template ptr<T>(vc) + w - 1;for (int uc = u0; uc != u1; uc += du){const int vp = vc - rv;const int up = uc - ru;const bool inside = vp >= 0 && vp < h && up >= 0 && up < w;uchar* _Lc = L.ptr<uchar>(vc, uc);uchar* _Lp = (uchar*)(L.data + vp * L.step.p[0] + up * L.step.p[1]); // for CV_DbgAssert avoidanceif (inside)updateCost(_census1[uc], _census2 - uc, _Lp, _Lc, uc, n, P1, P2);elseupdateCost(_census1[uc], _census2 - uc, _Lc, uc, n);}}
}

updataCost计算聚合值

template <typename T>
static inline void updateCost(T census1, const T* census2, const uchar* Lp, uchar* Lc, int u, int n, int P1, int P2)
{   uchar minLp = std::numeric_limits<uchar>::max();for (int d = 0; d < n; d++)minLp = std::min(minLp, Lp[d]);uchar _P1 = P1 - minLp;for (int d = 0; d < n; d++){const uchar MC = u - d >= 0 ? HammingDistance(census1, census2[d]) : DEFAULT_MC;const uchar Lp0 = Lp[d] - minLp;const uchar Lp1 = d > 0 ? Lp[d - 1] + _P1 : 0xFF;const uchar Lp2 = d < n - 1 ? Lp[d + 1] + _P1 : 0xFF;const uchar Lp3 = P2;Lc[d] = static_cast<uchar>(MC + min4(Lp0, Lp1, Lp2, Lp3));}
}

WTA获得视差，对视差图做中值滤波。再采用左右图检查，找出非法视差。

629     if (param_.pathType == SCAN_4PATH)
630         calcDisparity<WTA4Path>(L, S, D1, D2, param_.uniquenessRatio);
631     else
632         calcDisparity<WTA8Path>(L, S, D1, D2, param_.uniquenessRatio);
633
634     if (param_.medianKernelSize > 0)
635     {
636         cv::medianBlur(D1, D1, param_.medianKernelSize);
637         cv::medianBlur(D2, D2, param_.medianKernelSize);
638     }
639
640     const int max12Diff = param_.max12Diff << DISP_SHIFT;
641     if (max12Diff >= 0)
642     {
643         LRConsistencyCheck(D1, D2, max12Diff);
644     }

显示结果，黑色部分为LR-check找出的非法视差。

与GroundTruth的区别，保存在results/errors_disp_img_0/ 里.

下一篇想写一个立体视觉在实际中的应用案例。初步想法是聊一聊6D lab提出的Stixel.

计算机视觉大型攻略 —— 立体视觉（4）立体匹配算法简介与SGM相关推荐

计算机视觉大型攻略 —— 特征与匹配（3）特征描述符
接上一篇文章.这篇写特征描述符.特征匹配算法在确定角点后,还需要使用描述符来描述这些角点. 本文参考书籍:Computer Vision: Algorithms and Applications, R ...
计算机视觉大型攻略 —— 视觉里程计(1) 综述
参考文献: [1] Visual Odometry Part I: The First 30 Years and Fundamentals, Friedrich Fraundorfer and Da ...
计算机视觉大型攻略 —— SLAM(2) Graph-based SLAM（基于图优化的算法)
前面介绍了基于EKF的SLAM算法.EKF算法由于状态向量,协方差矩阵的大小随着特征点(路标)的增长而迅速增长,导致其不太适合大场景的应用.本文描述基于图优化的SLAM算法.目前由于SLAM图的稀疏性 ...
计算机视觉大型攻略 —— CUDA(2)执行模型
Professional CUDA C Programming[1]是一本不错的入门书籍,虽说命名为"Professional",但实际上非常适合入门阅读.他几乎涵盖了所有理论部分 ...
opencv计算机视觉编程攻略第2版,OpenCV计算机视觉编程攻略（第2版）pdf
摘要 1. 50多个知识点的案例解读,全面掌握基础知识与进阶内容 2. 学习OpenCV重要的图像操作类和函数 3. 初学者和从业者即查即用的工具书 4. 掌握计算机视觉与图像处理的基础知识与概念 O ...
定向光流直方图是什么_OpenCV计算机视觉编程攻略(第3版)
OpenCV计算机视觉编程攻略(第3版) 本书特色第3版简介内容速览阅读须知读者对象小标题排版规范读者反馈客户支持下载代码下载书中的彩色图片勘误举报盗版疑难解答电子书 1 ...
OpenCV计算机视觉编程攻略之行人检测
OpenCV计算机视觉编程攻略之行人检测,OpenCV 提供了一个基于HOG 和SVM且经过训练的行人检测器,可以用这个SVM 分类器以不同尺度的窗口扫描图像,在完整的图像中检测特定物体. 原图如下: ...
OpenCV计算机视觉编程攻略之提取图片轮廓-使用Canny函数
OpenCV计算机视觉编程攻略之提取图片轮廓-使用Canny函数,很方便..代码如下: #include <vector> #include <iostream> #inclu ...
OpenCV计算机视觉编程攻略之生成椒盐噪声实现
OpenCV计算机视觉编程攻略(第3版)P21的访问像素值,生成椒盐噪声实现. 运行结果图片,截图如下: 看书留下记录,代码如下: #include <random> #include & ...
万象物语找回服务器,万象物语新手大型攻略服务器、初始号的选择和新手前期需要做的事说明...
万象物语新手应该怎么玩?前期应该做什么呢?这里手机乐园寻隐者不遇小编来教教大家,咱们往下看! 万象物语新手大型攻略一.入坑须知不管是这个游戏哪方面吸引了你决定入坑,想长期玩下去并享受到游戏的乐趣, ...

计算机视觉大型攻略 —— 立体视觉（4）立体匹配算法简介与SGM