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

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

之前有兩篇博客簡要講過OpenCV3.4中的兩種立體匹配算法效果比較。以及利用視差圖合成新視點。里面用到的匹配圖像對是OpenCV自帶校正好的圖像對。而目前大多數(shù)立體匹配算法使用的都是標(biāo)準(zhǔn)測試平臺提供的標(biāo)準(zhǔn)圖像對，比如著名的有如下兩個：
MiddleBury

KITTI

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

主要分四個部分講解：

攝像機標(biāo)定（包括內(nèi)參和外參）

雙目圖像的校正（包括畸變校正和立體校正）

立體匹配算法獲取視差圖，以及深度圖

利用視差圖，或者深度圖進(jìn)行虛擬視點的合成

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

一、攝像機標(biāo)定

1.內(nèi)參標(biāo)定

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

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

使用OpenCV3.4+VS2015對手機進(jìn)行內(nèi)參標(biāo)定。標(biāo)定結(jié)果如下，手機鏡頭不是魚眼鏡頭，因此使用普通相機模型標(biāo)定即可：

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

2.外參標(biāo)定

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

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

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

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

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

③ 利用特征點對pts1和pts2，以及內(nèi)參矩陣camK，解算出本質(zhì)矩陣E：

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

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

cv::Mat R1, R2; cv::decomposeEssentialMat(E, R1, R2, t); R = R1.clone(); t = -t.clone();

二、雙目圖像的校正

1. 畸變校正

畸變校正前面已經(jīng)介紹過，利用畸變系數(shù)進(jìn)行畸變校正即可，下面說一下立體校正。

2. 立體校正

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

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

② 得到上述參數(shù)后，就可以使用下面的API進(jìn)行對極線校正操作了，并將校正結(jié)果保存到本地：

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

對極線校正結(jié)果如下所示，查看對極線校正結(jié)果是否準(zhǔn)確，可以通過觀察若干對應(yīng)點是否在同一行上粗略估計得出：

三、立體匹配

1. SGBM算法獲取視差圖

立體校正后的左右兩幅圖像得到后，匹配點是在同一行上的，可以使用OpenCV中的BM算法或者SGBM算法計算視差圖。由于SGBM算法的表現(xiàn)要遠(yuǎn)遠(yuǎn)優(yōu)于BM算法，因此采用SGBM算法獲取視差圖。SGBM中的參數(shù)設(shè)置如下：

int numberOfDisparities = ((imgSize.width / 8)+ 15) & -16; cv::Ptr sgbm = cv::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::MODE_HH); else if (alg == STEREO_SGBM) sgbm->setMode(cv::MODE_SGBM); else if (alg == STEREO_3WAY) sgbm->setMode(cv::MODE_SGBM_3WAY);sgbm->compute(imgL, imgR, disp);

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

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

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

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

如果將右視差圖不可靠視差值也設(shè)置為0，則如下：

至此，左視差圖和右視差圖遙相呼應(yīng)。

2. 視差圖空洞填充

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

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

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

3. 視差圖轉(zhuǎn)換為深度圖

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

depth = ( f * baseline) / disp

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

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

/* 函數(shù)作用：視差圖轉(zhuǎn)深度圖 輸入： dispMap ----視差圖，8位單通道，CV_8UC1 K ----內(nèi)參矩陣，float類型 輸出： depthMap ----深度圖，16位無符號單通道，CV_16UC1 */void disp2Depth(cv::Mat dispMap, cv::Mat &depthMap, cv::Mat K) { int type = dispMap.type(); float fx = K.at(0, 0); float fy = K.at(1, 1); float cx = K.at(0, 2); float cy = K.at(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米的深度，足夠了。

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

轉(zhuǎn)換后的深度圖如下：

左深度圖 右深度圖

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

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

視差圖到深度圖完成。

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

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

步驟如下：

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

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

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

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

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::zeros(height,width, CV_64F); 7 cv::Mat ptsMap = cv::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 // 積分區(qū)間 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 }