?實戰(zhàn) | 在真實場景實現(xiàn)雙目立體匹配獲取深度圖
點擊左上方藍字關注我們

目前大多數(shù)立體匹配算法使用的都是標準測試平臺提供的標準圖像對,比如著名的有如下兩個:
MiddleBury: http://vision.middlebury.edu/stereo/;
KITTI:http://www.cvlibs.net/datasets/kitti/eval_scene_flow.php?benchmark=stereo
但是對于想自己嘗試拍攝雙目圖片進行立體匹配獲取深度圖,進行三維重建等操作的童鞋來講,要做的工作是比使用校正好的標準測試圖像對要多的。因此博主覺得有必要從用雙目相機拍攝圖像開始,捋一捋這整個流程。
主要分四個部分講解:
攝像機標定(包括內(nèi)參和外參)
雙目圖像的校正(包括畸變校正和立體校正)
立體匹配算法獲取視差圖,以及深度圖
利用視差圖,或者深度圖進行虛擬視點的合成
注:如果沒有雙目相機,可以使用單個相機平行移動拍攝,外參可以通過攝像機自標定算出。我用自己的手機拍攝,拍攝移動時盡量保證平行移動。
一、攝像機標定
1.內(nèi)參標定
攝像機內(nèi)參反映的是攝像機坐標系到圖像坐標系之間的投影關系。攝像機內(nèi)參的標定使用張正友標定法,簡單易操作,具體原理請拜讀張正友的大作《A Flexible New Technique for Camera Calibration》。當然網(wǎng)上也會有很多資料可供查閱,MATLAB 有專門的攝像機標定工具包,OpenCV封裝好的攝像機標定API等。
攝像機的內(nèi)參包括,fx, fy, cx, cy,以及畸變系數(shù)[k1,k2,p1,p2,k3],詳細就不贅述。我用手機對著電腦拍攝各個角度的棋盤格圖像,棋盤格圖像如圖所示:

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

圖像分辨率為:3968 x 2976。上面標定結果順序依次為fx, fy, cx, cy, k1, k2, p1, p2, k3, 保存到文件中供后續(xù)使用。
2.外參標定
攝像機外參反映的是攝像機坐標系和世界坐標系之間的旋轉R和平移T關系。如果兩個相機的內(nèi)參均已知,并且知道各自與世界坐標系之間的R1、T1和R2,T2,就可以算出這兩個相機之間的Rotation和Translation,也就找到了從一個相機坐標系到另一個相機坐標系之間的位置轉換關系。攝像機外參標定也可以使用標定板,只是保證左、右兩個相機同時拍攝同一個標定板的圖像。外參一旦標定好,兩個相機的結構就要保持固定,否則外參就會發(fā)生變化,需要重新進行外參標定。
那么手機怎么保證拍攝同一個標定板圖像并能夠保持相對位置不變,這個是很難做到的,因為后續(xù)用來拍攝實際測試圖像時,手機的位置肯定會發(fā)生變化。因此我使用外參自標定的方法,在拍攝實際場景的兩張圖像時,進行攝像機的外參自標定,從而獲取當時兩個攝像機位置之間的Rotation和Translation。
比如:我拍攝這樣兩幅圖像,以后用來進行立體匹配和虛擬視點合成的實驗。





cv::Mat E = cv::findEssentialMat(tmpPts1, tmpPts2, camK, CV_RANSAC);cv::Mat R1, R2;cv::decomposeEssentialMat(E, R1, R2, t);R = R1.clone();t = -t.clone();
cv::stereoRectify(camK, D, camK, D, imgL.size(), R, -R*t, R1, R2, P1, P2, Q);
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);

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);
sgbm->setMinDisparity(-numberOfDisparities);sgbm->setNumDisparities(numberOfDisparities);sgbm->compute(imgR, imgL, disp);disp = abs(disp);




depth = ( f * baseline) / disp
/*函數(shù)作用:視差圖轉深度圖輸入: 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<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; //基線距離65mmif (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);}}



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;}}}// 積分區(qū)間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);}}
END
整理不易,點贊鼓勵一下吧↓
