精华内容
下载资源
问答
  • 多视图三维重建
    2022-02-13 20:16:33

    MVS

    MVG

    更多相关内容
  • 研究数据驱动的多视图三维重建技术,包括如何从视频或图像数据中获取摄像机的内部参数,外部参数和场景三维几何信息,有效地解决大规模三维重建中的光照一致性和几何一致性等问题,为当今应用需求提供了一套完备的解决...
  • 基于多视图三维重建出有颜色的点云
  • #资源达人分享计划#
  • ReconstructionReconstruction是一个用于视角三维重建的项目,该项目参考的工作完成,使用10台佳能D1200相机、6盏闪光灯同步采集,最终由算法计算得到人脸三维模型。项目有VS2010和VS2015两个版本。项目依赖 静态...
  • 项目、论文地址:在公众号「计算机视觉工坊」,后台回复「pixelNeRF」,即可直接下载。概述作者提出了pixelNeRF,一个只需要输入单张或张图像,就能得到连续场景表示的学习框架。由...

    项目、论文地址:在公众号「计算机视觉工坊」,后台回复「pixelNeRF」,即可直接下载。

    概述

    作者提出了pixelNeRF,一个只需要输入单张或多张图像,就能得到连续场景表示的学习框架。由于现存的构建神经辐射场【1】的方法涉及到独立优化每个场景的表示,这需要许多校准的视图和大量的计算时间,因此作者引入了一种新的网络架构。实验结果表明,在所有情况下,pixelNeRF在新视图合成和单图像三维重建方面都优于当前最先进的工作。

    简介

    该项目主要研究的问题是如何从一个稀疏的输入视图集中合成这个场景的新视图,在可微神经渲染出现之前,这个长期存在的问题一直没有得到进展。同时,最近的神经渲染场NeRF通过编码体积密度和颜色,在特定场景的新视图合成方面表现出很好的效果。虽然NeRF可以渲染非常逼真的新视图,但它通常是不切实际的,因为它需要大量的位姿图像和冗长的场景优化。

    在这篇文章中,作者对上述方法进行了改进,与NeRF网络不使用任何图像特征不同的是,pixelNeRF将与每个像素对齐的空间图像特征作为输入。这种图像调节允许框架在一组多视图图像上进行训练,学习场景先验,然后从一个或几个输入图像中合成视图,如下图所示。

    PixelNeRF具有很多特点:首先,Pixel可以在多视图图像的数据集上面进行训练,而不需要任何额外的监督;其次,PixelNeRF预测输入图像的摄像机坐标系中的NeRF表示,而不是标准坐标系,这是泛化看不见的场景和物体类别的必要条件,因为在有多个物体的场景中,不存在明确的规范坐标系;第三,它是完全卷积的,这允许它保持图像和输出3D表示之间的空间对齐;最后,PixelNeRF可以在测试时合并任意数量的输入视图,且不需要任何优化。

    相关工作

    新视图合成:这是一个长期存在的问题,它需要从一组输入视图中构建一个场景的新视图。尽管现在有很多工作都已经取得了逼真的效果,但是存在比较多的问题,例如需要密集的视图和大量的优化时间。

    其他方法通过学习跨场景共享的先验知识,从单个或少数输入视图进行新的视图合成,但是这些方法使用2.5D表示,因此它们能够合成的摄像机运动范围是有限的。在这项工作中,作者提出了PixelNeRF,能够直接从相当大的基线中合成新视图。

    基于学习的三维重建:得益于深度学习的发展,单视图或多视图的三维重建也得到快速的发展。问题是,很多表示方法都需要3D模型进行监督,尽管多视图监督限制更小、更容易获取,其中的很多方法也需要物体的mask。相比之下,PixelNeRF可以单独从图像中训练,允许它应用到含有两个对象的场景而不需要修改。

    以观察者为中心的三维重建:对于3D学习任务,可以在以观察者为中心的坐标系(即视图空间)或以对象为中心的坐标系(即规范空间)中进行预测。大多数现存的方法都是在规范空间中进行预测,虽然这使得学习空间规律更加容易,但是会降低不可见对象和具有多个对象场景的预测性能。PixelNeRF在视图空间中操作,这在【2】中已经被证明可以更好地重建看不见的对象类别,并且不鼓励对训练集的记忆。下表是PixelNeRF和其他方法的对比:

    背景介绍:NeRF

    NeRF【1】将场景编码为颜色和密度的连续体积辐射场f。特别地,对于一个3D点x和

    观察方向单位向量d,f返回微分密度σ和RGB颜色c:f(x, d) = (σ, c)。体积辐射场可以通过下面的函数渲染成2D图像:

    其中T(t)处理遮挡。对于具有姿态P的目标视图,相机光线可以参数化为r(t)=o+td,o为光线原点(相机中心)。沿着相机光线在预定义的深度边界[tn,tf]之间计算积分。在实践中,这种积分是通过沿每个像素射线采样点的数值求积来近似的。

    然后,将摄影机光线r的渲染像素值与对应的真实像素值C(r)进行比较,最后的loss定义如下:

    其中R(P)是具有目标姿态的所有相机光线的集合。虽然NeRF实现了最新的视图合成,但它是一种基于优化的方法,每个场景必须单独优化,场景之间没有知识共享。这种方法不仅耗时,而且在单个或极稀疏视图的限制下,无法利用任何先验知识来加速重建或完成形状。

    基于图像的NeRF

    为了克服上面提到的关于NeRF的问题,作者提出了一种基于空间图像特征的NeRF结构。该模型由两个部分组成:一个完全卷积的图像编码器E(将输入图像编码为像素对齐的特征网格)和一个NeRF网络f(给定一个空间位置及其对应的编码特征,输出颜色和密度)。

    单视图pixelNeRF:首先固定坐标系为输入图像的视图空间,并在这个坐标系中指定位置和摄像机光线。给定场景的输入图像I,首先提取出它的特征量W=E(I)。然后,对于相机光线上的一个点x,通过使用已知的内参,将x投影到图像坐标π(x)上,然后在像素特征之间进行双线性插值来提取相应的图像特征向量W(π(x))。最后把图像特征连同位置和视图方向(都在输入视图坐标系统中)传递到NeRF网络:

    其中γ()是x上的位置编码。

    合并多个视图:多个视图提供了有关场景的附加信息,并解决了单视图固有的三维几何歧义。作者扩展了该模型,不同于现有的在测试时只使用单个输入视图的方法,它允许在测试时有任意数量的视图。

    在有多个输入视图的情况下,只假设相对的相机姿态是已知的,为了便于解释,可以为场景任意固定一个世界坐标系。把输入图像记为I,其相关联的摄像机记为P=[R t]。对于新的目标摄影机光线,将视图方向为d的点x转换到每个输入视图i的坐标系,转换如下:

    为了获得输出的密度和颜色,作者独立地处理每个视图坐标帧中的坐标和相应的特征,并在NeRF网络中聚合视图。将NeRF网络的初始层表示为f1,它分别处理每个输入视图空间中的输入,并将最终层表示为f2,它处理聚合视图。

    和单视图类似,作者将每个输入图像编码成特征体积W(i)=E(I(i))。对于点x(i),在投影图像坐标π(x(i))处从特征体W(i)中提取相应的图像特征,然后将这些输入传递到f1,以获得中间向量:

    最后用平均池化算子ψ将中间向量V(i)聚合并传递到最后一层f2,得到预测的密度和颜色:

    效果和对比

    特定类别的单视图重建

    特定类别的双视图重建

    特定类别的单视图和双视图重建结果对比

    参考文献:

    【1】Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik,Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng. Nerf: Representing scenes as neural radiance fields for view synthesis. In Eur. Conf. Comput. Vis., 2020

    【2】Daeyun Shin, Charless Fowlkes, and Derek Hoiem. Pixels, voxels, and views: A study of shape representations for single view 3d object shape prediction. In IEEE Conf. Comput.Vis. Pattern Recog., 2018.

    备注:作者系我们「3D视觉从入门到精通」知识星球特邀嘉宾:一个超干货的3D视觉学习社区

    本文仅做学术分享,如有侵权,请联系删文。

    下载1

    在「计算机视觉工坊」公众号后台回复:深度学习,即可下载深度学习算法、3D深度学习、深度学习框架、目标检测、GAN等相关内容近30本pdf书籍。

    下载2

    在「计算机视觉工坊」公众号后台回复:计算机视觉,即可下载计算机视觉相关17本pdf书籍,包含计算机视觉算法、Python视觉实战、Opencv3.0学习等。

    下载3

    在「计算机视觉工坊」公众号后台回复:SLAM,即可下载独家SLAM相关视频课程,包含视觉SLAM、激光SLAM精品课程。

    重磅!计算机视觉工坊-学习交流群已成立

    扫码添加小助手微信,可申请加入3D视觉工坊-学术论文写作与投稿 微信交流群,旨在交流顶会、顶刊、SCI、EI等写作与投稿事宜。

    同时也可申请加入我们的细分方向交流群,目前主要有ORB-SLAM系列源码学习、3D视觉CV&深度学习SLAM三维重建点云后处理自动驾驶、CV入门、三维测量、VR/AR、3D人脸识别、医疗影像、缺陷检测、行人重识别、目标跟踪、视觉产品落地、视觉竞赛、车牌识别、硬件选型、深度估计、学术交流、求职交流等微信群,请扫描下面微信号加群,备注:”研究方向+学校/公司+昵称“,例如:”3D视觉 + 上海交大 + 静静“。请按照格式备注,否则不予通过。添加成功后会根据研究方向邀请进去相关微信群。原创投稿也请联系。

    ▲长按加微信群或投稿

    ▲长按关注公众号

    3D视觉从入门到精通知识星球:针对3D视觉领域的视频课程(三维重建系列三维点云系列结构光系列手眼标定相机标定orb-slam3等视频课程)、知识点汇总、入门进阶学习路线、最新paper分享、疑问解答五个方面进行深耕,更有各类大厂的算法工程人员进行技术指导。与此同时,星球将联合知名企业发布3D视觉相关算法开发岗位以及项目对接信息,打造成集技术与就业为一体的铁杆粉丝聚集区,近2000星球成员为创造更好的AI世界共同进步,知识星球入口:

    学习3D视觉核心技术,扫描查看介绍,3天内无条件退款

     圈里有高质量教程资料、可答疑解惑、助你高效解决问题

    觉得有用,麻烦给个赞和在看~  

    展开全文
  • 基于多视图的几何的三维重建的具体研究方法,从方面考虑构成三维重建
  • 基于图像点特征的多视图三维重建

    千次阅读 2017-12-25 17:02:34
    基于图像点特征的多视图三维重建大致分成四步:一.图像特征点描述 ;二.图像特征点匹配; 三.基于图像的稀疏三维重建;四.视图稠密匹配与三维重建 一.图像特征点描述: 为描述所检测到的图像特征点,通常提取特征...

    基于图像点特征的多视图三维重建大致分成四步:一.图像特征点描述 ;二.图像特征点匹配;

    三.基于图像的稀疏三维重建;四.多视图稠密匹配与三维重建


    一.图像特征点描述:

    为描述所检测到的图像特征点,通常提取特征点周围的一块区域,并利用特征向量来描述该区域作为特征描述符。在某些特定拍摄条件下(如已矫正的立体像对),像素对应在各图像中的坐标仅在水平方向上有差异。在这种情况下,可使用搜索窗口内的像素差异平方和(Sum of squared differences,SSD),归一化互相关性(hormalized cross-correlation,NCC)等误差测度来计算特征点周围区域的相似度


    二.图像特征点匹配

    一种方法:特征匹配旨在确立特征之间的对应关系,最简单的是直接比较描述符之间的欧式距离,然后选取距离小于一定阈值的特征匹配作为图像对应。

            另一种方法:通过在特征空间中搜索最近邻来确定匹配关系。由于有的特征可能不存在对应,因此人需要设置阈值来剔除误匹配。在有充足训练数据的情况下,不同特征的阈值可通过学习获得;在训练数据不可得的情况下,通常采用距离比测试来剔除大部分误匹配。

    在高维特征空间中搜索候选匹配时,广泛使用的索引结构是多维搜索树。Muja和Lowe对这类方法进行了对比研究,并提出一种基于层次K均值树的优先搜索策略。在高维近似最近邻搜索问题中,多随机k-d树的性能最优。为了进一步剔除候选匹配中的误匹配,需要采用几何验证。

    具体来说,通常在随机采样一致(random sample consensus,RANSAC)算法中估计两视图间的基本矩阵,通过计算图像匹配与基本矩阵的拟合程度来剔除误匹配。对于各图像中的每个特征点,采用深度优先搜索找到它在其他图像中的对应,由此获得的多个视图间的对应称为特征轨迹。


    三.基于图像的稀疏三维重建

    1) 相机定标:一是利用定标物体上对各三维坐标点与其在图像平面上的投影坐标之间的约束关系来求解相机参数(可获得较高精度的相机参数估计);二是自定标。该方法不仅不需要定标物体,还不需要任何场景的先验知识,顾客满足更多实际应用的需求。绝对二次曲线(AC)及其对偶(DAC)成像的不变性,证明了两幅图像之间存在两个形如kruppa方程的二次非线性约束,利用三幅图像可实现内部参数不变的相机自定标。但,kruppa方程的解对噪声非常敏感,且解的个数随图像数量的增加以指数方式增长,限制了其在许多实际场合的应用。

    2) SFM:在相机内部参数已知的情况下,同时恢复相机的外部参数和三维场景结构的过程。为缓解SFM计算效率低下问题,第一种方案是利用现代图形处理单元GPU或者计算机集群的强大计算能力加速SFM中从图像匹配(SIFT)到集束优化(Bundle Adjustment)等各处理阶段。第二种方案是从数据组织和重建策略上改进SFM技术。


    四.多视图稠密匹配与三维重建

    第一类算法首先在三维体积中计算某种代价函数,然后从此三维体积中抽取物体的表面模型。此外,许多算法定义一个体积马尔可夫随机场(MRF),并采用最大流或者多路图割方法来提取最优的三维模型表面。

    第二类算法通过迭代的表面进化来最小化某种代价函数。这类方法中几何体的表现形式可为体素(voxel),水平集(level-set)和多边形网格(mesh)等形式。空间雕刻(space carving)及其变种也属于这一类算法,它们从初始体积重渐进剔除不一致的体素。

    第三类算法利用基于图像空间的方法计算一系列深度图(通常每幅输入图像对应一副深度图)。为确保所有深度图所表达的三维场景的一致性,该类方法在深度图上施加一致性约束,或者在后处理阶段融合多幅深度图信息创建完整的三维场景。

    第四类算法直接通过图像特征提取和匹配来获得图像对应,然后建立场景的稠密或者准稠密三维点云,最后从稠密三维点云中重建三维表面模型(目前总体性能最好)。


    展开全文
  • 使用OpenCV实现了图像的三维重建。 使用VS2015开发,程序运行后会读取images目录下的图片进行重建。 重建完成后,可以运行Viewer下的SfMViewer.exe查看重建结果。 详见博客 ...
  • 研究了基于多视图立体视觉的三维重建方法.给定一组待重建物体的图像,利用传统的立体视觉技术计算每幅图像的深度信息;然后设计一个基于体素划分的模型融合方法,将图像的深度信息融合为一个完整的三维初始模型;随后,对...
  • 三维重建源代码

    2016-10-19 11:55:57
    对于三维重建,有很大帮助,初学者可以学习
  • 多视图几何三维重建的基本原理: 从两个或者个视点观察同一景物,已获得在个不同的视角下对景物的张感知图像,运用三角测量的基本原理计算图像像素间位置偏差,获得景物的三维深度信息,这一个过程与人类...

    多视图几何三维重建的基本原理:

    从两个或者多个视点观察同一景物,已获得在多个不同的视角下对景物的多张感知图像,运用三角测量的基本原理计算图像像素间位置偏差,获得景物的三维深度信息,这一个过程与人类观察外面的世界的过程是一样的。

    SfM:

    SfM的全称为Structure from Motion,即通过相机的移动来确定目标的空间和几何关系,是三维重建的一种常见方法。它只需要普通的RGB摄像头,因此成本更低廉,且受环境约束较小,在室内和室外均能使用。但是,SfM背后需要复杂的理论和算法做支持,在精度和速度上都还有待提高,所以目前成熟的商业应用并不多。

    本文将实现一个简易的增量式SfM系统。

    书籍:计算机视觉中的多视图几何:原书第2版/(澳)理查德\cdot哈特利 北京:机械工业出版社,2019.8

    代码参考:OpenCV实现SfM(二):多目三维重建_arag2009的学习专栏-CSDN博客

    数据集来源:MVS Data Set – 2014 | DTU Robot Image Data Sets

    对极几何&基础矩阵

     

     

    特征点提取和匹配

    从上面的分析可知,要求取两个相机的相对关系,需要两幅图像中的对应点,这就变成的特征点的提取和匹配问题。

    对于图像差别较大的情况,推荐使用SIFT特征,因为SIFT对旋转、尺度、透视都有较好的鲁棒性。(如果差别不大,可以考虑其他更快速的特征,比如SURF、ORB等。 ),然后采用knn最近邻算法进行匹配,距离计算采用欧氏距离。

    由于OpenCV将SIFT包含在了扩展部分中,所以官网上下载的版本是没有SIFT的,

    为此需要到GitHub - opencv/opencv_contrib: Repository for OpenCV's extra modules下载扩展包opencv_contrib,并按照里面的说明重新编译OpenCV。

    编译opencv的时候要通过cmake -DOPENCV_EXTRA_MODULES_PATH将opencv_contrib的module编译进来,一定要保持两者的版本一致

    opencv链接:GitHub - opencv/opencv: Open Source Computer Vision Library

    重点说明的是,匹配结果往往有很多误匹配,为了排除这些错误,这里使用了Ratio Test方法,即使用KNN算法寻找与该特征最匹配的2个特征,若第一个特征的匹配距离与第二个特征的匹配距离之比小于某一阈值,就接受该匹配,否则视为误匹配。

    void match_features(
    	cv::Mat & query,
    	cv::Mat & train,
    	std::vector<cv::DMatch>& matches,
    	std::vector<std::vector<cv::DMatch>>& knn_matches
    )
    {
    	cv::BFMatcher matcher(cv::NORM_L2);
    	matcher.knnMatch(query, train, knn_matches, 2);
    
    	//获取满足Ratio Test的最小匹配的距离
    	float min_dist = FLT_MAX;
    	for (int r = 0; r < knn_matches.size(); ++r)
    	{
    		//Ratio Test
    		if (knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance)
    			continue;
    
    		float dist = knn_matches[r][0].distance;
    		if (dist < min_dist) min_dist = dist;
    	}
    	matches.clear();
    	for (size_t r = 0; r < knn_matches.size(); ++r)
    	{
    		//排除不满足Ratio Test的点和匹配距离过大的点
    		if (
    			knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance ||
    			knn_matches[r][0].distance > 5 * cv::max(min_dist, 10.0f)
    			)
    			continue;
    		//保存匹配点
    		matches.push_back(knn_matches[r][0]);
    	}
    
    }

    之后可以参考《计算机视觉中的多视图几何》第11章 基本矩阵估计 算法 11.4

    std::vector<std::vector<cv::DMatch>> knn_matches;
    match_features(this->leftimage->descriptor, this->rightimage->descriptor, m_matchs,knn_matches);//上文的匹配函数
    
    auto& kp1 = this->leftimage->keypoints;
    auto& kp2 = this->rightimage->keypoints;
    cv::Mat m_Fundamental;
    ransace_refine_matchs(kp1, kp2, m_matchs,m_Fundamental);//这里采用基于RANSAC的基础矩阵估计进行匹配点对的refine
    
    int last_num = 0;
    int next_num = 0;
    double T = 10;
    cv::Mat x1, x2;
    cv::Point2f p1, p2;
    double d;
    int index_kp2;
    do {
    
        last_num = m_matchs.size();
    
        std::vector<int> label1(kp1.size(), -1);
        std::vector<int> label2(kp2.size(), -1);
        for (int i = 0; i < m_matchs.size(); i++) {
            label1[m_matchs[i].queryIdx] = 1;
            label2[m_matchs[i].trainIdx] = 1;
        }
    
        for (int i = 0; i < knn_matches.size(); i++) {
            if (label1[i] < 0) {
    
                index_kp2 = knn_matches[i][0].trainIdx;
    
                if (label2[index_kp2] < 0) {
                    p1 = kp1[knn_matches[i][0].queryIdx].pt;
                    p2 = kp2[index_kp2].pt;
    
                    x1 = (cv::Mat_<double>(3, 1) << p1.x, p1.y, 1);
                    x2 = (cv::Mat_<double>(1, 3) << p2.x, p2.y, 1);
    
                    x1 = x2*m_Fundamental*x1;
                    d = abs(x1.at<double>(0, 0));
                    if (d < T) {
                        m_matchs.push_back(knn_matches[i][0]);
                        label1[i] = 1;
                        label2[index_kp2] = 1;
                    }
                }
    
            }
    
        }
        ransace_refine_matchs(kp1, kp2, m_matchs,m_Fundamental);//这里采用基于RANSAC的基础矩阵估计进行匹配点对的refine(再一次)
        next_num = m_matchs.size();
        std::cout <<"loop: "<< next_num - last_num << "\n";
    } while (next_num - last_num > 0);//如果匹配点对数目比上一次少则结束循环
    

    其中的基于RANSAC的基础矩阵估计,采用opencv的cv::findFundamentalMat,并设置参数method=cv::FM_RANSAC

    void ransace_refine_matchs(
    	std::vector<MyKeyPoint>& kp1, std::vector<MyKeyPoint>& kp2, 
    	std::vector<cv::DMatch>& matchs,
    	cv::Mat& m_Fundamental)
    {
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < m_matchs.size(); i++) {
    		_p1.push_back(kp1[m_matchs[i].queryIdx].pt);
    		_p2.push_back(kp2[m_matchs[i].trainIdx].pt);
    	}
    
    	std::vector<uchar> m_RANSACStatus;
    
    	m_Fundamental = cv::findFundamentalMat(_p1, _p2, m_RANSACStatus, cv::FM_RANSAC);//
    
    	std::vector<cv::DMatch> matchs_copy = matchs;
    	matchs.clear();
    
    	for (int i = 0; i < matchs_copy.size(); ++i) {
    		if (m_RANSACStatus[i] != 0) matchs.push_back(matchs_copy[i]);
    	}
    
    }

    双视图的三维重建

    bool find_transform(cv::Mat & K, std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, cv::Mat & R, cv::Mat & T, cv::Mat & mask)
    {
    	//根据内参矩阵获取相机的焦距和光心坐标(主点坐标)
    	double focal_length = 0.5*(K.at<double>(0) + K.at<double>(4));
    	cv::Point2d principle_point(K.at<double>(2), K.at<double>(5));
    
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < p1.size(); i++) {
    		_p1.push_back(p1[i].pt);
    		_p2.push_back(p2[i].pt);
    	}
    
    	//根据匹配点求取本征矩阵,使用RANSAC,进一步排除失配点
    	cv::Mat E = cv::findEssentialMat(_p1, _p2, focal_length, principle_point, cv::RANSAC, 0.999, 1.0, mask);
    	if (E.empty()) return false;
    
    	double feasible_count = cv::countNonZero(mask);
    	std::cout << (int)feasible_count << " -in- " << p1.size() << std::endl;
    	//对于RANSAC而言,outlier数量大于50%时,结果是不可靠的
    	if (feasible_count <= 15 || (feasible_count / p1.size()) < 0.6)
    		return false;
    
    	//分解本征矩阵,获取相对变换
    	int pass_count = cv::recoverPose(E, _p1, _p2, R, T, focal_length, principle_point, mask);
    
    	//同时位于两个相机前方的点的数量要足够大
    	if (((double)pass_count) / feasible_count < 0.7)
    		return false;
    	return true;
    }

     

    参考:https://blog.csdn.net/weixin_42587961/article/details/97241162#Ax0V_556

    这里我们用opencv的cv::triangulatePoints进行三角测量

    void reconstruct(cv::Mat & K, cv::Mat & R1, cv::Mat & T1, cv::Mat & R2, cv::Mat & T2,
    	std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, std::vector<MyPoint3f>& structure)
    {
    
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < p1.size(); i++) {
    		_p1.push_back(p1[i].pt);
    		_p2.push_back(p2[i].pt);
    	}
    
    	//两个相机的投影矩阵[R T],triangulatePoints只支持float型
    	cv::Mat proj1(3, 4, CV_32FC1);
    	cv::Mat proj2(3, 4, CV_32FC1);
    
    	R1.convertTo(proj1(cv::Range(0, 3), cv::Range(0, 3)), CV_32FC1);
    	T1.convertTo(proj1.col(3), CV_32FC1);
    
    	R2.convertTo(proj2(cv::Range(0, 3), cv::Range(0, 3)), CV_32FC1);
    	T2.convertTo(proj2.col(3), CV_32FC1);
    
    	cv::Mat fK;
    	K.convertTo(fK, CV_32FC1);
    	proj1 = fK * proj1;
    	proj2 = fK * proj2;
    	//三角重建
    	cv::Mat s;
    	cv::triangulatePoints(proj1, proj2, _p1, _p2, s);
    
    	structure.clear();
    	structure.reserve(s.cols);
    	for (int i = 0; i < s.cols; ++i)
    	{
    		cv::Mat_<float> col = s.col(i);
    		col /= col(3);  //齐次坐标,需要除以最后一个元素才是真正的坐标值
    		structure.push_back(MyPoint3f(cv::Point3f(col(0), col(1), col(2)), p1[i].color));
    	}
    }

    多视图的三维重建(增量式SfM)

    我们知道,两个相机之间的变换矩阵(旋转矩阵和平移向量)可以通过估计本质矩阵以及分解本质矩阵来实现。

    那么设相机1的坐标系为世界坐标系,现在加入第三幅视图,如何确定第三个相机(后面称为相机3)到相机1坐标系的变换矩阵呢?

    最简单的想法,就是沿用双目重建的方法,即在第三幅图像和第一幅图像之间提取特征点匹配然后三角测量重建。

    那么加入第四幅、第五幅,乃至更多呢?随着图像数量的增加,新加入的图像与第一幅图像的差异可能越来越大,特征点的提取变得异常困难,这时就不能再沿用双目重建的方法了。

    现在我们采用PnP的方法,该方法根据空间中的点与图像中的点的对应关系,求解相机在空间中的位置和姿态。也就是说,如果我知道相机1坐标系中的一些三维空间点,也知道这些空间点和相机3图像的2维点匹配,那么PnP就可以解出相机3和相机1之间的变换矩阵(旋转矩阵和平移向量)。

    多视图的三维重建的一个大致流程:

    1. 使用双目重建的方法,对头两幅图像进行重建,这样就得到了一些空间中的点,

    2. 加入第3幅图像后,使其与第2幅图像进行特征匹配,这些匹配点中,肯定有一部分也是图像1与图像2之间的匹配点,也就是说,这些匹配点中有一部分的空间坐标是已知的相机1坐标系下的空间点,同时又知道这些点在第3幅图像中的二维像素坐标

    3. 求解PnP所需的信息都有了。由于空间点的坐标都是相机1坐标系下的,所以由PnP求出的相机位置和姿态也是相机1坐标系下的,即相机3到相机1的变换矩阵。

    4. 利用相机3和相机1的变换矩阵继续三角测量重建更多的点云,并融合到已有的点云中

    后续加入更多的图像也是这个道理,并且要进行点云的融合。

    完整代码

    GlobalVTKUtils.h

    #pragma once
    #ifndef GLOBALVTKUTILS_H
    #define GLOBALVTKUTILS_H
    #include<vector>
    #include<array>
    class GlobalVTKUtils {
    public:
    	GlobalVTKUtils() {};
    	~GlobalVTKUtils() {};
    	static void ShowPoints(std::vector<std::array<float,3>> points,std::vector<std::array<float,3>> colors = std::vector<std::array<float,3>>());
    };
    #endif // !SHOWVTK_H

    GlobalVTKUtils.cpp

    #include"GlobalVTKUtils.h"
    #include <vtkActor.h>
    #include <vtkRenderer.h>
    #include <vtkRenderWindow.h>
    #include <vtkRenderWindowInteractor.h>
    #include <vtkProperty.h>
    #include <vtkInteractorStyleTrackball.h>
    #include <vtkInteractorStyleTrackballCamera.h>
    #include <vtkPoints.h>
    #include <vtkPolyData.h>
    #include <vtkCellArray.h>
    #include <vtkPolyDataMapper.h>
    #include <vtkSmartPointer.h>
    #include <vtkSphereSource.h>
    #include <vtkXMLPolyDataWriter.h>
    #include <vtkLookupTable.h>
    #include <vtkPointData.h>
    #include <vtkUnsignedCharArray.h>
    
    #include "vtkFloatArray.h"  //浮点型数组类
    #include "vtkAutoInit.h" 
    VTK_MODULE_INIT(vtkRenderingOpenGL2); // VTK was built with vtkRenderingOpenGL2
    VTK_MODULE_INIT(vtkInteractionStyle);
    VTK_MODULE_INIT(vtkRenderingFreeType);
    
    
    void GlobalVTKUtils::ShowPoints(std::vector<std::array<float, 3>> points, std::vector<std::array<float, 3>> colors)
    {
    	vtkSmartPointer<vtkRenderer> renderer = vtkSmartPointer< vtkRenderer>::New();
    	vtkSmartPointer<vtkRenderWindow> renderWindow = vtkSmartPointer<vtkRenderWindow>::New();
    	renderWindow->AddRenderer(renderer);
    	vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor = vtkSmartPointer<vtkRenderWindowInteractor>::New();
    	vtkSmartPointer<vtkInteractorStyleTrackballCamera> istyle = vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
    	renderWindowInteractor->SetRenderWindow(renderWindow);
    	vtkSmartPointer<vtkInteractorStyleTrackballCamera> style = vtkSmartPointer<vtkInteractorStyleTrackballCamera>::New();
    	renderWindowInteractor->SetInteractorStyle(style);//设置交互器风格 
    
    
    
    	std::vector<int> label(points.size(), -1);
    	double x1, y1, z1;
    	double x2, y2, z2;
    	double T = 0.1;
    	double value;
    	for (int i = 0; i < points.size(); i++) {
    
    		double minvalue = 1e10;
    		int min_index = -1;
    		for (int j = 0; j < points.size(); j++) {
    			if (j != i) {
    				x1 = points[i][0];y1 = points[i][1];z1 = points[i][2];
    				x2 = points[j][0];y2 = points[j][1];z2 = points[j][2];
    				value = std::sqrt(std::pow(x1 - x2, 2) + std::pow(y1 - y2, 2) + std::pow(z1 - z2, 2));
    				
    				if (value < minvalue) {
    					minvalue = value;
    					min_index = j;
    					if (minvalue < T)break;
    				}
    
    			}
    		}
    		//std::cout << i << " " << minvalue << std::endl;
    		if (min_index != -1 && minvalue>T) {
    			label[i] = 1;
    		}
    	}
    
    	std::vector<std::array<float, 3>> copy = points;
    	std::vector<std::array<float, 3>> copy_c = colors;
    	points.clear();
    	colors.clear();
    	for (int i = 0; i < copy.size(); i++) {
    		if (label[i] < 0) {
    			points.push_back(copy[i]);
    			colors.push_back(copy_c[i]);
    		}
    	}
    
    
    	//MAD去除离群点
    	double mean_pt[3];
    	for (int i = 0; i < points.size(); i++) {
    		auto pt = points[i];
    		for (int j = 0; j < 3; j++) {
    			mean_pt[j] += pt[j];
    		}
    	}
    	for (int j = 0; j < 3; j++) {
    		mean_pt[j] /= points.size();
    	}
    
    	double mad = 0;
    	for (int i = 0; i < points.size(); i++) {
    		auto pt = points[i];
    		double x, y, z;
    		x = pt[0];y = pt[1];z = pt[2];
    		mad += std::sqrt(std::pow(x - mean_pt[0], 2) + std::pow(y - mean_pt[1], 2) + std::pow(z - mean_pt[2], 2));
    	}
    	mad/=points.size();
    
    
    	vtkIdType idtype;
    	double x, y, z;
    	vtkPoints *vpoints = vtkPoints::New();
    	vtkCellArray *vcells = vtkCellArray::New();
    
    	vtkSmartPointer<vtkUnsignedCharArray> ptColor = vtkSmartPointer<vtkUnsignedCharArray>::New();
    	ptColor->SetNumberOfComponents(3);
    	bool has_color = true;
    	if (colors.empty())has_color = false;
    
    	for (int i = 0; i < points.size(); i ++) {
    		auto pt = points[i];
    		
    		x = pt[0];y = pt[1];z = pt[2];
    		//std::cout << pt[0] << " " << pt[1] << " " << pt[2] << "\n";
    		if (std::sqrt(std::pow(x - mean_pt[0], 2) + std::pow(y - mean_pt[1], 2) + std::pow(z - mean_pt[2], 2)) > 2*mad) {
    			continue;//不显示离群点
    		}
    		
    		idtype = vpoints->InsertNextPoint(x, y, z);
    		vcells->InsertNextCell(1, &idtype);
    		if (has_color) {
    			auto color = colors[i];
    			ptColor->InsertNextTuple3(color[2],color[1],color[0]);
    		}
    		
    		
    	}
    	// 渲染机制未知,需要同时设置点坐标与点坐标对应的verts
    	// verts中的id必须与点坐标对应
    	vtkPolyData *polyData = vtkPolyData::New();
    	polyData->SetPoints(vpoints);
    	polyData->SetVerts(vcells);
    	if (has_color) {
    		polyData->GetPointData()->SetScalars(ptColor);
    	}
    
    	//下面为正常的可视化流程,可设置的点云颜色、大小等已注释
    	vtkPolyDataMapper *mapper = vtkPolyDataMapper::New();
    	mapper->SetInputData(polyData);
    
    	vtkActor *actor = vtkActor::New();
    	actor->SetMapper(mapper);
    	actor->GetProperty()->SetPointSize(4);
    
    	renderer->AddActor(actor);
    
    	renderer->SetBackground(0, 0, 0);//设置背景颜色
    	renderWindow->Render();
    	renderWindow->SetSize(800, 800);//设置窗口大小
    	renderWindowInteractor->Start();
    
    	vtkSmartPointer<vtkXMLPolyDataWriter> writer =
            vtkSmartPointer<vtkXMLPolyDataWriter>::New();
        writer->SetFileName("polyData.vtp");
        writer->SetInputData(polyData);
        writer->Write();
    
    }
    

    MultiViewStereoReconstructor.h

    #pragma once
    #ifndef MULTIVIEWSTEREORECONSTRUCTOR
    #define MULTIVIEWSTEREORECONSTRUCTOR
    
    
    #include<vector>
    #include<opencv2/opencv.hpp>
    #include<array>
    #include<memory>
    
    class MyKeyPoint :public cv::KeyPoint {
    public:
    	MyKeyPoint() {};
    	MyKeyPoint(cv::Point2f _pt, std::vector<float> _color = std::vector<float>(), 
    		float _size = 1, float _angle = -1, float _response = 0, int _octave = 0, int _class_id = -1)
    		: cv::KeyPoint(_pt,_size,_angle,_response,_octave,_class_id),color(_color){};
    	MyKeyPoint(cv::KeyPoint _keypoint) :cv::KeyPoint(_keypoint) {};
    	MyKeyPoint(cv::KeyPoint _keypoint,cv::Vec3b _color) :cv::KeyPoint(_keypoint) {
    		setColor(_color);
    	};
    	MyKeyPoint(cv::KeyPoint _keypoint,std::vector<float> _color) :cv::KeyPoint(_keypoint),color(_color) {};
    	void setColor(cv::Vec3b _color);
    	~MyKeyPoint() {};
    
    	bool operator==(MyKeyPoint &that);
    
    public:
    	std::vector<float> color;
    
    };
    
    class MyPoint3f :public cv::Point3f {
    public:
    	MyPoint3f() {};
    	MyPoint3f(cv::Point3f _p) :cv::Point3f(_p) {};
    	MyPoint3f(cv::Point3f _p,cv::Vec3b _color) :cv::Point3f(_p) {
    		setColor(_color);
    	};
    	MyPoint3f(cv::Point3f _p,std::vector<float> _color) :cv::Point3f(_p),color(_color) {};
    	void setColor(cv::Vec3b _color);
    public:
    	std::vector<float> color;
    };
    
    class UnitView {
    public:
    	UnitView(cv::Mat& _image) {
    		image = _image.clone();
    	};
    	~UnitView() {};
    
    	bool extract_feature(cv::Ptr<cv::Feature2D> _executor);
    public:
    	cv::Mat image;
    	cv::Mat descriptor;
    	std::vector<MyKeyPoint> keypoints;
    	cv::Mat rotation;//旋转矩阵
    	cv::Mat motion;
    };
    
    class PointCloudModel {
    public:
    	PointCloudModel() {};
    	~PointCloudModel() {};
    public:
    	std::vector<std::vector<int>> correspond_struct_idx;
    	std::vector<MyPoint3f> structure;
    	cv::Mat K;
    	cv::Mat distCoeffs;
    };
    
    class MyMatch {
    
    public:
    	MyMatch(int _type=1): match_type(_type){};
    	~MyMatch() {};
    
    	bool SetUp();
    	bool Build();
    	bool init_structure();
    	void match_features(cv::Mat& query, cv::Mat& train, 
    		std::vector<cv::DMatch>& matches,
    		std::vector<std::vector<cv::DMatch>>& knn_matches);
    	void ransace_refine_matchs(std::vector<MyKeyPoint>& kp1, std::vector <MyKeyPoint>& kp2, 
    		std::vector<cv::DMatch>& matchs,cv::Mat& m_Fundamental);
    	bool find_transform(cv::Mat& K, std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, cv::Mat& R, cv::Mat& T, cv::Mat& mask);
    	void maskout_points(std::vector<MyKeyPoint>& p1, cv::Mat& mask);
    	void reconstruct(
    		cv::Mat& K, 
    		cv::Mat& R1, cv::Mat& T1, 
    		cv::Mat& R2, cv::Mat& T2, 
    		std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, 
    		std::vector<MyPoint3f>& structure);
    	void get_objpoints_and_imgpoints(
    		std::vector<cv::DMatch>& matches,
    		std::vector<int>& last_struct_indices,
    		std::vector<MyPoint3f>& global_structure,
    		std::vector<MyKeyPoint>& key_points,
    		std::vector<cv::Point3f>& object_points,
    		std::vector<cv::Point2f>& image_points);
    	void fusion_structure(
    		std::vector<cv::DMatch>& matches,
    		std::vector<int>& last_struct_indices,
    		std::vector<int>& next_struct_indices, 
    		std::vector<MyPoint3f>& global_structure,
    		std::vector<MyPoint3f>& local_structure);
    	cv::Mat get_disparity(cv::Mat& img1,cv::Mat& img2);
    
    
    public:
    	
    	std::shared_ptr<UnitView> leftimage;
    	std::shared_ptr<UnitView> rightimage;
    	std::shared_ptr<PointCloudModel> global_model;
    	int match_type = 1;//0表示这个匹配是第一张图像和第二张图像的匹配,否则为1; 默认为1;
    
    private:
    	std::vector<MyPoint3f> m_structure;
    	std::vector<cv::DMatch> m_matchs;
    	std::vector<MyKeyPoint> m1;
    	std::vector<MyKeyPoint> m2;
    };
    
    class MultiViewStereoReconstructor {
    
    public:
    	MultiViewStereoReconstructor() {};
    
    	~MultiViewStereoReconstructor() {};
    
    	void SetK(cv::Mat& _K) { K = _K.clone(); };
    	void SetDistCoeffs(cv::Mat _distCoeffs) { distCoeffs = _distCoeffs.clone(); };
    
    	void SetData(std::vector<cv::Mat>& _images);
    
    	bool Execute();
    
    	void Visualize();
    
    
    
    public:
    	std::shared_ptr<PointCloudModel> global_model;
    	std::vector<std::shared_ptr<UnitView>> images;
    	std::vector<MyMatch> matchs;
    	cv::Mat K;
    	cv::Mat distCoeffs;
    };
    
    
    #endif // !MULTIVIEWSTEREORECONSTRUCTOR
    
    

    MultiViewStereoReconstructor.cpp

    #pragma once
    #include "MultiViewStereoReconstructor.h"
    #include <opencv2/xfeatures2d/nonfree.hpp>
    #include "GlobalVTKUtils.h"
    #include <omp.h>
    #include <thread>
    
    void MyKeyPoint::setColor(cv::Vec3b _color)
    {
    	this->color.resize(3);
    	for (int i = 0; i < 3; i++)this->color[i] = _color[i];
    }
    
    bool MyKeyPoint::operator==(MyKeyPoint & that)
    {
    	return this->pt == that.pt;
    }
    
    void MyPoint3f::setColor(cv::Vec3b _color)
    {
    	this->color.resize(3);
    	for (int i = 0; i < 3; i++)this->color[i] = _color[i];
    }
    
    
    bool UnitView::extract_feature(cv::Ptr<cv::Feature2D> _executor)
    {
    
    	std::vector<cv::KeyPoint> _keypoints;
    	_executor->detectAndCompute(this->image, cv::noArray(), _keypoints, this->descriptor);
    
    	if (_keypoints.size() <= 10) return false;
    
    	for (auto& key : _keypoints) {
    		this->keypoints.push_back(MyKeyPoint(key, this->image.at<cv::Vec3b>(key.pt.y, key.pt.x)));
    	}
    
    }
    
    
    
    bool MyMatch::SetUp()
    {
    	//std::cout << "MyMatch::SetUp\n";
    	std::vector<std::vector<cv::DMatch>> knn_matches;
    	match_features(this->leftimage->descriptor, this->rightimage->descriptor, m_matchs,knn_matches);
    
    	auto& kp1 = this->leftimage->keypoints;
    	auto& kp2 = this->rightimage->keypoints;
    	cv::Mat m_Fundamental;
    	ransace_refine_matchs(kp1, kp2, m_matchs,m_Fundamental);
    	
    
    #if 1
    	int last_num = 0;
    	int next_num = 0;
    	double T = 10;
    	cv::Mat x1, x2;
    	cv::Point2f p1, p2;
    	double d;
    	int index_kp2;
    	do {
    
    		last_num = m_matchs.size();
    
    		std::vector<int> label1(kp1.size(), -1);
    		std::vector<int> label2(kp2.size(), -1);
    		for (int i = 0; i < m_matchs.size(); i++) {
    			label1[m_matchs[i].queryIdx] = 1;
    			label2[m_matchs[i].trainIdx] = 1;
    		}
    
    		for (int i = 0; i < knn_matches.size(); i++) {
    			if (label1[i] < 0) {
    
    				index_kp2 = knn_matches[i][0].trainIdx;
    
    				if (label2[index_kp2] < 0) {
    					p1 = kp1[knn_matches[i][0].queryIdx].pt;
    					p2 = kp2[index_kp2].pt;
    
    					x1 = (cv::Mat_<double>(3, 1) << p1.x, p1.y, 1);
    					x2 = (cv::Mat_<double>(1, 3) << p2.x, p2.y, 1);
    
    					x1 = x2*m_Fundamental*x1;
    					d = abs(x1.at<double>(0, 0));
    					if (d < T) {
    						m_matchs.push_back(knn_matches[i][0]);
    						label1[i] = 1;
    						label2[index_kp2] = 1;
    					}
    				}
    				
    			}
    			
    		}
    		ransace_refine_matchs(kp1, kp2, m_matchs,m_Fundamental);
    		next_num = m_matchs.size();
    		std::cout <<"loop: "<< next_num - last_num << "\n";
    	} while (next_num - last_num > 0);
    
    #endif
    
    	for (int i = 0; i < m_matchs.size(); i++) {
    
    		m1.push_back(kp1[m_matchs[i].queryIdx]);
    		m2.push_back(kp2[m_matchs[i].trainIdx]);
    	}
    
    
    
    	//std::cout << "leftimage->image.size(): " << this->leftimage->image.size() << "\n";
    	//std::cout << "leftimage->descriptor.size(): " << this->leftimage->descriptor.size() << "\n";
    	//std::cout << "rightimage->descriptor.size(): " << this->rightimage->descriptor.size() << "\n";
    	//std::cout << "leftimage->keypoints.size(): " << this->leftimage->keypoints.size() << "\n";
    	//std::cout << "rightimage->keypoints.size(): " << this->rightimage->keypoints.size() << "\n";
    	//std::cout << "m1.size(): " << m1.size() << "\n";
    	//std::cout << "m2.size(): " << m2.size() << "\n";
    
    
    #if 0
    	int gap = 10;
    	std::vector<MyKeyPoint> &key1 = kp1;
    	std::vector<MyKeyPoint> &key2 = kp2;
    
    	std::vector<cv::KeyPoint> tkey1;
    	std::vector<cv::KeyPoint> tkey2;
    	std::vector<cv::DMatch> tMatches;
    	int InlinerCount = 0;
    	for (int i = 0; i < m_matchs.size(); i += gap) {
    		tkey1.push_back(key1[m_matchs[i].queryIdx]);
    		tkey2.push_back(key2[m_matchs[i].trainIdx]);
    		cv::DMatch match;
    		match.queryIdx = InlinerCount;
    		match.trainIdx = InlinerCount;
    		InlinerCount++;
    		tMatches.push_back(match);
    	}
    
    	cv::Mat img_matches;
    	cv::drawMatches(
    		this->leftimage->image, tkey1, this->rightimage->image, tkey2,
    		tMatches, img_matches, cv::Scalar::all(-1), cv::Scalar::all(-1),
    		std::vector<char>(), cv::DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS);
    	cv::imwrite("FASTResult.jpg", img_matches);
    	cv::imshow("Match0", img_matches);
    
    	cv::waitKey(0);
    #endif
    
    
    
    
    
    	//std::cout << "done MyMatch::SetUp\n";
    	return true;
    }
    
    void MyMatch::match_features(
    	cv::Mat & query,
    	cv::Mat & train,
    	std::vector<cv::DMatch>& matches,
    	std::vector<std::vector<cv::DMatch>>& knn_matches
    )
    {
    	//std::vector<std::vector<cv::DMatch>> knn_matches;
    	cv::BFMatcher matcher(cv::NORM_L2);
    	matcher.knnMatch(query, train, knn_matches, 2);
    
    	//std::cout << knn_matches.size() << "\n";
    	//std::cout <<"knn_matches: "<< knn_matches.size() << "\n";
    	//for (int i = 0; i < knn_matches.size(); i++) {
    	//	std::cout << i << " " << knn_matches[i][0].queryIdx << " " << knn_matches[i][0].trainIdx<<" "<<knn_matches[i][0].distance << "\n";
    	//	std::cout << i << " " << knn_matches[i][1].queryIdx << " " << knn_matches[i][1].trainIdx<<" "<<knn_matches[i][1].distance << "\n";
    	//	std::cout << "\n\n";
    	//}
    	//获取满足Ratio Test的最小匹配的距离
    	float min_dist = FLT_MAX;
    	for (int r = 0; r < knn_matches.size(); ++r)
    	{
    		//Ratio Test
    		if (knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance)
    			continue;
    
    		float dist = knn_matches[r][0].distance;
    		if (dist < min_dist) min_dist = dist;
    	}
    	matches.clear();
    	for (size_t r = 0; r < knn_matches.size(); ++r)
    	{
    		//排除不满足Ratio Test的点和匹配距离过大的点
    		if (
    			knn_matches[r][0].distance > 0.6*knn_matches[r][1].distance ||
    			knn_matches[r][0].distance > 5 * cv::max(min_dist, 10.0f)
    			)
    			continue;
    		//保存匹配点
    		matches.push_back(knn_matches[r][0]);
    	}
    
    }
    void MyMatch::ransace_refine_matchs(
    	std::vector<MyKeyPoint>& kp1, std::vector<MyKeyPoint>& kp2, 
    	std::vector<cv::DMatch>& matchs,
    	cv::Mat& m_Fundamental)
    {
    	//auto& kp1 = this->leftimage->keypoints;
    	//auto& kp2 = this->rightimage->keypoints;
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < m_matchs.size(); i++) {
    		_p1.push_back(kp1[m_matchs[i].queryIdx].pt);
    		_p2.push_back(kp2[m_matchs[i].trainIdx].pt);
    	}
    	//cv::Mat m_Fundamental;
    	std::vector<uchar> m_RANSACStatus;
    
    	m_Fundamental = cv::findFundamentalMat(_p1, _p2, m_RANSACStatus, cv::FM_RANSAC);//
    
    	std::vector<cv::DMatch> matchs_copy = matchs;
    	matchs.clear();
    
    	for (int i = 0; i < matchs_copy.size(); ++i) {
    		if (m_RANSACStatus[i] != 0) matchs.push_back(matchs_copy[i]);
    	}
    
    }
    
    
    bool MyMatch::find_transform(cv::Mat & K, std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, cv::Mat & R, cv::Mat & T, cv::Mat & mask)
    {
    	//根据内参矩阵获取相机的焦距和光心坐标(主点坐标)
    	double focal_length = 0.5*(K.at<double>(0) + K.at<double>(4));
    	cv::Point2d principle_point(K.at<double>(2), K.at<double>(5));
    
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < p1.size(); i++) {
    		_p1.push_back(p1[i].pt);
    		_p2.push_back(p2[i].pt);
    	}
    
    	//根据匹配点求取本征矩阵,使用RANSAC,进一步排除失配点
    	cv::Mat E = cv::findEssentialMat(_p1, _p2, focal_length, principle_point, cv::RANSAC, 0.999, 1.0, mask);
    	if (E.empty()) return false;
    
    	double feasible_count = cv::countNonZero(mask);
    	std::cout << (int)feasible_count << " -in- " << p1.size() << std::endl;
    	//对于RANSAC而言,outlier数量大于50%时,结果是不可靠的
    	if (feasible_count <= 15 || (feasible_count / p1.size()) < 0.6)
    		return false;
    
    	//分解本征矩阵,获取相对变换
    	int pass_count = cv::recoverPose(E, _p1, _p2, R, T, focal_length, principle_point, mask);
    
    	//同时位于两个相机前方的点的数量要足够大
    	if (((double)pass_count) / feasible_count < 0.7)
    		return false;
    	return true;
    }
    
    void MyMatch::maskout_points(std::vector<MyKeyPoint>& p1, cv::Mat & mask)
    {
    	std::vector<MyKeyPoint> p1_copy = p1;
    	p1.clear();
    
    	for (int i = 0; i < mask.rows; ++i) {
    		if (mask.at<uchar>(i) > 0) p1.push_back(p1_copy[i]);
    	}
    }
    
    void MyMatch::reconstruct(cv::Mat & K, cv::Mat & R1, cv::Mat & T1, cv::Mat & R2, cv::Mat & T2,
    	std::vector<MyKeyPoint>& p1, std::vector<MyKeyPoint>& p2, std::vector<MyPoint3f>& structure)
    {
    
    	std::vector<cv::Point2f> _p1, _p2;
    	for (int i = 0; i < p1.size(); i++) {
    		_p1.push_back(p1[i].pt);
    		_p2.push_back(p2[i].pt);
    	}
    
    
    	//两个相机的投影矩阵[R T],triangulatePoints只支持float型
    	cv::Mat proj1(3, 4, CV_32FC1);
    	cv::Mat proj2(3, 4, CV_32FC1);
    
    	R1.convertTo(proj1(cv::Range(0, 3), cv::Range(0, 3)), CV_32FC1);
    	T1.convertTo(proj1.col(3), CV_32FC1);
    
    	R2.convertTo(proj2(cv::Range(0, 3), cv::Range(0, 3)), CV_32FC1);
    	T2.convertTo(proj2.col(3), CV_32FC1);
    
    	cv::Mat fK;
    	K.convertTo(fK, CV_32FC1);
    	proj1 = fK * proj1;
    	proj2 = fK * proj2;
    	//三角重建
    	cv::Mat s;
    	cv::triangulatePoints(proj1, proj2, _p1, _p2, s);
    
    	structure.clear();
    	structure.reserve(s.cols);
    	for (int i = 0; i < s.cols; ++i)
    	{
    		cv::Mat_<float> col = s.col(i);
    		col /= col(3);  //齐次坐标,需要除以最后一个元素才是真正的坐标值
    		structure.push_back(MyPoint3f(cv::Point3f(col(0), col(1), col(2)), p1[i].color));
    	}
    }
    
    bool MyMatch::init_structure()
    {
    	cv::Mat & K = global_model->K;
    
    	std::cout << "MyMatch::init_structure\n";
    	cv::Mat R, T;//旋转矩阵和平移向量
    	cv::Mat mask;//mask中大于零的点代表匹配点,等于零代表失配点
    
    	find_transform(K, m1, m2, R, T, mask);
    	maskout_points(m1, mask);//去除野值
    	maskout_points(m2, mask);//去除野值
    
    	cv::Mat R0 = cv::Mat::eye(3, 3, CV_64FC1);
    	cv::Mat T0 = cv::Mat::zeros(3, 1, CV_64FC1);
    	reconstruct(K, R0, T0, R, T, m1, m2, m_structure);
    
    	this->leftimage->rotation = R0;
    	this->leftimage->motion = T0;
    	this->rightimage->rotation = R;
    	this->rightimage->motion = T;
    
    
    	auto& correspond_struct_idx = global_model->correspond_struct_idx;
    	correspond_struct_idx.clear();
    	correspond_struct_idx.resize(2);
    	correspond_struct_idx[0].resize(this->leftimage->keypoints.size(), -1);
    	correspond_struct_idx[1].resize(this->rightimage->keypoints.size(), -1);
    
    	//填写头两幅图像的结构索引
    	int idx = 0;
    	for (int i = 0; i < m_matchs.size(); ++i)
    	{
    		if (mask.at<uchar>(i) == 0)
    			continue;
    		correspond_struct_idx[0][m_matchs[i].queryIdx] = idx;
    		correspond_struct_idx[1][m_matchs[i].trainIdx] = idx;
    		++idx;
    	}
    	//m_correspond_struct_idx = correspond_struct_idx[1];
    
    	//初始化,直接赋值
    	global_model->structure = m_structure;
    
    	return true;
    }
    
    
    void MyMatch::get_objpoints_and_imgpoints(
    	std::vector<cv::DMatch>& matches,
    	std::vector<int>& last_struct_indices,
    	std::vector<MyPoint3f>& global_structure,
    	std::vector<MyKeyPoint>& key_points,
    	std::vector<cv::Point3f>& object_points,
    	std::vector<cv::Point2f>& image_points)
    {
    	object_points.clear();
    	image_points.clear();
    
    	for (int i = 0; i < matches.size(); ++i)
    	{
    		int query_idx = matches[i].queryIdx;
    		int train_idx = matches[i].trainIdx;
    
    		int struct_idx = last_struct_indices[query_idx];
    		if (struct_idx < 0) continue;
    
    		object_points.push_back(global_structure[struct_idx]);
    		image_points.push_back(key_points[train_idx].pt);
    	}
    }
    
    void MyMatch::fusion_structure(
    	std::vector<cv::DMatch>& matches,
    	std::vector<int>& last_struct_indices,
    	std::vector<int>& next_struct_indices,
    	std::vector<MyPoint3f>& global_structure,
    	std::vector<MyPoint3f>& local_structure)
    {
    
    	for (int i = 0; i < matches.size(); ++i)
    	{
    		int query_idx = matches[i].queryIdx;
    		int train_idx = matches[i].trainIdx;
    
    		int struct_idx = last_struct_indices[query_idx];
    		if (struct_idx >= 0) //若该点在空间中已经存在,则这对匹配点对应的空间点应该是同一个,索引要相同
    		{
    			next_struct_indices[train_idx] = struct_idx;
    			continue;
    		}
    		//若该点在空间中已经存在,将该点加入到结构中,且这对匹配点的空间点索引都为新加入的点的索引
    		global_structure.push_back(local_structure[i]);
    		last_struct_indices[query_idx] = next_struct_indices[train_idx] = global_structure.size() - 1;
    	}
    
    
    
    }
    
    cv::Mat MyMatch::get_disparity(cv::Mat & img1, cv::Mat & img2)
    {
    	// Compute disparity
    	cv::Mat disparity;
    	cv::Ptr<cv::StereoMatcher> pStereo = cv::StereoSGBM::create(0, 16, 5);
    	
    	pStereo->compute(img1, img2, disparity);
    
    	return disparity;
    }
    
    
    
    bool MyMatch::Build()
    {
    	cv::Mat& K = global_model->K;
    	//std::cout << "MyMatch::Build\n";
    	if (this->match_type == 0) {
    		init_structure();
    	}
    	else {
    		//增量方式重建剩余的图像
    		cv::Mat r, R, T;
    		std::vector<cv::Point3f> object_points;
    		std::vector<cv::Point2f> image_points;
    
    		auto& global_structure = global_model->structure;
    		auto& global_struct_idx = global_model->correspond_struct_idx;
    		auto& last_struct_idx = global_struct_idx[global_struct_idx.size() - 1];
    		//获取第i幅图像中匹配点对应的三维点,以及在第i+1幅图像中对应的像素点
    		get_objpoints_and_imgpoints(
    			m_matchs,
    			last_struct_idx,
    			global_structure,
    			this->rightimage->keypoints,
    			object_points,
    			image_points
    		);
    
    		//求解变换矩阵
    		cv::solvePnPRansac(object_points, image_points, K, cv::noArray(), r, T);
    		//将旋转向量转换为旋转矩阵
    		cv::Rodrigues(r, R);
    
    		this->rightimage->rotation = R;
    		this->rightimage->motion = T;
    
    		reconstruct(
    			K, this->leftimage->rotation, this->leftimage->motion, R, T,
    			m1, m2, m_structure);
    
    
    		std::vector<int> next_struct_indices(this->rightimage->keypoints.size(), -1);
    
    		//将新的重建结果与之前的融合
    		fusion_structure(
    			m_matchs,
    			last_struct_idx,
    			next_struct_indices,
    			global_structure,
    			m_structure);
    
    		global_struct_idx.push_back(next_struct_indices);
    	}
    
    
    
    	return true;
    }
    
    
    
    
    void MultiViewStereoReconstructor::SetData(std::vector<cv::Mat>& _images) {
    	for (auto& im : _images) {
    
    		std::shared_ptr<UnitView> im_ptr = std::shared_ptr<UnitView>(new UnitView(im));
    		//UnitView* im_ptr = new UnitView(im);
    		images.push_back(im_ptr);
    	}
    	cv::Ptr<cv::Feature2D> _executor = cv::xfeatures2d::SIFT::create(0, 3, 0.04, 10);
    
    #pragma omp parallel for num_threads(2*omp_get_num_procs()-1)
    	for (int i = 0; i < images.size(); i++) {
    		printf("thread: %d\n", std::this_thread::get_id());
    		images[i]->extract_feature(_executor);
    	}
    
    	matchs.resize(images.size() - 1);
    }
    
    
    
    
    
    
    bool MultiViewStereoReconstructor::Execute() {
    	if (K.empty()) {
    		std::cout << "must set K matrix!\n";
    		return false;
    	}
    
    	global_model = std::shared_ptr<PointCloudModel>(new PointCloudModel());
    	global_model->K = K.clone();
    	global_model->distCoeffs = distCoeffs.clone();
    
    	//#pragma omp parallel for num_threads(2*omp_get_num_procs()-1)
    	for (int i = 0; i < images.size() - 1; i++) {
    		MyMatch &match = matchs[i];
    		if (i == 0) {
    			match.match_type = 0;
    		}
    		match.leftimage = images[i];
    		match.rightimage = images[i + 1];
    		match.global_model = global_model;
    
    		//std::cout << "Matching images " << i << " - " << i + 1 << std::endl;
    		std::printf("thread: %d Matching images %d - %d\n", std::this_thread::get_id(), i, i + 1);
    		match.SetUp();
    
    	}
    
    	for (int i = 0; i < images.size() - 1; i++) {
    		matchs[i].Build();
    	}
    
    };
    
    
    
    void MultiViewStereoReconstructor::Visualize()
    {
    	auto& global_structure = global_model->structure;
    	std::vector < std::array<float, 3>> tpoints;
    	std::vector < std::array<float, 3>> tcolors;
    	for (size_t i = 0; i < global_structure.size(); ++i)
    	{
    		std::array<float, 3> tp;
    		tp[0] = global_structure[i].x;
    		tp[1] = global_structure[i].y;
    		tp[2] = global_structure[i].z;
    		tpoints.push_back(tp);
    
    		std::array<float, 3> tc;
    		tc[0] = global_structure[i].color[0];
    		tc[1] = global_structure[i].color[1];
    		tc[2] = global_structure[i].color[2];
    		tcolors.push_back(tc);
    	}
    	std::cout <<"points.size(): "<< tpoints.size() << "\n";
    	GlobalVTKUtils::ShowPoints(tpoints, tcolors);
    	//GlobalVTKUtils::ShowPoints(tpoints);
    }
    

    main.cpp

    
    #include"MultiViewStereoReconstructor.h"
    #include <io.h>
    #include <chrono>
    #ifndef INIT_TIMER2
    #define INIT_TIMER2(name)     \
      std::chrono::high_resolution_clock::time_point timerStart##name = std::chrono::high_resolution_clock::now();
    
    #endif // !INIT_TIMER2
    
    #ifndef STOP_TIMER2
    #define STOP_TIMER2(name,title)                                                                                               \
      std::cout << "RUNTIME of " << title << ": "                                                                           \
                << std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::high_resolution_clock::now() - timerStart##name).count()\
                << " ms " << std::endl;
    #endif // !STOP_TIMER2
    
    void get_files(const std::string & path, std::vector<std::string>& files)
    {
    	std::cout << "[jhq] get_files" << std::endl << std::endl;
    	//文件句柄
    	long long hFile = 0;
    	//文件信息,_finddata_t需要io.h头文件
    	struct _finddata_t fileinfo;
    	std::string p;
    	if ((hFile = _findfirst(p.assign(path).append("/*").c_str(), &fileinfo)) != -1) {
    		do {
    			if (strcmp(fileinfo.name, ".") != 0 && strcmp(fileinfo.name, "..") != 0) {
    				files.push_back(p.assign(path).append(fileinfo.name));
    			}
    		} while (_findnext(hFile, &fileinfo) == 0);
    		_findclose(hFile);
    	}
    }
    
    
    int main()
    {
    	INIT_TIMER2(MultiViewStereoReconstructor)
    
    	std::vector<std::string> img_names;
    	get_files("./Viewer16/", img_names);
    	
    
    
    
    	cv::Mat K(cv::Matx33d(
    		2889.84, 0, 825.91,
    		0, 2881.08, 602.73,
    		0, 0, 1));
    	cv::Mat distCoeffs = (cv::Mat_<float>(1, 5) << -0.06122164316121147, -1.243003308575242, 0.001462736730339558, -0.0001684641589496723, 1.749846418870851);
    
    	std::vector<cv::Mat> images;
    
    	for (auto it = img_names.begin(); it != img_names.end(); ++it){
    		std::cout << *it << std::endl;
    #if 1
    		cv::Mat im = cv::imread(*it);
    		cv::Size image_size = im.size();
    		cv::Mat mapx = cv::Mat(image_size, CV_32FC1);
    		cv::Mat mapy = cv::Mat(image_size, CV_32FC1);
    		cv::Mat R =cv::Mat::eye(3, 3, CV_32F);
    
    		cv::initUndistortRectifyMap(K, distCoeffs, R, K,image_size, CV_32FC1, mapx, mapy);//用来计算畸变映射
    		cv::Mat newimage;
    		cv::remap(im, newimage, mapx, mapy, cv::INTER_LINEAR);//把求得的映射应用到图像上
    
    		images.push_back(newimage);
    #else
    		images.push_back(cv::imread(*it));
    #endif
    	}
    
    
    	MultiViewStereoReconstructor mvs;
    	mvs.SetK(K);
    	mvs.SetDistCoeffs(distCoeffs);
    	mvs.SetData(images);
    	mvs.Execute();
    
    	STOP_TIMER2(MultiViewStereoReconstructor,"MultiViewStereoReconstructor")
    
    
    	mvs.Visualize();
    
    
    	
    	return 0;
    }

    展开全文
  • 多视图几何中的三维重建

    千次阅读 2020-09-21 11:05:27
    1)B站多视图几何三维重建视频讲解:https://www.bilibili.com/video/BV1Sj411f73e 2)武汉大学摄影测量与遥感专业博士李迎松的CSDN: https://blog.csdn.net/rs_lys/article/list/1 涉及的内容主要是 sfm, Patch...
  • 在完成两视图三维重建之后,接下来就是进行视图重建。视图重建的难点在于如何确定后一个相机到世界坐标系的变换矩阵。
  • 目视觉三维重建研究综述

    千次阅读 热门讨论 2020-11-20 21:22:15
    目视觉三维重建研究综述 目视觉三维重建研究综述 三维重建在视觉方面具有很高的研究价值,在机器人视觉导航、智能车环境感知系统以及虚拟现实中被广泛应用. 1.三维重建方法 从整体上来看,三维重建技术主要通过...
  • 上期文章介绍了用于三维重建的深度学习框架MVSNet[1],这也是如今比较主流的深度估计的神经网络框架。框架的原理按照双目立体匹配框架步骤:匹配代价构造、匹配代价累积、深度估计和深度图优化四个步骤。使用过...
  • 为了方便大家了解基于多视图立体的三维重建技术,更重要的是能亲手利用开源数据集或者自己采集的影像跑一遍流程,进而对整个流程更为熟悉,本文整理了近年来几种经典的基于传统方法和基于深度学习方法的三维重建...
  • 【实例简介】《计算机视觉中的多视图几何》匹配源代码,matlab中关于三维重建的源代码【实例截图】【核心代码】matlab中关于三维重建的源代码,《计算机视觉中的多视图几何》匹配源代码├── vgg_examples│ ├── ...
  • 针对图像序列的高效、稳定三维重建问题,提出一种基于独立三视图三维重建方法。首先以三视图为重建单元,叠加三视图约束,实现基础矩阵的鲁棒性计算,然后通过三角测量恢复得到三视图对应的三维场景局部三维点云。...
  • ‍点击上方“3D视觉工坊”,选择“星标”干货第一时间送达作者丨Longway来源丨三维重建学习笔记本人从2018年开始接触基于深度学习的单视图三维重建算法,中间因为考研等原因休息了一年,...
  • 基于传统多视图几何的三维重建算法 基于深度学习的三维重建算法 总地来说,尽管目前传统的三维重建算法依旧占据研究的主要部分,但是越来越的研究者开始关注于用CNN探索三维重建,或者说...
  • CT头颅三维重建

    2018-04-19 18:12:17
    在MATLAB R2016a版本下运行,里面包含头颅的CT图像,运行结果可以显示出头颅的立体图像以及三视图
  • 第二种基于深度图融合的方法在前面已经介绍过:基于图像的三维重建——深度图计算方法1-SGM/tSGM(9)以及基于图像的三维重建——深度图计算方法2-PatchMatch(10),这两种都是通过计算得到的不同视角的深度图进行融合...
  • 三维重建开源项目汇总

    千次阅读 2021-07-10 00:55:48
    作者丨次方AIRX@知乎来源丨https://zhuanlan.zhihu.com/p/141946874编辑丨计算机视觉工坊1、Meshroom ⭐4,474Meshroom是一款基于...
  • 1. 基于图像的图像3D重建 传统上首先使用 ... 然后,此输出用作Multi-View Stereo(多视图立体)的输入,以恢复场景冲密集表示。 path/to/project/sparse ---------------- 包含所有重建组件的稀疏模型 path/to/proj...
  • 基于多视图三维重建中特征提取与特征匹配并行化研究_冯兵 基于多视图三维重建中特征提取与特征匹配并行化研究_冯兵 基于多视图三维重建中特征提取与特征匹配并行化研究_冯兵
  • 基于深度学习的视觉三维重建研究总结 1.三维重建意义 三维重建作为环境感知的关键技术之一,可用于自动驾驶、虚拟现实、运动目标监测、行为分析、安防监控和重点人群监护等。现在每个人都在研究识别,但识别只是...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 8,160
精华内容 3,264
关键字:

多视图三维重建