【論文筆記】On Building an Accurate Stereo Matching System on Graphics Hardware

——同事推薦的一篇paper，想把它和SGBM opencv源碼學習筆記寫在一起，是因為感覺pipeline真的很像，可能他們都用了scanline optimization的原因吧 : )

前言

1. 文中指出，在middleBury上排名前10的，基本都是那種將多種重量級算法(CC, CA, optimization, refine)串聯起來，效果倒是好了，需要使用trick或是設計一些復雜的數據結構才能用GPU加速，在移動設備上使用，性能堪憂；

2. 作者將幾篇輕量級的模塊組合在了一起，可以在GPU上進行加速，效果也還不錯，使得在移動設備上使用成為了可能。其實，pipline也比較正常：AD Census + CBCA + Scanline Optimization + 常規Refine，分別引用自以下文獻：

• AD Census:
  X. Sun, X. Mei, S. Jiao, M. Zhou, and H. Wang. Stereo matching with reliable disparity propagation. In Proc. 3DIMPVT, pages 132–139, 2011.
• CBCA:
  K. Zhang, J. Lu, and G. Lafruit. Cross-based local stereo matching using orthogonal integral images. IEEE TCSVT, 19(7):1073–1079, 2009
• Scanline Optimization:
  H. Hirschm¨uller. Stereo processing by semiglobal matching and mutual information. IEEE TPAMI, 30(2):328–341, 2008.

摘要

paper鏈接

paper

源碼鏈接

code

核心思想

代價計算(Cost Computation, CC): AD + Census，即：the absolute differences (AD) measure and the census transform

代價聚合(Cost Aggregation, CA): CBCA，即基於crossed窗口的代價聚合方法

視差求解(disparity computation): Scanline optimization + WTA，很常規，引用於SGM的掃描線優化 + winner takes all

后處理(post process, refine)：區域投票 + 插值 + 邊緣優化 + 亞像素增強

貢獻點

作者構造了上述的pipeline，並基於cuda實現了GPU版本

相關研究

算法

AD Census

題外話：

R. Zabih and J. Woodfill. Non-parametric local transforms for computing visual correspondence. In Proc. ECCV, pages 151–158, 1994.

這篇paper里提到了一些代價計算的方法，包括AD, SGBM里邊用的BT，Census等方法，空的時候可以去讀讀；

H. Hirschm¨uller and D. Scharstein. Evaluation of stereo matching costs on images with radiometric differences.
IEEE TPAMI, 31(9):1582–1599, 2009.

H. Hirschm¨uller 大牛做的一個代價計算特征的評估，結論如下：

census shows the best overall results in local and global stereo matching methods

Census思想：

對於左圖中的點 \(\vec{p}= I_{left}(x, y)\)，在視差層級 \(d\) 下，\(C_{census}(\vec{p}, d)\) 定義為以左圖中點 \(\vec{p}\)為中心，整個鄰域窗口內的census編碼和以右圖中點為中心\(\vec{pd} = I_{right}(x-d, y)\)，鄰域窗口的census編碼之間的hamming距離，https://blog.csdn.net/qq_40313712/article/details/86349363 是一個比較簡介明了的解釋。

census優點——將鄰域點像素值和當前點的相對關系編碼起來，對噪聲以及畸變引入的outliers有較好的魯棒性

Census encodes local image structures with relative orderings of the pixel intensities other than the intensity values themselves, and therefore tolerates outliers due to radiometric changes and image noise

census不足——難以區分鄰域點和當前點較為相似的區域，諸如: 重復紋理場景；

However, this asset could also introduce matching ambiguities in image regions with repetitive or similar local structures.

為解決Census的問題，提出結合AD來使用：

AD定義如下：

\[C_{AD}(\vec{p}, d) = \frac{1}{3}\sum_{i=R, G, B}|I_i^{left}(\vec{p}) - I_i^{right}(\vec{pd})| \tag{1} \]

AD Census的定義如下：

\[C_{AD Census} = \rho(C_{census}(\vec{p}, d), \lambda_{census}) + \rho(C_{AD}(\vec{p}, d), \lambda_{AD}) \tag{2} \]

其中, \(\rho(c, \lambda)\)的定義如下：

\(\rho(c, \lambda) = 1- exp(-\frac{c}{\lambda}) \tag{3}\)

一個簡單對比如下圖：

CBCA

代價聚合可看做對代價立方體做refine的一個過程，也就是傳說中的聚合，用於降低噪聲對后面視差計算的影響，CBCA全稱是Cross-based Cost Aggregation，本質上就是一種cross型窗口的均值濾波，思想如下：

從上圖可以總結CBCA可分為兩步：

1. cross窗口構造
2. 窗口內的均值濾波

cross窗口構造

cross窗口，也就是十字窗口，會有自己的上下左右界限，左邊界限到當前點之間的距離，也就是文中說的"left arm"，上、下以及右亦如此；

每一個點都會有自己的cross窗口，用於實現下面的代價聚合；

現在問題來了，上下左右界限怎么確定，便是cross窗口構造要解決的核心問題。

縱觀可利用的資源，有以下幾個信息可以使用起來：

• 當前點和鄰域點之間的顏色距離 \(D_c(\vec{p_1}, \vec{p}) = max_{i=R,G,B}|I_i(\vec{p_1}) - I_i(\vec{p})|\)
• 當前點和鄰域點之間的空間距離 \(D_s(\vec{p_1}, \vec{p}) = |\vec{p1}-\vec{p}|\)

判斷一個鄰域點是否在當前點的cross窗口內，作者定義如下幾個准則：

1. \(D_c(\vec{p_1}, \vec{p}) < \tau_1\) 同時 \(D_c(\vec{p_1}, \vec{p_1} + (1, 0)) < \tau_1\)
2. \(D_s(\vec{p_1}, \vec{p}) < L1\)
3. 如果 \(L2<D_s(\vec{p_1}, \vec{p}) < L1\) 需滿足：\(D_c(\vec{p_1}, \vec{p}) < \tau_2\)

1是保證不跨邊緣做聚合，需要當前點和鄰域點之間的顏色距離小於\(\tau_1\)，同時，該鄰域點和它往外一點之間的顏色距離也要小於\(\tau_1\);

2以及3，是為了cover住無紋理場景足夠的平滑性，通過設定一個較大的\(L_1\)來對無紋理場景保證有足夠大的窗口，如果窗口半徑大於預設的半徑\(L2\)，則需滿足\(D_c(\vec{p_1}, \vec{p}) < \tau_2\)才可以。

窗口內的均值濾波

一些先驗知識：

• 均值濾波基於積分圖來實現可以達到O(n)的復雜度，並實現radius-free；

• 分離為一維均值濾波行列方向交替進行，會比直接算二維均值濾波快；

在上面構造好cross窗口后，構造代價立方體的積分圖，根據上一步計算得到的窗口半徑即可實現cross型窗口的代價聚合；

由於cross型窗口均值濾波並非常規均值濾波，在經歷行列交替的一維濾波后，通常會存在條狀artifacts，因此，便有了作者下邊的這部分迭代操作：

• 做4次迭代，一次迭代包含行方向和列方向；
• 1,3次迭代先行后列
• 2,4次迭代先列后行

引用一句作者的評價：

our improved method are presented in Figure 3, which shows that the enhanced cross construction rules and aggregation strategy can produce more accurate results in large textureless regions and near depth discontinuities.

​ 對於一個研究濾波起家的人來看，這簡直就是傳說中的保邊平滑性嘛，哈哈

Scanline Optimization

掃描線優化是從SGM原文里邊引用的，核心思想可借用SGM原文的公式來表述：

\[E(D)=\sum_\vec{p}(C(\vec{p}, D_p) + \sum_{\vec{q}\isin{N_p}}P1T[|D_p-D_q|=1]+\sum_{\vec{q}\isin{N_p}}P2T[|D_p-D_q|>1])\tag{4} \]

為了得到一個全局最優結果，通常會構建一個包含數據項和平滑項的能量函數，通過最小化能量函數值來求得最優視差解。代價立方體本身存在一個代價，對於鄰域內的像素點和當前點視差相差1的點做P1的懲罰，大於1的點做P2的懲罰。因此，第一項即為數據項，后面兩項則為平滑項。

但由於在立體匹配中，代價立方體維度都會很大，想要求解這個最優化問題，直接用傳統方法，內存和耗時都會很大。於是，SGM原文提出可使用動態規划求出近似解，代價函數簡化成了下式：

\[C_\vec{r}=C_1(\vec{p}, d) + min(C_\vec{r}(\vec{p}-\vec{r}, d), C_\vec{r}(\vec{p}-\vec{r}, d \pm 1) + P1, min_kC_\vec{r}(\vec{p} - \vec{r} , k) + P2) - min_kC_\vec{r}(\vec{p}-\vec{r}, k)\tag{5} \]

定義 \(D_1 = d_c(\vec{p}, \vec{p}-\vec{r})\)，\(D_2 = D_c(\vec{pd}, \vec{pd}-\vec{r})\), 通過\(D_1\)和\(D_2\)的值來控制\(P_1, P_2\)的值：

1. \(P_1 = \Pi_1, P_2 = \Pi_2, if D_1 < \tau_{SO}, D_2 < \tau_{SO}\)
2. \(P_1 = \Pi_1 / 4, P_2 = \Pi_2 / 4, if D_1 < \tau_{SO}, D_2 > \tau_{SO}\)
3. \(P_1 = \Pi_1 / 4, P_2 = \Pi_2 / 4, if D_1 > \tau_{SO}, D_2 < \tau_{SO}\)
4. \(P_1 = \Pi_1 / 10, P_2 = \Pi_2 / 10, if D_1 > \tau_{SO}, D_2 > \tau_{SO}\)

最終的代價立方體結果為: \(C_2(\vec{p}, d) = \frac{1}{4}\sum_\vec{r}C_\vec{r}(\vec{p}, d)\)

最終的視差圖用WTA即可得到，即：對每個像素點，取代價最小的那個視差值作為當前點視差

Post-process & Refine

經由WTA計算得到的視差圖通常在遮擋區域存在着無效視差，同時也會存在一些噪聲點甚至是大面積出錯的情況，作者用了以下幾個模塊來做refine：

1. Outlier Detection

   • 這部分主要是基於左右視差圖的一致性檢測來提取出視差不准確區域，簡單而粗暴的設一個閾值，左右圖匹配點對應的左右視差大於某個閾值，則該點即為視差不准確點；

   • 視差不准確點可分為兩種：1）遮擋點；2）誤匹配點

   • 對於上述兩種點的區分，作者使用的是SGM原文里的方法：
     1） 在極線上搜索，若在右視差圖上找到相交的點，則該點為誤匹配點

     2） 在極線上右視差圖上沒有搜索到相交的點，則該點為遮擋點

2. Iterative Region Voting

   對於上一步檢測出來的outlier，統計其cross窗口中可靠點的視差直方圖\(H_\vec{p}\)，分為\(d_{max}+1\)個bins，對於可靠點數最多的bin即為\(d_\vec{p}^*\)，記所有可靠像素點數為\(S_\vec{p} = \sum_{d=0}^{d=d_{max}}H_\vec{p}(d)\)，當\(H_\vec{p}\)滿足以下條件時：

   \(S_\vec{p} > \tau_S, \frac{H_\vec{p}(d_\vec{p}^*)}{S_\vec{p}} > \tau_H\)

   outlier處的視差由\(d_\vec{p}^*\)替代，五次迭代后，錯誤區域會少很多

3. Proper Interpolation

   經由上一步之后，錯誤點會少很多。這一步會對誤匹配點&遮擋點區別對待；

   搜尋outliter周圍最相近的16個方向上的可靠點；

   對於遮擋點，選取其中最小的視差最為其視差值；

   對於誤匹配點，選取和其最相似的那個點的視差作為其視差值；

4. Depth Discontinuity Adjustment

   Step 1: 檢測視差圖的邊緣

   Step 2: 判斷邊界點兩邊點的代價值，如果低於當前點的代價值，則調整當前邊界點的視差至邊界兩邊的那個點

5. Sub-pixel Enhancement

   ——基於quadratic polynomial interpolation的亞像素增強方法

后面再跟一個3x3的中值濾波

refine各個stage帶來的收益：

實驗

測試平台

PC with Core2Duo 2.20GHz CPU and NVIDIA GeForce GTX 480 graphics card

數據集

Middlebury

參數設置

性能

	Tsukuba	Venus	Teddy	Cones
CPU	2.5s	4.5s	15s	15s
GPU	0.016s	0.032s	0.095s	0.094s

tips

調參很重要

后記——SGBM OpenCV源碼學習筆記

在opencv源碼的位置:

opencv-4.0.0_src\opencv-4.0.0\modules\calib3d\src\stereosgbm.cpp

推薦博客：

https://blog.csdn.net/wwp2016/article/details/86080722

https://zhuanlan.zhihu.com/p/53060518

pipeline

pipeline很符合常規立體匹配pipeline：代價計算 -> 代價聚合 -> 視差計算 -> 視差后處理，和上面paper的邏輯極為相近，可以串起來做相關優化和嘗試，所以把他們寫在了一篇文檔里。

函數調用與參數解析

SGBM opencv源碼將上述pipeline封裝為兩個接口供外界調用：1）通過帶參數設置的構造函數設置參數；2）算法核心；

    StereoSGBMImpl()
    {
        params = StereoSGBMParams();
    }

    StereoSGBMImpl( int _minDisparity, int _numDisparities, int _SADWindowSize,
                    int _P1, int _P2, int _disp12MaxDiff, int _preFilterCap,
                    int _uniquenessRatio, int _speckleWindowSize, int _speckleRange,
                    int _mode )
    {
        params = StereoSGBMParams( _minDisparity, _numDisparities, _SADWindowSize,
                                   _P1, _P2, _disp12MaxDiff, _preFilterCap,
                                   _uniquenessRatio, _speckleWindowSize, _speckleRange,
                                   _mode );
        // _minDisparity: 最小視差，默認為0; 對於實際應用中可能存在負視差的場景，可適當是設置為負視差，e.g. -3；
        
        // _numDisparities: 視差層級，需設置為16的倍數，默認為16；
       	//				  注意：輸出視差圖的視差層級為(16 * _minDisparity + _numDisparities)；
        
        // _SADWindowSize: SAD窗口直徑，代價聚合部分，形同boxfilter的窗口直徑，默認為5；
        
        // _P1: 鄰域點視差和當前點相差1的懲罰量，即上面eq4的P1，默認為2；
        
        // _P2: 鄰域點視差和當前點相差大於1的懲罰量，即上面eq4的P2， 默認：max(params.P2 > 0 ? params.P2 : 5, P1+1);
        
        // _disp12MaxDiff: LR Check中不匹配點檢測的閾值，默認：1；
        
        // _preFilterCap: SGBM 代價計算是基於BT特征的，這個值是將每個匹配點對的匹配代價值抑制在_preFilterCap - 2 * _preFilterCap范圍內；
        //                 注意：該值會影響到匹配代價的數量級，P1， P2的調參需要和該值協同修改；
        //                 BT的大概思想：
        //                 if k < 256*4 - _preFilterCap, 則CC = _preFilterCap；
        //                 if 256*4 - _preFilterCap < k < 256*4 + _preFilterCap, 則CC = k - 256*4 + _preFilterCap；
        //                 if k > 256*4 + _preFilterCap, 則CC = 2 * _preFilterCap;
        
        // _uniquenessRatio: 唯一匹配檢測閾值。要求除了bestDisp前后各一個視差之外，其余視差值對應的S必須大於minS * 1.x，默認為10, 大於1.1；
        
        // _speckleWindowSize & _speckleRange: filterSpeckle的參數，用於對視差對做斑點去除的函數;
        
        // _mode: 
        //       - MODE_SGBM_3WAY:
        //       - MODE_HH4:
        //       - MODE_HH：默認模式
    }

    void compute( InputArray leftarr, InputArray rightarr, OutputArray disparr ) CV_OVERRIDE
    {
        CV_INSTRUMENT_REGION();

        Mat left = leftarr.getMat(), right = rightarr.getMat();
        CV_Assert( left.size() == right.size() && left.type() == right.type() &&
                   left.depth() == CV_8U );

        disparr.create( left.size(), CV_16S );
        Mat disp = disparr.getMat();

        if(params.mode==MODE_SGBM_3WAY)
            computeDisparity3WaySGBM( left, right, disp, params, buffers, num_stripes );
        else if(params.mode==MODE_HH4)
            computeDisparitySGBM_HH4( left, right, disp, params, buffer );
        else
            computeDisparitySGBM( left, right, disp, params, buffer );

        medianBlur(disp, disp, 3);

        if( params.speckleWindowSize > 0 )
            filterSpeckles(disp, (params.minDisparity - 1)*StereoMatcher::DISP_SCALE, params.speckleWindowSize,
                           StereoMatcher::DISP_SCALE*params.speckleRange, buffer);
    }

一些代碼細節

SGBM的三個模式

1. MODE_SGBM_3WAY
2. MODE_SGBM_HH
3. MODE_SGBM_HH4

SGBM的內存管理

核心算法思想

1. 基於BT的代價計算

   ——Birchfield-Tomasi metric

   from BT org paper

   最初是在A Pixel Dissimilarity Measure That Is Insensitive to Image Sampling中由Stan Birchfield and Carlo Tomasi提出，旨在提升相似性度量對image sampling的魯棒，BT在該paper的Section 2中提出了BT的基本思路，在Section 3從理論上和實驗中證明了該算法的有效性，具體的內容沒有很清楚，但對BT的主要思想做了下梳理，內容如下：

   對於兩幅連續的立體對\(I_L(x_L)\), \(I_R(x_R)\)，經過硬件處理、圖像采樣后，左圖和右圖匹配點的灰度值可能會出現一些偏差，而BT的思想就是基於簡單的線性插值，搜索\(x+1, x-1\)的\(I_{min}，I_{max}\)來協同表征當前位置下左右視圖的相似度：

數學表達如下：

設\(I_R^- \equiv \hat{I_R}(x_R - \frac{1}{2}) = \frac{1}{2}(I_R(x_R) + I_R(x_R - 1))\)， \(I_R^+ \equiv \hat{I_R}(x_R + \frac{1}{2}) = \frac{1}{2}(I_R(x_R) + I_R(x_R + 1))\)

設\(I_{min} = \min{I_R^-, I_R^+, I_R(x_R)}\), \(I_{max} = \max{I_R^-, I_R^+, I_R(x_R)}\)

則cost為\(\bar{d}(x_L, x_R, I_L, I_R) = \max{\{0, I_L(x_L) - I_{max}, I_{min} - I_L(x_L)\}}\)

下圖為BT vs AD的一個具體的例子：

SGBM中的BT

2. 基於Scanline Optimization的代價優化

   ——代價聚合鄰域點的取值方式

調用demo

#include "stdafx.h"
#include "opencv2/opencv.hpp
using namespace std;
using namespace cv;
int _tmain(int argc, _TCHAR* argv[])
{
	Mat left = imread("imgL.jpg", IMREAD_GRAYSCALE);
	Mat right = imread("imgR.jpg", IMREAD_GRAYSCALE);
	Mat disp;
	int mindisparity = 0;
	int ndisparities = 64;  
	int SADWindowSize = 11; 
	//SGBM
	cv::Ptr<cv::StereoSGBM> sgbm = cv::StereoSGBM::create(mindisparity, ndisparities, SADWindowSize);
	int P1 = 8 * left.channels() * SADWindowSize* SADWindowSize;
	int P2 = 32 * left.channels() * SADWindowSize* SADWindowSize;
	sgbm->setP1(P1);
	sgbm->setP2(P2);
	sgbm->setPreFilterCap(15);
	sgbm->setUniquenessRatio(10);
	sgbm->setSpeckleRange(2);
	sgbm->setSpeckleWindowSize(100);
	sgbm->setDisp12MaxDiff(1);
	//sgbm->setMode(cv::StereoSGBM::MODE_HH);
	sgbm->compute(left, right, disp);
	disp.convertTo(disp, CV_32F, 1.0 / 16);                //除以16得到真實視差值
	Mat disp8U = Mat(disp.rows, disp.cols, CV_8UC1);       //顯示
	normalize(disp, disp8U, 0, 255, NORM_MINMAX, CV_8UC1);
	imwrite("results/SGBM.jpg", disp8U);
	return 0;
}