OpenCV实现SfM(二):双目三维重建[通俗易懂]

OpenCV实现SfM(二):双目三维重建[通俗易懂]使用OpenCV3.0实现双目三维重建,原理清晰,实践有效。

大家好,又见面了,我是你们的朋友全栈君。

注意:本文中的代码必须使用OpenCV3.0或以上版本进行编译,因为很多函数是3.0以后才加入的。
目录:

文章目录

#极线约束与本征矩阵

在三维重建前,我们先研究一下同一点在两个相机中的像的关系。假设在世界坐标系中有一点 p p p,坐标为 X X X,它在1相机中的像为 x 1 x_1 x1,在2相机中的像为 x 2 x_2 x2(注意 x 1 x_1 x1 x 2 x_2 x2为齐次坐标,最后一个元素是1),如下图。
这里写图片描述
X X X到两个相机像面的垂直距离分别为 s 1 s_1 s1 s 2 s_2 s2,且这两个相机具有相同的内参矩阵 K K K,与世界坐标系之间的变换关系分别为 [ R 1    T 1 ] [R_1\ \ T_1] [R1  T1] [ R 2    T 2 ] [R_2\ \ T_2] [R2  T2],那么我们可以得到下面两个等式
s 1 x 1 = K ( R 1 X + T 1 ) s 2 x 2 = K ( R 2 X + T 2 ) s_1x_1 = K(R_1X + T_1) \\ s_2x_2 = K(R_2X + T_2) s1x1=K(R1X+T1)s2x2=K(R2X+T2)
由于K是可逆矩阵,两式坐乘K的逆,有
s 1 K − 1 x 1 = R 1 X + T 1 s 2 K − 1 x 2 = R 2 X + T 2 s_1K^{-1}x_1 = R_1X + T_1 \\ s_2K^{-1}x_2 = R_2X + T_2 s1K1x1=R1X+T1s2K1x2=R2X+T2
K − 1 x 1 = x 1 ′ K^{-1}x_1 = x_1^{‘} K1x1=x1 K − 1 x 2 = x 2 ′ K^{-1}x_2 = x_2^{‘} K1x2=x2,则有
s 1 x 1 ′ = R 1 X + T 1 s 2 x 2 ′ = R 2 X + T 2 s_1x_1^{‘} = R_1X + T_1 \\ s_2x_2^{‘} = R_2X + T_2 s1x1=R1X+T1s2x2=R2X+T2
我们一般称 x 1 ′ x_1^{‘} x1 x 2 ′ x_2^{‘} x2为归一化后的像坐标,它们和图像的大小没有关系,且原点位于图像中心。
由于世界坐标系可以任意选择,我们将世界坐标系选为第一个相机的相机坐标系,这时 R 1 = I ,   T 1 = 0 R_1 = I,\ T_1 = 0 R1=I, T1=0。上式则变为
s 1 x 1 ′ = X s 2 x 2 ′ = R 2 X + T 2 s_1x_1^{‘} = X \\ s_2x_2^{‘} = R_2X + T_2 s1x1=Xs2x2=R2X+T2
将第一式带入第二式,有
s 2 x 2 ′ = s 1 R 2 x 1 ′ + T 2 s_2x_2^{‘} = s_1R_2x_1^{‘} + T_2 s2x2=s1R2x1+T2
x 2 ′ x_2^{‘} x2 T 2 T_2 T2都是三维向量,它们做外积(叉积)之后得到另外一个三维向量 T 2 ^ x 2 ′ \widehat{T_2}x_2^{‘} T2
x2
(其中 T 2 ^ \widehat{T_2} T2
为外积的矩阵形式, T 2 ^ x 2 ′ \widehat{T_2}x_2^{‘} T2
x2
代表 T 2 × x 2 ′ T_2\times x_2^{‘} T2×x2),且该向量垂直于 x 2 ′ x_2^{‘} x2 T 2 T_2 T2,再用该向量对等式两边做内积,有
0 = s 1 ( T 2 ^ x 2 ′ ) T R 2 x 1 ′ 0 = s_1(\widehat{T_2}x_2^{‘})^TR_2x_1^{‘} 0=s1(T2
x2)TR2x1


x 2 ′ T 2 ^ R 2 x 1 ′ = 0 x_2^{‘}\widehat{T_2}R_2x_1^{‘} = 0 x2T2
R2x1=
0

E = T 2 ^ R 2 E = \widehat{T_2}R_2 E=T2
R2

x 2 ′ E x 1 ′ = 0 x_2^{‘}Ex_1^{‘} = 0 x2Ex1=0
可以看出,上式是同一点在两个相机中的像所满足的关系,它和点的空间坐标、点到相机的距离均没有关系,我们称之为极线约束,而矩阵 E E E则称为关于这两个相机的本征矩阵。如果我们知道两幅图像中的多个对应点(至少5对),则可以通过上式解出矩阵 E E E,又由于 E E E是由 T 2 T_2 T2 R 2 R_2 R2构成的,可以从E中分解出 T 2 T_2 T2 R 2 R_2 R2
如何从 E E E中分解出两个相机的相对变换关系(即 T 2 T_2 T2 R 2 R_2 R2),背后的数学原理比较复杂,好在OpenCV为我们提供了这样的方法,在此就不谈原理了。

#特征点提取与匹配
从上面的分析可知,要求取两个相机的相对关系,需要两幅图像中的对应点,这就变成的特征点的提取和匹配问题。对于图像差别较大的情况,推荐使用SIFT特征,因为SIFT对旋转、尺度、透视都有较好的鲁棒性。如果差别不大,可以考虑其他更快速的特征,比如SURF、ORB等。
本文中使用SIFT特征,由于OpenCV3.0将SIFT包含在了扩展部分中,所以官网上下载的版本是没有SIFT的,为此需要到这里下载扩展包,并按照里面的说明重新编译OpenCV(哎~真麻烦,-_-!)。如果你使用其他特征,就不必为此辛劳了。
下面的代码负责提取图像特征,并进行匹配。

void extract_features(
	vector<string>& image_names,
	vector<vector<KeyPoint>>& key_points_for_all,
	vector<Mat>& descriptor_for_all,
	vector<vector<Vec3b>>& colors_for_all
	)
{ 
   
	key_points_for_all.clear();
	descriptor_for_all.clear();
	Mat image;

	//读取图像,获取图像特征点,并保存
	Ptr<Feature2D> sift = xfeatures2d::SIFT::create(0, 3, 0.04, 10);
	for (auto it = image_names.begin(); it != image_names.end(); ++it)
	{ 
   
		image = imread(*it);
		if (image.empty()) continue;

		vector<KeyPoint> key_points;
		Mat descriptor;
		//偶尔出现内存分配失败的错误
		sift->detectAndCompute(image, noArray(), key_points, descriptor);

		//特征点过少,则排除该图像
		if (key_points.size() <= 10) continue;

		key_points_for_all.push_back(key_points);
		descriptor_for_all.push_back(descriptor);

		vector<Vec3b> colors(key_points.size());
		for (int i = 0; i < key_points.size(); ++i)
		{ 
   
			Point2f& p = key_points[i].pt;
			colors[i] = image.at<Vec3b>(p.y, p.x);
		}
		colors_for_all.push_back(colors);
	}
}

void match_features(Mat& query, Mat& train, vector<DMatch>& matches)
{ 
   
	vector<vector<DMatch>> knn_matches;
	BFMatcher matcher(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 * max(min_dist, 10.0f)
			)
			continue;

		//保存匹配点
		matches.push_back(knn_matches[r][0]);
	}
}

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

得到匹配点后,就可以使用OpenCV3.0中新加入的函数findEssentialMat()来求取本征矩阵了。得到本征矩阵后,再使用另一个函数对本征矩阵进行分解,并返回两相机之间的相对变换R和T。注意这里的T是在第二个相机的坐标系下表示的,也就是说,其方向从第二个相机指向第一个相机(即世界坐标系所在的相机),且它的长度等于1。

bool find_transform(Mat& K, vector<Point2f>& p1, vector<Point2f>& p2, Mat& R, Mat& T, Mat& mask)
{ 
   
	//根据内参矩阵获取相机的焦距和光心坐标(主点坐标)
	double focal_length = 0.5*(K.at<double>(0) + K.at<double>(4));
	Point2d principle_point(K.at<double>(2), K.at<double>(5));

	//根据匹配点求取本征矩阵,使用RANSAC,进一步排除失配点
	Mat E = findEssentialMat(p1, p2, focal_length, principle_point, RANSAC, 0.999, 1.0, mask);
	if (E.empty()) return false;

	double feasible_count = countNonZero(mask);
	cout << (int)feasible_count << " -in- " << p1.size() << endl;
	//对于RANSAC而言,outlier数量大于50%时,结果是不可靠的
	if (feasible_count <= 15 || (feasible_count / p1.size()) < 0.6)
		return false;

	//分解本征矩阵,获取相对变换
	int pass_count = recoverPose(E, p1, p2, R, T, focal_length, principle_point, mask);

	//同时位于两个相机前方的点的数量要足够大
	if (((double)pass_count) / feasible_count < 0.7)
		return false;

	return true;
}

#三维重建
现在已经知道了两个相机之间的变换矩阵,还有每一对匹配点的坐标。三维重建就是通过这些已知信息还原匹配点在空间当中的坐标。在前面的推导中,我们有
s 2 x 2 = K ( R 2 X + T 2 ) s_2x_2 = K(R_2X + T_2) s2x2=K(R2X+T2)
这个等式中有两个未知量,分别是 s 2 s_2 s2 X X X。用 x 2 x_2 x2对等式两边做外积,可以消去 s 2 s_2 s2,得
0 = x 2 ^ K ( R 2 X + T 2 ) 0 = \widehat{x_2}K(R_2X+T_2) 0=x2
K(R2X+
T2)

整理一下可以得到一个关于空间坐标X的线性方程
x 2 ^ K R 2 X = − x 2 ^ K T 2 \widehat{x_2}KR_2X = -\widehat{x_2}KT_2 x2
KR2X=
x2
KT2

上面的方程不能直接取逆求解,因此化为其次方程
x 2 ^ K ( R 2    T ) ( X 1 ) = 0 \widehat{x_2}K(R_2\ \ T)\left(\begin{matrix}X \\ 1\end{matrix}\right) = 0 x2
K(R2  T)(X1)=
0

用SVD求X左边矩阵的零空间,再将最后一个元素归一化到1,即可求得X。其几何意义相当于分别从两个相机的光心作过 x 1 x_1 x1 x 2 x_2 x2的延长线,延长线的焦点即为方程的解,如文章最上方的图所示。由于这种方法和三角测距类似,因此这种重建方式也被称为三角化(triangulate)。OpenCV提供了该方法,可以直接使用。

void reconstruct(Mat& K, Mat& R, Mat& T, vector<Point2f>& p1, vector<Point2f>& p2, Mat& structure)
{ 
   
	//两个相机的投影矩阵[R T],triangulatePoints只支持float型
	Mat proj1(3, 4, CV_32FC1);
	Mat proj2(3, 4, CV_32FC1);

	proj1(Range(0, 3), Range(0, 3)) = Mat::eye(3, 3, CV_32FC1);
	proj1.col(3) = Mat::zeros(3, 1, CV_32FC1);

	R.convertTo(proj2(Range(0, 3), Range(0, 3)), CV_32FC1);
	T.convertTo(proj2.col(3), CV_32FC1);

	Mat fK;
	K.convertTo(fK, CV_32FC1);
	proj1 = fK*proj1;
	proj2 = fK*proj2;

	//三角化重建
	triangulatePoints(proj1, proj2, p1, p2, structure);
}

#测试
我用了下面两幅图像进行测试
这里写图片描述

得到了着色后的稀疏点云,是否能看出一点轮廓呢?!

这里写图片描述
这里写图片描述

图片中的两个彩色坐标系分别代表两个相机的位置。
在接下来的文章中,会将相机的个数推广到任意多个,成为一个真正的SfM系统。

关于源代码的使用
代码是用VS2013写的,OpenCV版本为3.0且包含扩展部分,如果不使用SIFT特征,可以修改源代码,然后使用官方未包含扩展部分的库。软件运行后会将三维重建的结果写入Viewer目录下的structure.yml文件中,在Viewer目录下有一个SfMViewer程序,直接运行即可读取yml文件并显示三维结构。

代码下载

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/151282.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

  • java如何实现封装_java如何实现封装

    java如何实现封装_java如何实现封装Java中类的封装是如何实现的封装是将对象的信息隐藏在对象内部,禁止外部程序直接访问对象内部的属性和方法。java封装类通过三个步骤实现:(1)修改属性的可见性,限制访问。(2)设置属性的读取方法。(3)在读取属性的方法中,添加对属性读取的限制。Java中什么叫封装呢?继承和多态都明白些,就是封装理解不上去,老师没关于这个问题,我想举一个例子:lz如果你接触过老的面向过程的编程,以前…

  • Softmax函数求导

    Softmax函数求导来源:https://blog.csdn.net/zt_1995/article/details/62227603其实整个推导,上面这个图片已经介绍得十分清楚了,但是仍有很多小步骤被省略掉了,我会补上详细的softmax求导的过程:…

  • 35 个 Java 代码性能优化总结 10-20「建议收藏」

    35 个 Java 代码性能优化总结 10-20

  • pycharmimport时找不到指定文件_pycharm系统找不到指定文件

    pycharmimport时找不到指定文件_pycharm系统找不到指定文件1、现象系统提示找不到指定的文件:Errorrunning’hello’:Cannotrunprogram"B:\pystudy\venv\Scripts\python.exe"(indirectory"\python-study"):CreateProcesserror=2,系统找不到指定的文件。2、原因原来的工程目录(B盘)下,保存了python的编…

  • java集合类框架的基本接口有哪些

    java集合类框架的基本接口有哪些转自:牛客网java集合类框架的基本接口有哪些?答:Collection:代表一组对象,每一个对象都是它的子元素Set:不包括重复元素的CollectionList:有顺序的Collection,并且可以包含重复元素Map:可以把键(key)映射到值(value)的对象,键不能重复下面是详细解释:转自:牛客网(一)总共有两大接口:Collecti

  • 最新手机号码归属地数据库(2017年4月1日)

    最新手机号码归属地数据库(2017年4月1日)2017年4月1日版近36万条记录celldb.cc最新号码归属地数据库手机号段数据库移动联通电信移动号段联通号段电信号段虚拟170号段171号段号码字段包括省市运营商邮编区号等信息移动号码:134135136137138139147150151152157158159178182183184187188联通…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号