真实场景的双目立体匹配（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，此时停止滤波，得到最终结果。
③ 多层次滤波考虑的是对于初始较大的空洞区域，需要参考更多的邻域值，如果采用较小的滤波窗口，不能够完全填充，而如果全部采用较大的窗口，则图像会被严重平滑。因此根据空洞的大小，不断调整滤波窗口。先用大窗口给所有空洞赋值，然后利用逐渐变成小窗口滤波覆盖原来的值，这样既能保证空洞能被填充上，也能保证图像不会被过度平滑。
空洞填充的函数代码如下，仅供参考：
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 }
