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

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

　　主要分四個部分講解：

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

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

一、攝像機標定

　　1.內參標定

　　攝像機內參反映的是攝像機座標系到影象座標系之間的投影關係。攝像機內參的標定使用張正友標定法，簡單易操作，具體原理請拜讀張正友的大作《A Flexible New Technique for Camera Calibration》。當然網上也會有很多資料可供查閱，MATLAB 有專門的攝像機標定工具包，OpenCV封裝好的攝像機標定API等。使用OpenCV進行攝像機標定的可以參考我的第一篇部落格：

　　攝像機的內參包括，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，此時停止濾波，得到最終結果。

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

 空洞填充的函式程式碼如下，僅供參考：

 1 void insertDepth32f(cv::Mat& depth)
 2 {
 3     const int width = depth.cols;
 4     const int height = depth.rows;
 5     float* data = (float*)depth.data;
 6     cv::Mat integralMap = cv::Mat::zeros(height, width, CV_64F);
 7     cv::Mat ptsMap = cv::Mat::zeros(height, width, CV_32S);
 8     double* integral = (double*)integralMap.data;
 9     int* ptsIntegral = (int*)ptsMap.data;
10     memset(integral, 0, sizeof(double) * width * height);
11     memset(ptsIntegral, 0, sizeof(int) * width * height);
12     for (int i = 0; i < height; ++i)
13     {
14         int id1 = i * width;
15         for (int j = 0; j < width; ++j)
16         {
17             int id2 = id1 + j;
18             if (data[id2] > 1e-3)
19             {
20                 integral[id2] = data[id2];
21                 ptsIntegral[id2] = 1;
22             }
23         }
24     }
25     // 積分割槽間
26     for (int i = 0; i < height; ++i)
27     {
28         int id1 = i * width;
29         for (int j = 1; j < width; ++j)
30         {
31             int id2 = id1 + j;
32             integral[id2] += integral[id2 - 1];
33             ptsIntegral[id2] += ptsIntegral[id2 - 1];
34         }
35     }
36     for (int i = 1; i < height; ++i)
37     {
38         int id1 = i * width;
39         for (int j = 0; j < width; ++j)
40         {
41             int id2 = id1 + j;
42             integral[id2] += integral[id2 - width];
43             ptsIntegral[id2] += ptsIntegral[id2 - width];
44         }
45     }
46     int wnd;
47     double dWnd = 2;
48     while (dWnd > 1)
49     {
50         wnd = int(dWnd);
51         dWnd /= 2;
52         for (int i = 0; i < height; ++i)
53         {
54             int id1 = i * width;
55             for (int j = 0; j < width; ++j)
56             {
57                 int id2 = id1 + j;
58                 int left = j - wnd - 1;
59                 int right = j + wnd;
60                 int top = i - wnd - 1;
61                 int bot = i + wnd;
62                 left = max(0, left);
63                 right = min(right, width - 1);
64                 top = max(0, top);
65                 bot = min(bot, height - 1);
66                 int dx = right - left;
67                 int dy = (bot - top) * width;
68                 int idLeftTop = top * width + left;
69                 int idRightTop = idLeftTop + dx;
70                 int idLeftBot = idLeftTop + dy;
71                 int idRightBot = idLeftBot + dx;
72                 int ptsCnt = ptsIntegral[idRightBot] + ptsIntegral[idLeftTop] - (ptsIntegral[idLeftBot] + ptsIntegral[idRightTop]);
73                 double sumGray = integral[idRightBot] + integral[idLeftTop] - (integral[idLeftBot] + integral[idRightTop]);
74                 if (ptsCnt <= 0)
75                 {
76                     continue;
77                 }
78                 data[id2] = float(sumGray / ptsCnt);
79             }
80         }
81         int s = wnd / 2 * 2 + 1;
82         if (s > 201)
83         {
84             s = 201;
85         }
86         cv::GaussianBlur(depth, depth, cv::Size(s, s), s, s);
87     }
88 }