精华内容
下载资源
问答
  • 多视图几何三维重建
    2021-07-19 09:20:00

    上期文章介绍了用于三维重建的深度学习框架MVSNet[1],这也是如今比较主流的深度估计的神经网络框架。框架的原理按照双目立体匹配框架步骤:匹配代价构造、匹配代价累积、深度估计和深度图优化四个步骤。使用过MVSNet的同学会发现,MVSNet使用3D的卷积神经网络对聚合后的代价体进行正则化,防止在学习过程中,受到低概率的错误匹配影响。

    但使用三维卷积神经网络(U-Net[2]),会造成非常大的GPU消耗,使得我们在使用过程中,受到一定的限制。同时,因为该正则化的模块,导致普通GPU单卡下无法训练和测试较高分辨率的影像集,也会影响深度估计范围和估计精度。

    图1 MVSNet代价体正则化

    针对该问题,本篇文章将介绍CVPR2019的R-MVSNet[3],并简单根据代码,介绍运行步骤和对应的问题。

    1、R-MVSNet

    R-MVSNet同样是香港科技大学姚遥等人在CVPR2019上提出的一种深度学习框架,它在MVSNet的基础上,解决了正则化过程中GPU消耗大、无法估计

    更多相关内容
  • 多视图几何三维重建的基本原理: 从两个或者个视点观察同一景物,已获得在个不同的视角下对景物的张感知图像,运用三角测量的基本原理计算图像像素间位置偏差,获得景物的三维深度信息,这一个过程与人类...

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

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

    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;
    }

    展开全文
  • MVS的深度学习实战系列到目前为止就先告一段落,本系列我们首先回顾了MVS的传统原理和传统方法,并介绍了COLMAP等软件的具体操作,之后,介绍了MVSNet、R-MVSNet和Cascade-MVSNet三个深度学习框架和三维重建的...

    点击上方“计算机视觉工坊”,选择“星标”

    干货第一时间送达

    MVSNet在2018年提出后,在估计深度图的应用中取得了非常好的结果。应用CNN于立体匹配的技术也使得传统的匹配效率整体提高。但是因为使用3D卷积神经网络进行深度正则化处理,所以即便在比较低的分辨率(900*600)下,也需要比较高的GPU消耗。针对该问题,该团队在CVPR2019上提出利用循环神经网络对3D代价体进行切片处理,大幅度减少GPU消耗,使得该网络框架不仅可以估计更大范围的场景,且估计精度更高。

    本篇文章仍将就MVSNet内存消耗大的问题,介绍CVPR2020的一篇文章:Cascade Cost Volume for High-Resolution Multi-View Stereo and Stereo Matching. 该文章沿用MVSNet深度估计的框架,具体创新在于改进Cost Volume的构造方式,使得利用深度学习估计深度时,在较低GPU消耗上估计高分辨率、大场景的深度。

    1、背景介绍

    基于深度学习的多视图立体,例如经典的MVSNet网络架构,通常会构造一个三维的代价体去回归场景的深度值,但MVSNet常受限于显存限制而无法对高分辨率的影像进行深度估计。

    在MVSNet框架的基础上,多种方法对显存增长问题提出了改进方案,上一篇文章我们介绍了R-MVSNet,该方法利用循环神经网络GRU,对三维代价体进行切片,这样不仅保留了靠前的深度和纹理信息,也减少了GPU的消耗,深度估计精度和深度估计范围要优于MVSNet,不同方法的比较结果可通过图1体现。

    图1 效果比较图

    不同于R-MVSNet,本篇文章(以下简称“Cascade-MVSNet”)则仍旧使用了MVSNet网络框架中的3D卷积神经网络对深度代价体进行正则化,但不同于其代价体的构造,Cascade-MVSNet利用链式代价体构造的策略,先估计较为粗糙的深度值,然后再进一步缩小深度估计范围,提高深度估计精度,实现了在较小的GPU消耗的条件下,得到较高分辨率和较高精度的深度图,经过稠密重建后,Cascade-MVSNet的结果也比之前所提到的方法要更为完整(图2)。

    图2 稠密重建结果比较

    2、代价体构造回顾

    沿用立体匹配的视差估计框架,通常深度学习方法都会构造一个沿深度方向的代价体,例如MVSNet借鉴平面扫描算法的原理,构造由沿不同深度而前视平行的平面组成一个相机椎体,然后每一个相机椎体经过采样,变成长宽一致的特征体,通过可微分化的单应性变换将不同视角下的特征扭曲到参考视角上,构成代价体。

    一般来说,代价体的长宽则是由输入影像分辨率提供的,深度范围则是通过稀疏重建后的结果提供先验条件。这样的构造方式,使得代价体正则化时,三维代价体的内存消耗会以三次指数的速率增长(图3)。

    图3 代价体构造示意图

    上一篇文章中R-MVSNet使用循环神经网络,对将代价体沿深度切成不同的深度图,利用GRU结构进行正则化处理,相比较MVSNet,能减少GPU的消耗。可是循环神经网络会涉及一个遗忘的过程,导致网络不能很好地保留像素周围的纹理信息,所以点云完整度不能得到很好地保留(图4)。

    图4 RNN代价体正则化过程

    3、Cascade-MVSNet

    为了解决信息保留和GPU消耗的两个问题,Cascade-MVSNet提出一种级联的代价体构造方法,并输出从粗到细的深度估计值。

    首先来看其整个网络的架构(图5)。

    图5 Cascade网络框架示意图

    可以看到,整个网络的结构还是沿袭了MVSNet框架(图6),还是以多视图的影像作为输入,然后经过可微分单应性变换形成不同视角下的特征体,在通过代价体构造,形成一个代价体,之后通过回归估计深度值。这里不再对MVSNet作更多的赘述,感兴趣的朋友可以回顾之前的文章:MVSNet。

    图6 MVSNet深度估计框架

    此篇文章不同的是:整个网络利用级联式代价体的构造策略。首先,原始输入影像,利用特征金字塔网络先对原始影像进行降采样,降低特征体的分辨率,使得可以拥有较为精确的深度估计范围。通过初始的MVSNet框架估计出低分辨率下的深度图后,进入下一阶段。

    下一阶段开始时,先进行上采样,然后以上一层的深度估计范围作为参考,确定改成的深度估计范围和深度估计间隔。最后输出一个较高分辨率的深度图。

    如图7,根据稀疏重建给予的先验深度范围,第一阶段的深度估计范围将包含整个场景。所以,会和低分辨率影像建立一个体量较小的代价体(图5)。在之后的阶段中,深度估计范围会进一步缩小。

    图7 不同阶段的深度估计范围变化

    除了深度估计范围的缩小,深度估计间隔也会进一步缩小,越小的深度估计间隔代表着越精细的深度估计精度。

    在第k个阶段,假设当前的深度估计范围Rk和深度平面估计间隔Ik,所以对应的深度估计平面数可以计算。当该阶段的图像分辨率固定后,一个更大的深度估计范围会产生更为精确的深度估计结果,但同时也会提升GPU的消耗。同时根据特征金字塔网络的特点,在级联式代价体构造过程中,每个阶段将按照上一阶段的两倍数量,即在每个阶段的图像的分辨率分别是之前的两倍。

    4、Loss的设置

    Cascade-MVSNet仍旧是监督学习下的网络结构,同时,因为仍旧沿用MVSNet的网络架构,所以Loss的构造和MVSNet的loss构造形式类似,不同的是,Cascade-MVSNet使用的是级联式的学习策略,所以Loss构造定义为:

    其中,Lk指的是第k个阶段的总Loss,λk则表示当前阶段的权重。一般来说,分辨率越高,设置的权重越大。

    5、Cascade-MVSNet实战操作

    首先,再次感谢Yaoyao(香港科技大学)给出已经预处理好的数据,感谢Alibaba集团提供的开源代码。因为Cascade-MVSNet的深度估计框架沿用MVSNet,所以其输入和MVSNet要求一致,这里不再重复说明,详细数据处理内容,请大家回顾实战系列-MVSNet。

    1)环境配置

    参考Alibaba的github主页中的installation,即可完成环境配置。

    (https://github.com/alibaba/cascade-stereo)

    2)深度估计

    环境配置结束后,需要打开CasMVSNet文件夹(图8),整个cascade-stereo包含多视图立体和双目立体,本篇文章只关注多视图立体方面的应用。所以进入CasMVNSet的文件夹。

    图8 CasMVSNet目录

    具体的配置这里不再重复,Github上的README.md介绍已经十分详尽。运行之后可以得到的结果如下:

    ·Scan10 数据库

                     原图                             b) 深度图                       c) 深度置信度图

    图9 深度图估计结果

    图10 稠密重建结果

    3)具体分析

    代码分析:

    Github上提供了一个shell的脚本,在终端直接运行即可,以下对代码进行简要介绍:

    #!/usr/bin/env bash
    TESTPATH="data/DTU/dtu_test_all"
    TESTLIST="lists/dtu/test.txt"
    CKPT_FILE=$1
    python test.py --dataset=general_eval --batch_size=1 --testpath=$TESTPATH --testlist=$TESTLIST --loadckpt $CKPT_FILE ${@:2}
    

    上述bash中,需要首先设置Test图片的位置,这里设置为下载好的数据集即可,TestList表示设置需要对那几个数据进行深度处理和重建,例如,我想对10,15,24号数据集重建,那Testlist中只要包含scan10, scan15, scan24即可。

    CKPT_FILE是训练好的网络的地址,通过官网可以下载预训练的模型。

    在test.py中,主要有几个flag需要注意,代码如下:

    parser.add_argument('--max_h', type=int, default=864, help='testing max h')
    parser.add_argument('--max_w', type=int, default=1152, help='testing max w')
    parser.add_argument('--filter_method', type=str, default='normal', choices=["gipuma", "normal"], help="filter method")
    parser.add_argument('--fusibile_exe_path', type=str, default='../fusibile/fusibile')
    parser.add_argument('--prob_threshold', type=float, default='0.9')
    parser.add_argument('--num_consistent', type=float, default='4')
    

    注意,可以改变测试分辨率的大小,但需要为32的倍数,且最好保证为16:9的长宽比例。其次,使用gipuma进行稠密重建时,需要预下载fusion这个库,这个在https://github.com/YoYo000/fusibile可以下载,然后在终端输入:

    Mkdir build
    Cd build
    Cmake ..
    Make
    

    之后将 --fusibile_exe_path的地址设置为可执行文件所在的地方。第三个要注意的flag是—prob_threshold,这个表示gipuma算法中利用置信度过滤深度图的阈值,如图11)所示,如果当前像素的置信度较高(c-上部),则说明当前像素的深度估计比较准确,若置信度较低,则说明当前像素上,深度估计有多个峰值,不能确定哪一个是最优深度,也表示当前像素估计的深度不准。因此,利用该阈值对深度图进行过滤,在稠密重建的时候,置信度较低的像素就不进行深度融合。

    图11 置信度解释示意图[2]

    结果分析:

    通过输出多个中间结果,对Cascade结果进行分析,首先是跑分结果,通过表1可得,利用Cascade的代价图构造方式在精度和完整度上都比之前的方法更优,比传统方法提升了近15个点,比R-MVSNet方法提升了近10个点。

    表1 跑分比较结果表

    其次,从深度图估计结果看(图12),深度图较为平滑、完整,且估计精度比较高,估计错误的大多数是背景或者白色无纹理地区。通过对当前深度图的点云映射(图13),当前网络在原图视角下可以保持较高的分辨率,当在meshlab中旋转视角,可以看到当前视角估计出来的稠密点云效果完整,主体部分基本正确,背景部分存在偏差是正常的体现,在整体融合的时候可以滤除掉。

    图12 scan15深度估计结果

    图13 单图恢复点云结果

    更多结果:

    图14 深度估计结果

    图15 稠密重建结果

    以上结果均由笔者亲自试验,效果能达到论文中描述的水平,且GPU消耗基本在7GB以内,可行性比较高。

    6、总结

    Cascade-MVSNet是一种利用级联式代价体构造,解决MVSNet在深度估计过程中的GPU消耗过大问题、点云完整度不高的深度学习方法。通过实验,该方法在现有的MVSNet框架下,能以较小的GPU消耗,得到较高精度的深度估计结果,同时也能保留较高的完整度。

    MVS的深度学习实战系列到目前为止就先告一段落,本系列我们首先回顾了MVS的传统原理和传统方法,并介绍了COLMAP等软件的具体操作,之后,介绍了MVSNet、R-MVSNet和Cascade-MVSNet三个深度学习框架和三维重建的pipeline,希望大家通过本系列的分享,加深对MVS三维重建的理解。

    附录

    开源数据集下载链接:在公众号计算机视觉工坊」,后台回复mvsnet,即可直接下载。

    我们需要下载的是用于测试的数据集,所以打开百度云链接时,点击 “mvsnet”,再点击“preprocessed_inputs”,下载其中的“dtu.zip”和“tankandtemples.zip”即可(图16)。

    图16 开源数据集下载

    感谢香港科技大学姚遥开源的数据集链接。

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

    下载1

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

    下载2

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

    下载3

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

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

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

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

    ▲长按加微信群或投稿

    ▲长按关注公众号

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

    展开全文
  • 点击上方“3D视觉工坊”,选择“星标”干货第一时间送达1. 概述MVS是一种从具有一定重叠度的多视图视角中恢复场景的稠密结构的技术,传统方法利用几何、光学一致性构造匹配代价,进行匹配代价...

    点击上方“3D视觉工坊”,选择“星标”

    干货第一时间送达

    1. 概述

    MVS是一种从具有一定重叠度的多视图视角中恢复场景的稠密结构的技术,传统方法利用几何、光学一致性构造匹配代价,进行匹配代价累积,再估计深度值。虽然传统方法有较高的深度估计精度,但由于存在在缺少纹理或者光照条件剧烈变化的场景中的错误匹配,传统方法的深度估计完整度还有很大的提升空间。

    近年来卷积神经网络已经成功被应用在特征匹配上,提升了立体匹配的精度。在这样的背景下,香港科技大学Yaoyao等人,在2018年提出了一种基于深度学习的端到端深度估计框架——MVSNet[1]。

    虽然这两年已经有更多的框架在精度和完整度上面超过了MVSNet,但笔者发现,无论是基于监督学习的R-MVSNet 、Cascade – MVSNet等网络,还是基于自监督学习的M3VSNet,核心网型设计都是在借鉴MVSNet而完成的,而且MVSNet也是比较早期且较为完整的三维重建深度学习框架,了解该框架的原理、数据IO与实际操作能加深对2020年以来各种新方法的理解。

    本篇的主要内容为MVSNet的深度估计框架介绍、深度估计原理及如何利用MVSNet进行三维重建的具体操作。

    (注意:本方法需要11GB以上的GPU,推荐使用至少NVIDIA 1080Ti,且同样建议使用Linux进行操作。驱动建议安装CUDA9.0以上版本。)

    2. MVSNet框架

    图1 MVSNet网络形状

    由上图所示,MVSNet是一种监督学习的方法,以一个参考影像和多张原始影像为输入,而得到参考影像深度图的一种端到端的深度学习框架。网络首先提取图像的深度特征,然后通过可微分投影变换构造3D的代价体,再通过正则化输出一个3D的概率体,再通过soft argMin层,沿深度方向求取深度期望,获得参考影像的深度图。

    3. 基于深度学习的深度估计原理

    按照立体匹配框架[2],基于M VSNet的深度图估计步骤如下:深度特征提取,构造匹配代价,代价累计,深度估计,深度图优化。

    深度特征提取。深度特征指通过神经网络提取的影像特征,相比传统SIFT、SURF的特征有更好的匹配精度和效率[3]。经过视角选择之后,输入已经配对的N张影像,即参考影像和候选集。首先利用一个八层的二维卷积神经网络(图2)提取立体像对的深度特征Fi,输出32通道的特征图.

    图2 特征采样的神经网络

    为防止输入的像片被降采样后语义信息的丢失,像素的临近像素之间的语义信息已经被编码到这个32通道的特征中,并且各个图像提取过程的网络是权值共享的。

    考虑到对亚像素的深度估计,以保证深度图平滑,该单应性矩阵是完全可以微分的。通过投影变换,N张影像可形成N个特征体(图3-b),这个特征体就是匹配代价的表示。

    代价累积。MVSNet的代价累积通过构造代价体实现的。代价体是一个由长、宽与参考影像长宽一样的代价图在深度方向连接而成的三维结构(图4-a),在深度维度,每一个单位表示一个深度值。其中,某一深度的代价图上面的像素表示参考影像同样的像素在相同深度处,与候选集影像的匹配代价。

    图4  代价累计结果

    当已知概率体时,最简单的方法可以获取参考影像的所有像素在不同深度的概率图,按照赢者通吃原则[6]直接估计深度图。

    然而,赢者通吃原则无法在亚像素级别估计深度,造成深度突变、不平滑情况。所以需要沿着概率体的深度方向,以深度期望值作为该像素的深度估计值,使得整个深度图中的不同部分内部较为平滑。

    深度图优化。原始代价体往往是含有噪声污染的,因此,为防止噪声使得网络过度拟合,MVSNet中使用基于多尺度的三维卷积神经网络进行代价体正则化,利用U-Net网络[7],对代价体进行降采样,并提取不同尺度中的上下文信息和临近像素信息,对代价体进行过滤(图5)。

    图5  MVSNet代价体正则化

    4. MVSNet实战操作

    首先感谢Yaoyao(香港科技大学)给出的开源代码和已经预处理好的模型和数据。其次,本节首先介绍MVSNet利用开源数据集中的数据进行深度估计和点云恢复,之后,回顾该网络数据IO相关的要求和数据文件夹结构,最后,介绍如何利用自己采集的数据进行深度估计和点云恢复。

    4.1 利用开源数据集进行深度估计

    1) 环境配置

    参考Yaoyao的github主页中installation即可完成环境配置。

    https://github.com/YoYo000/MVSNet

    2) 数据整理

    在文末分享的百度云盘中下载数据集preprocessed_inputs/dtu.zip和预训练好的网络models/tf_model_19307.zip。将tf_model解压,在其中的3DCNNs/中获得训练好的模型 model.ckpt-100000.data-00000-of-00001。

    将test.py 中的pretrained_model 地址改为tf_model中3DCNNs的地址。

    图7 更改预训练模型的地址

    解压下载好的dtu.zip,到用于深度估计的数据集。以scan10为例,该文件夹的结构如下图。

    图8 scan10文件夹结构

    3) 深度估计

    运行代码:

    python test.py --dense_folder TEST_DATA_FOLDER --regularization '3DCNNs' --max_w 1152 --max_h 864 --max_d 192 --interval_scale 1.06
    

    注意:

    • flag --dense_folder 要设定为scan10的地址。

    • flag –regularization表示正则化代价体的方式,MVSNet中使用的是3D卷积神经网络。

    • 图片的大小可以按照GPU的大小变更参数,但是需要时32的整数倍(特征提取时2D神经网络要求图像是32的整数倍)

    • Max_d 和interval_scale 建议先按照默认的要求,在使用自己的数据时,我们会给出调整的方式

    4) 查看深度估计结果

    深度估计结束后,会在dense_folder文件夹中生成文件夹depths_mvsnet/。利用visualize.py 即查看相关的深度图和概率图(图9),输入代码:

    python visualized.py XXXXXXX.pfm

    图9 MVSNet深度图估计结果

    5) 点云生成

    MVSNet只能生成深度图,需要借助其他点云融合方法。这里Yaoyao给出了一种方法——fusible.可以直接利用生成的深度图和RGB图片,生成稠密点云。

    • 安装fusible

    git clone https://github.com/YoYo000/fusibile

    cd fusible

    mkdir build && cmake .

    make .

    • 运行代码

    python depthfusion.py --dense_folder TEST_DATA_FOLDER --fusibile_exe_path FUSIBILE_EXE_PATH --prob_threshold 0.3
    

    注意:每一个深度图在估计的时候,都有一个对应的概率图,及表示该像素沿深度方向不同的概率,所以flag –prob_threhold 表示用于过滤深度的概率阈值,即如果当前深度图的某个像素的深度对应的概率低于该阈值,则在点云融合的时候,该像素的深度不被使用,即该像素不会被投影到三维空间点。这样可以对稠密点云进行过滤优化。

    在dense_folder中,之后会生成新的文件夹points_mvsnet。在该文件夹下,找到consisitency_Check-time/final3d_model.ply. 即为融合之后的三维点云,利用meshlab打开即可。图10是笔者实验过后的稠密点云。

    图10 深度图融合之后的稠密点云

    6) 流程回顾

    环境配置 -> 整理输入数据 -> 配准好flag,运行代码 -> 安装深度融合软件 -> 生成稠密点云

    4.2 数据IO要求及数据配置回顾

    MVSNet的整个流程完整、简单,再回顾一下整个流程中文件夹的结构。

    1) Input

    输入时,有一个初始的文件夹(--dense_folder),存放了两个文件夹cams,images 和一个pair.txt.(图8以scan10为例)

    其中images文件夹,包含了49张照片,是放置畸变校正之后RGB图像的地方, cams中也包含了49个.txt文件,每一个txt文件就是对应序号图片的相机内参数和外参数. 第三个pair.txt 保存了,当每一个照片为参考影像时, 与其最相似的前十张相片(具体的排序方式是:在稀疏重建过程中, 两张相片共视的三维点数, 共视的三维点越多,则排名越高.)

    图11 pair.txt文件

    逐行解释如下:

    10 选取排名较高的前十张影像 10号 得分 2346.41分 1号, 得分2036.53 9号 得分1243.89。所以,对于0号影像来说,筛选前5个类似影像为10号,1号,9号,12号,11号,从数据集找出对应的RGB图片可观察筛选情况(图12)。

    图12 0号(ref) 10号 1号 9号 12号 11号 影像

    如果将MVSNet暂时看成黑箱, 那么输入到MVSNet的input文件夹,结构一定要遵从上述结构(图8)。

    2) Output

    在深度估计结束后,  会在dense_folder下新建一个depths_mvsnet, 里面包含了生成的深度图,点云生成后,还有一个文件夹points_mvsnet,里面包含场景的三维点云。

    图13 output文件夹结构

    4.3 利用自采数据进行深度估计

    将MVSNet看做一个黑箱, 需要将自采数据整理成4-1中dense_folder的文件夹结构,即需要每一个包含经过畸变矫正的RGB图片的images文件夹、包含每一个相片的内外参数的cams文件夹,和一个pair.txt文件。

    一般的自采数据,可获得RGB图片,所以,需要再利用稀疏重建的原理获得相机的内外参数和相关中间输出。这里提供了一个完整的数据处理pipeline,第一,先利用COLMAP进行稀疏重建,第二,导出相机文件和稀疏场景的信息,第三利用Yaoyao提供的python代码,将COLMAP结果直接转换为MVSNet的输入格式。具体步骤如下:

    1) 稀疏重建

    利用COLMAP进行稀疏重建,根据上一篇文章, 稀疏重建结束后,可得到了如下的场景。

    图14 稀疏重建结果

    2) 图像纠正

    点击“reconstruction” 中的“dense reconstruction”. 打开稠密重建窗口,进行图像畸变纠正。选择dense_folder为worksapce地址,点击“undistortion”进行图像去畸变。(原理这里不再赘述)

    图15 图像纠正窗口

    再返回到dense_folder下,可以看到当前的文件夹结构(图16)。这里的sparse是自动创建的保存相机数据和稀疏场景结构的文件夹。

    图16 纠正之后的文件结构

    3) 导出相机参数和中间数据

    点击COLMAP中“file”的 “export data as TXT”, 并以2)中的sparse/文件夹作为导出地址,导出三个文件:camera.txt, images.txt, points3D.txt.

    图17 相机数据和中间数据(选中部分)

    • camera.txt是保存相机内参数的文件

    • images.txt是保存图片外参数和图片二维特征与三维空间点对应信息的文件

    • points3D.txt是保存三维空间点在世界坐标系下坐标、RGB值以及在各个影像上的轨迹(track)

    4) 结果转换

    感谢Yaoyao的结果转换代码(colmap2mvsnet.py),可以将COLMAP的结果直接转化成MVSNet输入需要的结果。执行代码为:

    python colmap2mvsnet.py –dense_folder our_dense_folder –max_d XXX –interval_scale XXX
    
    • Flag –dense_folder 需要以我们的dense_folder作为地址,里面必须包含sparse文件夹

    • --max_d 表示,最大估计的离散深度采样区间数,因为MVSNet是按照平面扫描原理进行深度估计的,所以深度是离散采样的,一般我们设定为192个深度采样区间。

    • --interval_scale表示每个深度区间的大小,默认为1.06(mm)。

    • 深度估计范围:已知深度最小值depth_min,则深度最大值满足:

    • 注意:我们需要估计自采数据的深度范围,已保证在深度采样区间内,能对目标场景进行有效的深度估计,举个例子,如果自采数据的深度范围为45cm – 80cm,那么我们用于深度估计的区间范围应该是35cm,从45 – 80cm。此时如果我们设定的深度区间为0 – 35cm,那么估计出来的深度图肯定是错误的。所以对于自采数据,大家需要尝试不同的深度区间,以找到合适的取值范围。

    • 在估计自采数据的深度范围时,需要修改Yaoyao的代码,需要的朋友们可以私戳笔者,获取修改之后的代码。

    5) 深度估计

    转换结束后,将dense_folder中的取出cams,images两个文件夹和pair.txt文件,放入一个新的文件夹中,该文件夹就可以作为MVSNet的输入。然后,参考4-1节即可利用自采数据进行深度估计和稠密重建。


    图18 自采数据运行结果

    原图,深度图,稠密点云

    6) 流程回顾

    利用自己采集的数据进行深度估计的步骤为:采集数据 -> 稀疏重建 -> COLMAP文件导出 -> 数据转换 -> 输入MVSNet  -> 输出深度图 -> 深度融合,恢复点云。

    5. 总结

    MVSNet是一种端到端的基于监督学习的深度估计网络,从2018年提出伊始,就作为一种比较流行的深度学习框架,其提出的深度代价匹配和构造方法能快速对MVS数据进行参考影像的深度估计,其精度相比较传统方法也有一定的提升。

    可是,其在正则化过程中仍旧消耗近11GB的内存,这使得许多人在使用MVSNet的时候受到限制,下一篇文章,将讲述如何利用循环神经网络和一种链式的Cost Volume构造方法减少MVSNet的GPU消耗,敬请期待。

    参考文献

    [1] Yao Yao, Luo Zixin, Li Shiwei, Fang Tian, Quan Long. MVSNet: Depth Inference for Unstructured Multi-View Stereo. European Conference on Computer Vision (ECCV)

    [2] Scharstein D , Szeliski R . A Taxonomy and Evaluation of Dense Two-Frame Stereo Correspondence Algorithms[J]. International Journal of Computer Vision, 2002, 47(1-3):7-42.

    [3] Han, X., Leung, T., Jia, Y., Sukthankar, R., Berg, A.C.: Matchnet: Unifying feature and metric learning for patch-based matching. Computer Vision and Pattern Recognition (CVPR) (2015)

    [4] Collins R T . A Space-Sweep Approach to True Multi-Image Matching[C] Computer Vision and Pattern Recognition, 1996. Proceedings CVPR '96, 1996 IEEE Computer Society Conference on. IEEE, 1996.

    [5] Yang R , Pollefeys M . Multi-Resolution Real-Time Stereo on Commodity Graphics Hardware[C] 2003 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2003. Proceedings. IEEE, 2003.

    [6] Furukawa Y , Hernández, Carlos. Multi-View Stereo: A Tutorial[J]. Foundations & Trends in Computer Graphics & Vision, 2015, 9(1-2):1-148.

    [7] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI) (2015)

    附录

    开源数据集下载链接:(感谢香港科技大学 姚遥开源的数据集链接)我们需要下载的是用于测试的数据集,在公众号「3D视觉工坊」后台回复「mvsnet」,即可获得下载链接。打开百度云链接时,点击 “mvsnet”,再点击“preprocessed_inputs”,下载其中的“dtu.zip”和“tankandtemples.zip”即可(图18)。

    图19  目标目录

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

    下载1

    在「3D视觉工坊」公众号后台回复:3D视觉即可下载 3D视觉相关资料干货,涉及相机标定、三维重建、立体视觉、SLAM、深度学习、点云后处理、多视图几何等方向。

    下载2

    在「3D视觉工坊」公众号后台回复:3D视觉github资源汇总即可下载包括结构光、标定源码、缺陷检测源码、深度估计与深度补全源码、点云处理相关源码、立体匹配源码、单目、双目3D检测、基于点云的3D检测、6D姿态估计汇总等。

    下载3

    在「3D视觉工坊」公众号后台回复:相机标定即可下载独家相机标定学习课件与视频网址;后台回复:立体匹配即可下载独家立体匹配学习课件与视频网址。

    重磅!3DCVer-学术论文写作投稿 交流群已成立

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

    同时也可申请加入我们的细分方向交流群,目前主要有3D视觉CV&深度学习SLAM三维重建点云后处理自动驾驶、CV入门、三维测量、VR/AR、3D人脸识别、医疗影像、缺陷检测、行人重识别、目标跟踪、视觉产品落地、视觉竞赛、车牌识别、硬件选型、学术交流、求职交流等微信群。

    一定要备注:研究方向+学校/公司+昵称,例如:”3D视觉 + 上海交大 + 静静“。请按照格式备注,可快速被通过且邀请进群。原创投稿也请联系。

    ▲长按加微信群或投稿

    ▲长按关注公众号

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

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

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

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

    展开全文
  • 为了方便大家了解基于多视图立体的三维重建技术,更重要的是能亲手利用开源数据集或者自己采集的影像跑一遍流程,进而对整个流程更为熟悉,本文整理了近年来几种经典的基于传统方法和基于深度学习方法的三维重建...
  • 多视图几何中的三维重建

    千次阅读 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...
  • 多视图几何三维重建

    千次阅读 2021-12-06 21:12:54
    基础知识 相机标定 相机成像原理:世界坐标系、相机坐标系、图像坐标系、像素坐标系之间的转换 基于单目视觉的三维重建算法综述 SIFT图像匹配算法 增量式SFM估计相机运动和三角测量进行重建 增量式SFM后的非线性优化...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,462
精华内容 1,784
关键字:

多视图几何三维重建