计算机视觉大型攻略 —— 立体视觉（4）立体匹配算法简介与SGM_立体匹配中半全局和全局的区别

本文是立体视觉部分的第四篇，立体匹配。主要介绍了立体匹配的算法思路，详细介绍了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 给出了立体匹配算法的几个步骤，

1. 计算匹配的代价(Matching cost)
2. 代价聚合(Cost aggregation)
3. 视差计算
4. 视差优化

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

代价计算

无论什么样的算法，第一步都需要定义像素和像素之间的相关性，也就是匹配代价(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) 算法为例，他的算法流程如下，

1. 通过计算灰度的平方差，计算像素所有d的代价（DSI）。
2. 在一个固定窗口内，为该像素的每个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}$ 
3. 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是像素集合。$\tiny \sum_{p\in R}C(p,d)$是数据项，定义了像素代价的匹配程度。$\tiny \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_FOR
            for (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_FOR
            for (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 avoidance
 
            if (inside)
                updateCost(_census1[uc], _census2 - uc, _Lp, _Lc, uc, n, P1, P2);
            else
                updateCost(_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.