真实场景的双目立体匹配（Stereo Matching）获取深度图详解
点击上方“小白学视觉”，选择加"星标"或“置顶”
重磅干货，第一时间送达双目立体匹配一直是双目视觉的研究热点，双目相机拍摄同一场景的左、右两幅视点图像，运用立体匹配匹配算法获取视差图，进而获取深度图。而深度图的应用范围非常广泛，由于其能够记录场景中物体距离摄像机的距离，可以用以测量、三维重建、以及虚拟视点的合成等。
之前有两篇博客简要讲过OpenCV3.4中的两种立体匹配算法效果比较。以及利用视差图合成新视点。里面用到的匹配图像对是OpenCV自带校正好的图像对。而目前大多数立体匹配算法使用的都是标准测试平台提供的标准图像对，比如著名的有如下两个：
MiddleBury
KITTI
但是对于想自己尝试拍摄双目图片进行立体匹配获取深度图，进行三维重建等操作的童鞋来讲，要做的工作是比使用校正好的标准测试图像对要多的。因此博主觉得有必要从用双目相机拍摄图像开始，捋一捋这整个流程。
主要分四个部分讲解：
摄像机标定（包括内参和外参）
双目图像的校正（包括畸变校正和立体校正）
立体匹配算法获取视差图，以及深度图
利用视差图，或者深度图进行虚拟视点的合成
注：如果没有双目相机，可以使用单个相机平行移动拍摄，外参可以通过摄像机自标定算出。我用自己的手机拍摄，拍摄移动时尽量保证平行移动。
一、摄像机标定
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_8UC1K ----内参矩阵，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);}
}
注：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 }
来源：博客园 作者：一度逍遥
链接：https://www.cnblogs.com/riddick/p/8486223.html
好消息！
小白学视觉知识星球
开始面向外开放啦
真实场景的双目立体匹配（Stereo Matching）获取深度图详解相关推荐
- 一文详解如何在真实场景的双目立体匹配（Stereo Matching）获取深度图
来源:博客园 作者:一度逍遥 双目立体匹配一直是双目视觉的研究热点,双目相机拍摄同一场景的左.右两幅视点图像,运用立体匹配匹配算法获取视差图,进而获取深度图.而深度图的应用范围非常广泛,由于其能够记录 ...
- php 获取流文件大小,php获取文件大小详解
本文主要和大家分享php获取文件大小详解,希望能帮助到大家. 通过filesize函数可以取得文件的大小,文件大小是以字节数表示的.$filename = '/data/webroot/usercod ...
- 荒岛求生html5母狼攻,荒岛求生各资源作用及获取方法详解
荒岛求生各资源作用及获取方法详解 2018-03-06 13:45:22来源:游戏下载编辑:野狐禅评论(0) 荒岛求生地图中有很多资源可以获取,它们各有什么作用呢?下面就为大家带来荒岛求生资源作用及获 ...
- 猛兽之地服务器维护,猛兽之地Roguelands全材料获取途径详解
以下就是小编为大家带来的猛兽之地Roguelands全材料获取途径详解: 蘑菇镇主要产出: 发光蘑菇 屎壳螂 少量怪兽抓 古老废墟主要产出: 甲壳碎片 辛辣种子 星空岩 折磨大陆主要产出: 星之果 少 ...
- 妖怪手帐获取服务器信息失败,妖怪手帐妖怪全获取途径详解分享
妖怪手帐妖怪怎么获得?游戏中小伙伴们需要操控妖怪来完成副本等操作,妖怪是小伙伴们的伙伴.那么如何获得妖怪呢?有哪些获取的方法途径?下面是小编带来的攻略解析,一起来关注下哦! 主要妖怪来源途径: 分为玉 ...
- 魅蓝e android版本,魅蓝E怎么刷机 root权限获取教程详解
魅蓝E怎么刷机 root权限获取教程详解 来源:www.18183.com作者:雾里看海时间:2016-08-15 魅蓝E怎么刷机?这边有几个步骤,大家别遗漏了哦,下面小编具体给大家讲讲魅蓝E怎么刷机 ...
- 超级跑跑服务器维修,梦幻西游手游神兽超级跑跑如何获取 超级泡泡获取方法详解...
梦幻西游手游最新价增加了新的神兽超级泡泡,有着超强的实力与超萌外表,成为许多玩家的追求,那么超级泡泡该如何获取呢,有什么方法呢,下面一起来看看吧. 一,兑换法 简单来说想要获得超级泡泡最直接快速的方法 ...
- 并发场景下的幂等问题——分布式锁详解
简介:本文从钉钉实人认证场景的一例数据重复问题出发,分析了其原因是因为并发导致幂等失效,引出幂等的概念.针对并发场景下的幂等问题,提出了一种实现幂等可行的方法论,结合通讯录加人业务场景对数据库幂等问题 ...
- java 查看文件属性_java File类获取文件属性详解
你知道java File类获取文件属性方法吗?下面的文章要给大家讲解的就是这个方面的内容,希望下面的内容可以对你有所帮助哦. 在Java中获取文件属性信息的第一步是先创建一个File类对象并指向一个已 ...
最新文章
- mysql 二进制日志
- 汇编:模拟C语言实现break与continue
- Android StageFrightMediaScanner源码解析
- ThinkPHP5路由图解
- 深入理解Delete（JavaScript）
- Canvas-图片旋转
- css clearfix_如何使用CSS清除浮点数（clearfix）？
- 买铅笔（洛谷-P1909）
- jquery：TypeError: $(...).on is not a function
- 设计灵感|没有素材时，如何用排版取胜？
- Android学习总结（3）——Handler深入详解
- 腾讯社招 —— 应用宝后端工程师-电话面试
- live2d内嵌html,live2d web端加载moc3模型
- 数学建模 ---斯皮尔曼相关系数
- STM32F4 之STM32CubeMx编程学习
- java 音频 网络传输_如何流式传输音频？
- Typora高亮颜色设置
- 夏季 肝病患者养生保健重点在哪 请详解
- Selenium 页面加载慢(一直转圈)
- 完全平方数（C语言，调用函数）