真實場景的雙目立體匹配（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、攝像機標定

　　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關係。若是兩個相機的內參均已知，而且知道各自與世界座標系之間的R一、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();

 

2、雙目圖像的校訂

　　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);

 　　對極線校訂結果以下所示，查看對極線校訂結果是否準確，能夠經過觀察若干對應點是否在同一行上粗略估計得出：

  

 

---

 

 

3、立體匹配

 　　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 }

View Code