真實場景的雙目立體匹配（Stereo Matching）獲取深度圖詳解

　　雙目立體匹配一直是雙目視覺的研究熱點，雙目相機拍攝同一場景的左、右兩幅視點圖像，運用立體匹配匹配算法獲取視差圖，進而獲取深度圖。而深度圖的應用范圍非常廣泛，由於其能夠記錄場景中物體距離攝像機的距離，可以用以測量、三維重建、以及虛擬視點的合成等。

　　之前有兩篇博客簡要講過OpenCV3.4中的兩種立體匹配算法效果比較：http://www.cnblogs.com/riddick/p/8318997.html 。以及利用視差圖合成新視點： http://www.cnblogs.com/riddick/p/7355353.html。里面用到的匹配圖像對是OpenCV自帶校正好的圖像對。而目前大多數立體匹配算法使用的都是標准測試平台提供的標准圖像對，比如著名的有如下兩個：
　　MiddleBury: http://vision.middlebury.edu/stereo/；

　　KITTI：http://www.cvlibs.net/datasets/kitti/eval_scene_flow.php?benchmark=stereo。

　　但是對於想自己嘗試拍攝雙目圖片進行立體匹配獲取深度圖，進行三維重建等操作的童鞋來講，要做的工作是比使用校正好的標准測試圖像對要多的。因此博主覺得有必要從用雙目相機拍攝圖像開始，捋一捋這整個流程。

　　主要分四個部分講解：

• 攝像機標定（包括內參和外參）
• 雙目圖像的校正（包括畸變校正和立體校正）
• 立體匹配算法獲取視差圖，以及深度圖
• 利用視差圖，或者深度圖進行虛擬視點的合成

　　

　　注：如果沒有雙目相機，可以使用單個相機平行移動拍攝，外參可以通過攝像機自標定算出。我用自己的手機拍攝，拍攝移動時盡量保證平行移動。

一、攝像機標定

　　1.內參標定

　　攝像機內參反映的是攝像機坐標系到圖像坐標系之間的投影關系。攝像機內參的標定使用張正友標定法，簡單易操作，具體原理請拜讀張正友的大作《A Flexible New Technique for Camera Calibration》。當然網上也會有很多資料可供查閱，MATLAB 有專門的攝像機標定工具包，OpenCV封裝好的攝像機標定API等。使用OpenCV進行攝像機標定的可以參考我的第一篇博客：http://www.cnblogs.com/riddick/p/6696858.html。里面提供有張正友標定法OpenCV實現的源代碼git地址，僅供參考。

　　攝像機的內參包括，fx, fy, cx, cy,以及畸變系數[k1,k2,p1,p2,k3]，詳細就不贅述。我用手機對着電腦拍攝各個角度的棋盤格圖像，棋盤格圖像如圖所示：

　　使用OpenCV3.4+VS2015對手機進行內參標定。標定結果如下，手機鏡頭不是魚眼鏡頭，因此使用普通相機模型標定即可：

　　圖像分辨率為：3968 x 2976。上面標定結果順序依次為fx, fy, cx, cy,   k1, k2, p1, p2, k3， 保存到文件中供后續使用。

　　2.外參標定

　　攝像機外參反映的是攝像機坐標系和世界坐標系之間的旋轉R和平移T關系。如果兩個相機的內參均已知，並且知道各自與世界坐標系之間的R1、T1和R2，T2，就可以算出這兩個相機之間的Rotation和Translation，也就找到了從一個相機坐標系到另一個相機坐標系之間的位置轉換關系。攝像機外參標定也可以使用標定板，只是保證左、右兩個相機同時拍攝同一個標定板的圖像。外參一旦標定好，兩個相機的結構就要保持固定，否則外參就會發生變化，需要重新進行外參標定。

　　那么手機怎么保證拍攝同一個標定板圖像並能夠保持相對位置不變，這個是很難做到的，因為后續用來拍攝實際測試圖像時，手機的位置肯定會發生變化。因此我使用外參自標定的方法，在拍攝實際場景的兩張圖像時，進行攝像機的外參自標定，從而獲取當時兩個攝像機位置之間的Rotation和Translation。

　　比如：我拍攝這樣兩幅圖像，以后用來進行立體匹配和虛擬視點合成的實驗。

  

　① 利用攝像機內參進行畸變校正，手機的畸變程度都很小，校正后的兩幅圖如下：

  

　　② 將上面兩幅畸變校正后的圖作為輸入，使用OpenCV中的光流法提取匹配特征點對，pts1和pts2，在圖像中畫出如下：

  

　　③ 利用特征點對pts1和pts2，以及內參矩陣camK，解算出本質矩陣E：

    cv::Mat E = cv::findEssentialMat(tmpPts1, tmpPts2, camK, CV_RANSAC);

　　④  利用本質矩陣E解算出兩個攝像機之間的Rotation和Translation，也就是兩個攝像機之間的外參。以下是OpenCV中API函數實現的，具體請參見API文檔：

    cv::Mat R1, R2;
    cv::decomposeEssentialMat(E, R1, R2, t);

    R = R1.clone();
    t = -t.clone();

 

二、雙目圖像的校正

　　1. 畸變校正

　　畸變校正前面已經介紹過，利用畸變系數進行畸變校正即可，下面說一下立體校正。

　　2. 立體校正

　　① 得到兩個攝像機之間的 Rotation和Translation之后，要用下面的API對兩幅圖像進行立體對極線校正，這就需要算出兩個相機做對極線校正需要的R和T，用R1,T1, R2, T2表示，以及透視投影矩陣P1，P2：

cv::stereoRectify(camK, D, camK, D, imgL.size(), R, -R*t,  R1, R2, P1, P2, Q);

 　　② 得到上述參數后，就可以使用下面的API進行對極線校正操作了，並將校正結果保存到本地：

　　cv::initUndistortRectifyMap(P1(cv::Rect(0, 0, 3, 3)), D, R1, P1(cv::Rect(0, 0, 3, 3)), imgL.size(), CV_32FC1, mapx, mapy);
    cv::remap(imgL, recImgL, mapx, mapy, CV_INTER_LINEAR);
    cv::imwrite("data/recConyL.png", recImgL);

    cv::initUndistortRectifyMap(P2(cv::Rect(0, 0, 3, 3)), D, R2, P2(cv::Rect(0, 0, 3, 3)), imgL.size(), CV_32FC1, mapx, mapy);
    cv::remap(imgR, recImgR, mapx, mapy, CV_INTER_LINEAR);
    cv::imwrite("data/recConyR.png", recImgR);

 　　對極線校正結果如下所示，查看對極線校正結果是否准確，可以通過觀察若干對應點是否在同一行上粗略估計得出：

  

 

---

 

 

三、立體匹配

 　　1. SGBM算法獲取視差圖

　　立體校正后的左右兩幅圖像得到后，匹配點是在同一行上的，可以使用OpenCV中的BM算法或者SGBM算法計算視差圖。由於SGBM算法的表現要遠遠優於BM算法，因此采用SGBM算法獲取視差圖。SGBM中的參數設置如下：

　　 int numberOfDisparities = ((imgSize.width / 8) + 15) & -16;
    cv::Ptr<cv::StereoSGBM> sgbm = cv::StereoSGBM::create(0, 16, 3);
    sgbm->setPreFilterCap(32);
    int SADWindowSize = 9;
    int sgbmWinSize = SADWindowSize > 0 ? SADWindowSize : 3;
    sgbm->setBlockSize(sgbmWinSize);
    int cn = imgL.channels();
    sgbm->setP1(8 * cn*sgbmWinSize*sgbmWinSize);
    sgbm->setP2(32 * cn*sgbmWinSize*sgbmWinSize);
    sgbm->setMinDisparity(0); sgbm->setNumDisparities(numberOfDisparities);
    sgbm->setUniquenessRatio(10);
    sgbm->setSpeckleWindowSize(100);
    sgbm->setSpeckleRange(32);
    sgbm->setDisp12MaxDiff(1);
    int alg = STEREO_SGBM;
    if (alg == STEREO_HH)
        sgbm->setMode(cv::StereoSGBM::MODE_HH);
    else if (alg == STEREO_SGBM)
        sgbm->setMode(cv::StereoSGBM::MODE_SGBM);
    else if (alg == STEREO_3WAY)
        sgbm->setMode(cv::StereoSGBM::MODE_SGBM_3WAY);
    sgbm->compute(imgL, imgR, disp);

　　默認計算出的是左視差圖，如果需要計算右視差圖，則將上面加粗的三條語句替換為下面前三條語句。由於視差值計算出來為負值，disp類型為16SC1，因此需要取絕對值，然后保存：

　　 sgbm->setMinDisparity(-numberOfDisparities); sgbm->setNumDisparities(numberOfDisparities); sgbm->compute(imgR, imgL, disp);

    disp = abs(disp);

　　SGBM算法得到的左、右視差圖如下，左視差圖的數據類型為CV_16UC1，右視差圖的數據類型為CV_16SC1 （SGBM中視差圖中不可靠的視差值設置為最小視差（mindisp-1）*16。因此在此例中，左視差圖中不可靠視差值設置為-16，截斷值為0；右視差圖中不可靠視差值設置為（-numberOfDisparities-1）*16，取絕對值后為（numberOfDisparities+1）*16，所以兩幅圖會有較大差別）：

左視差圖（不可靠視差值為0）                                                      右視差圖（不可靠視差值為 （numberOfDisparities+1）*16 ）

  

　　如果將右視差圖不可靠視差值也設置為0，則如下

　　至此，左視差圖和右視差圖遙相呼應。

　　2. 視差圖空洞填充

　　視差圖中視差值不可靠的視差大多數是由於遮擋引起，或者光照不均勻引起。既然牛逼如SGBM也覺得不可靠，那與其留着做個空洞，倒不如用附近可靠的視差值填充一下。

　　空洞填充也有很多方法，在這里我檢測出空洞區域，然后用附近可靠視差值的均值進行填充。填充后的視差圖如下：

填充后左視差圖                                                                             填充后右視差圖

  

 

---

 

　　3. 視差圖轉換為深度圖

　　視差的單位是像素（pixel），深度的單位往往是毫米（mm）表示。而根據平行雙目視覺的幾何關系（此處不再畫圖推導，很簡單），可以得到下面的視差與深度的轉換公式：

 depth = ( f * baseline) / disp

 　　上式中，depth表示深度圖；f表示歸一化的焦距，也就是內參中的fx； baseline是兩個相機光心之間的距離，稱作基線距離；disp是視差值。等式后面的均已知，深度值即可算出。

　　

　　在上面我們用SGBM算法獲取了視差圖，接下來轉換為深度圖，函數代碼如下：

/*
函數作用：視差圖轉深度圖
輸入：
　　dispMap ----視差圖，8位單通道，CV_8UC1
　　K       ----內參矩陣，float類型
輸出：
　　depthMap ----深度圖，16位無符號單通道，CV_16UC1
*/

void disp2Depth(cv::Mat dispMap, cv::Mat &depthMap, cv::Mat K)
{
    int type = dispMap.type();

    float fx = K.at<float>(0, 0);
    float fy = K.at<float>(1, 1);
    float cx = K.at<float>(0, 2);
    float cy = K.at<float>(1, 2);
    float baseline = 65; //基線距離65mm if (type == CV_8U)
    {
        const float PI = 3.14159265358;
        int height = dispMap.rows;
        int width = dispMap.cols;

        uchar* dispData = (uchar*)dispMap.data;
        ushort* depthData = (ushort*)depthMap.data;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                int id = i*width + j;
                if (!dispData[id])  continue;  //防止0除
                depthData[id] = ushort( (float)fx *baseline / ((float)dispData[id]) );
            }
        }
    }
    else
    {
        cout << "please confirm dispImg's type!" << endl;
        cv::waitKey(0);
    }
}

 

　　注：png的圖像格式可以保存16位無符號精度，即保存范圍為0-65535，如果是mm為單位，則最大能表示約65米的深度，足夠了。

　　上面代碼中我設置深度圖的精度為CV_16UC1，也就是ushort類型，將baseline設置為65mm，轉換后保存為png格式即可。如果保存為jpg或者bmp等圖像格式，會將數據截斷為0-255。所以保存深度圖，png格式是理想的選擇。（如果不是為了獲取精確的深度圖，可以將baseline設置為1，這樣獲取的是相對深度圖，深度值也是相對的深度值）

 　　轉換后的深度圖如下：

左深度圖                                                                                        右深度圖 

  

　　空洞填充后的深度圖，如下：

左深度圖（空洞填充后）                                                                右深度圖（空洞填充后）

  

　　視差圖到深度圖完成。

 

　　注：視差圖和深度圖中均有計算不正確的點，此文意在介紹整個流程，不特別注重算法的優化，如有大神望不吝賜教。

---

 

附：視差圖和深度圖的空洞填充

 　　步驟如下：

　　① 以視差圖dispImg為例。計算圖像的積分圖integral，並保存對應積分圖中每個積分值處所有累加的像素點個數n（空洞處的像素點不計入n中，因為空洞處像素值為0，對積分值沒有任何作用，反而會平滑圖像）。

　　② 采用多層次均值濾波。首先以一個較大的初始窗口去做均值濾波（積分圖實現均值濾波就不多做介紹了，可以參考我之前的一篇博客），將大區域的空洞賦值。然后下次濾波時，將窗口尺寸縮小為原來的一半，利用原來的積分圖再次濾波，給較小的空洞賦值（覆蓋原來的值）；依次類推，直至窗口大小變為3x3，此時停止濾波，得到最終結果。

　　③ 多層次濾波考慮的是對於初始較大的空洞區域，需要參考更多的鄰域值，如果采用較小的濾波窗口，不能夠完全填充，而如果全部采用較大的窗口，則圖像會被嚴重平滑。因此根據空洞的大小，不斷調整濾波窗口。先用大窗口給所有空洞賦值，然后利用逐漸變成小窗口濾波覆蓋原來的值，這樣既能保證空洞能被填充上，也能保證圖像不會被過度平滑。

 

 空洞填充的函數代碼如下，僅供參考：

void insertDepth32f(cv::Mat& depth)
{
    const int width = depth.cols;
    const int height = depth.rows;
    float* data = (float*)depth.data;
    cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
    cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
    double* integral = (double*)integralMap.data;
    int* ptsIntegral = (int*)ptsMap.data;
    memset(integral, 0, sizeof(double) * width * height);
    memset(ptsIntegral, 0, sizeof(int) * width * height);
    for (int i = 0; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 0; j < width; ++j)
        {
            int id2 = id1 + j;
            if (data[id2] > 1e-3)
            {
                integral[id2] = data[id2];
                ptsIntegral[id2] = 1;
            }
        }
    }
    // 積分區間
    for (int i = 0; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 1; j < width; ++j)
        {
            int id2 = id1 + j;
            integral[id2] += integral[id2 - 1];
            ptsIntegral[id2] += ptsIntegral[id2 - 1];
        }
    }
    for (int i = 1; i < height; ++i)
    {
        int id1 = i * width;
        for (int j = 0; j < width; ++j)
        {
            int id2 = id1 + j;
            integral[id2] += integral[id2 - width];
            ptsIntegral[id2] += ptsIntegral[id2 - width];
        }
    }
    int wnd;
    double dWnd = 2;
    while (dWnd > 1)
    {
        wnd = int(dWnd);
        dWnd /= 2;
        for (int i = 0; i < height; ++i)
        {
            int id1 = i * width;
            for (int j = 0; j < width; ++j)
            {
                int id2 = id1 + j;
                int left = j - wnd - 1;
                int right = j + wnd;
                int top = i - wnd - 1;
                int bot = i + wnd;
                left = max(0, left);
                right = min(right, width - 1);
                top = max(0, top);
                bot = min(bot, height - 1);
                int dx = right - left;
                int dy = (bot - top) * width;
                int idLeftTop = top * width + left;
                int idRightTop = idLeftTop + dx;
                int idLeftBot = idLeftTop + dy;
                int idRightBot = idLeftBot + dx;
                int ptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] - (ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
                double sumGray = integral[idRightBot] + integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
                if (ptsCnt <= 0)
                {
                    continue;
                }
                data[id2] = float(sumGray / ptsCnt);
            }
        }
        int s = wnd / 2 * 2 + 1;
        if (s > 201)
        {
            s = 201;
        }
        cv::GaussianBlur(depth, depth, cv::Size(s, s), s, s);
    }
}