双目立体匹配流程详解
原文链接 https://www.cnblogs.com/riddick/p/8486223.html
真实场景的双目立体匹配(Stereo Matching)获取深度图详解
KITTI:http://www.cvlibs.net/datasets/kitti/eval_scene_flow.php?benchmark=stereo。
但是对于想自己尝试拍摄双目图片进行立体匹配获取深度图,进行三维重建等操作的童鞋来讲,要做的工作是比使用校正好的标准测试图像对要多的。因此博主觉得有必要从用双目相机拍摄图像开始,捋一捋这整个流程。
注:如果没有双目相机,可以使用单个相机平行移动拍摄,外参可以通过摄像机自标定算出。我用自己的手机拍摄,拍摄移动时尽量保证平行移动。
摄像机的内参包括,fx, fy, cx, cy,以及畸变系数[k1,k2,p1,p2,k3],详细就不赘述。我用手机对着电脑拍摄各个角度的棋盘格图像,棋盘格图像如图所示:
使用OpenCV3.4+VS2015对手机进行内参标定。标定结果如下,手机镜头不是鱼眼镜头,因此使用普通相机模型标定即可:
图像分辨率为:3968 x 2976。上面标定结果顺序依次为fx, fy, cx, cy, k1, k2, p1, p2, k3, 保存到文件中供后续使用。
比如:我拍摄这样两幅图像,以后用来进行立体匹配和虚拟视点合成的实验。
① 利用摄像机内参进行畸变校正,手机的畸变程度都很小,校正后的两幅图如下:
② 将上面两幅畸变校正后的图作为输入,使用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();
畸变校正前面已经介绍过,利用畸变系数进行畸变校正即可,下面说一下立体校正。
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);
对极线校正结果如下所示,查看对极线校正结果是否准确,可以通过观察若干对应点是否在同一行上粗略估计得出:
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);
左视差图(不可靠视差值为0) 右视差图(不可靠视差值为 (numberOfDisparities+1)*16 )
视差图中视差值不可靠的视差大多数是由于遮挡引起,或者光照不均匀引起。既然牛逼如SGBM也觉得不可靠,那与其留着做个空洞,倒不如用附近可靠的视差值填充一下。
空洞填充也有很多方法,在这里我检测出空洞区域,然后用附近可靠视差值的均值进行填充。填充后的视差图如下:
视差的单位是像素(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米的深度,足够了。
注:视差图和深度图中均有计算不正确的点,此文意在介绍整个流程,不特别注重算法的优化,如有大神望不吝赐教。
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 }
双目立体匹配流程详解相关推荐
- 双目立体匹配步骤详解
点击上方"3D视觉工坊",选择"星标" 干货第一时间送达 来源:https://blog.csdn.net/rs_lys/article/details/833 ...
- python 立体匹配_双目立体匹配步骤详解
点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 来源: https://blog.csdn.net/rs_lys/ ...
- matlab相关系数影像匹配_双目立体匹配步骤详解
作者:李迎松
- [双目视差] 立体匹配步骤详解
文章目录 立体匹配步骤详解 Step1 匹配代价计算 Step2 代价聚合 Step3 视差计算 Step4 视差优化 立体匹配步骤详解 Step1 匹配代价计算 匹配代价计算的目的是衡量待匹配像素与 ...
- 跨境电商三单对碰三单申报流程详解
跨境电商三单对碰三单申报流程详解 概要:三单申报是指"电子订单.电子运单.支付凭证". 1.电子订单: 适合申报企业类型"电商企业.电商交易平台.电商境内代理企业&quo ...
- Android事件流程详解
Android事件流程详解 网络上有不少博客讲述了android的事件分发机制和处理流程机制,但是看过千遍,总还是觉得有些迷迷糊糊,因此特地抽出一天事件来亲测下,向像我一样的广大入门程序员详细讲述an ...
- 基于spark mllib_Spark高级分析指南 | 机器学习和分析流程详解(下)
- 点击上方"中国统计网"订阅我吧!- 我们在Spark高级分析指南 | 机器学习和分析流程详解(上)快速介绍了一下不同的高级分析应用和用力,从推荐到回归.但这只是实际高级分析过程 ...
- View的绘制-draw流程详解
目录 作用 根据 measure 测量出的宽高,layout 布局的位置,渲染整个 View 树,将界面呈现出来. 具体分析 以下源码基于版本27 DecorView 的draw 流程 在<Vi ...
- 杂志订阅管理系统c++_电池管理系统BMS功能安全开发流程详解
点击上面 "电动知家"可以订阅哦! BMS功能安全开发流程详解 BMS和ISO26262 - BMS & ISO26262简介 BMS即Battery Management ...
最新文章
- 异常详细信息: System.Web.HttpException: 请求在此上下文中不可用
- 全球及中国管道运输行业建设发展与投资战略规划报告2022版
- MM模块部分名词解释
- c++大作业迷宫游戏 规定时间内完成_小学生做作业磨蹭的7个原因及对策!太准了~...
- 编译3.0的linux内核,Ubuntu 编译 Linux 3.0-rc4 内核
- Oracle非重要文件恢复,redo、暂时文件、索引文件、password文件
- Linux系统下I/O操作讲解,深入了解实战高级I/O编程
- Zune支持哪些格式?
- 反爬虫策略分析及处理
- 阿里 离线数据同步工具 DataX 初试
- VS2005 安装WTL
- 动画 | 什么是红黑树?(与2-3树等价)
- 天龙八部荣耀版体验服服务器未响应,《天龙八部荣耀版》创新竖版手游官网-合区来了!体验服合区测试解析...
- 火山PC自绘高级表格及超级列表框
- 玲珑杯Unity开发心得——进度条界面(异步加载游戏场景)
- 《C语言程序设计》江宝钏主编-习题5-3-动态最大值!!!
- 【算法详解-数学】(1)φ的基本知识
- cannot create temp dir for unpacking extensions
- ansys 定义变量(关联模型)
- pac文件提取服务器,[工具使用] privoxy 实现 PAC 请求过滤