OpenCV实现SfM(四):Bundle Adjustment_aipiano的博客-程序员资料

技术标签: 三维  计算机视觉  opencv  

Bundle Adjustment

在上一篇文章中,成功将三维重建扩展到了任意数量的图像,但是,随着图像的增多,累计误差会越来越大,从而影响最终的重建效果。要解决这个问题,需要用到Bundle Adjustment(下文简称BA)。
BA本质上是一个非线性优化算法,先来看看它的原型
min ⁡ x ∑ i ρ i ( ∣ ∣ f i ( x i 1 , x i 2 , . . . , x i k ) ∣ ∣ 2 ) \min_x \sum_i{\rho_i(||f_i(x_{i1}, x_{i2}, ..., x_{ik})||^2)} xminiρi(fi(xi1,xi2,...,xik)2)
其中 x x x是我们需要优化的参数, f f f一般称为代价函数(Cost Function), ρ \rho ρ为损失函数(Loss Function)。其中 f f f的返回值可能是一个向量,因此总的代价取该向量的2-范数。
对于三维重建中的BA,代价函数往往是反向投影误差,比如我们需要优化的参数有相机的内参(焦距、光心、畸变等)、外参(旋转和平移)以及点云,设图像 i i i的内参为 K i K_i Ki,外参为 R i R_i Ri T i T_i Ti,点云中某一点的坐标为 P j P_j Pj,该点在 i i i图像中的像素坐标为 p j i p_j^i pji,则可以写出反向投影误差
f ( K i , R i , T i , P j ) = π ( K i [ R i    T i ] P j ) − p j i f(K_i, R_i, T_i, P_j)=\pi(K_i[R_i\ \ T_i]P_j) - p_j^i f(Ki,Ri,Ti,Pj)=π(Ki[Ri  Ti]Pj)pji
上式中的 P j P_j Pj p j i p_j^i pji均为齐次坐标,其中 π \pi π为投影函数,有 π ( p ) = ( p x / p z ,   p y / p z ,   1 ) \pi(p)=(p_x/p_z,\ p_y/p_z,\ 1) π(p)=(px/pz, py/pz, 1).
而损失函数 ρ \rho ρ的目的是为了增强算法的鲁棒性,使得算法不易受离群点(Outliers)的影响,常见的有Huber函数、Tukey函数等,这些函数的图像如下
Loss Functions

若不使用损失函数,即 ρ ( x ) = x \rho(x)=x ρ(x)=x,那么就如上图中的黑色曲线,代价(Cost)随着误差以二次幂的速度增长,也许一个误差较大的点,就能左右算法的优化方向。其他的损失函数,可以发现,随着误差的增大,要么代价的增长是趋于线性的(Huber),要么干脆趋于不变(Tukey),这样就降低了误差较大的点对总代价的影响。损失函数实际上就是代价的重映射过程。

Ceres Solver

如何求解BA?总体思想是使用梯度下降,比如高斯-牛顿迭代、Levenberg-Marquardt算法等,由于BA还有自己的一些特殊性,比如稀疏性,在实现时还有很多细节需要处理,在此就不细说了。好在现在有很多用于求解非线性最小二次问题的库,文中使用的就是Google的一个开源项目——Ceres Solver.
Ceres Solver专为求解此类问题进行了大量的优化,有很高的效率,尤其在大规模问题上,其优势更加明显。并且,Ceres内置了一些常用的函数,比如对坐标的旋转以及各类损失函数,使其在开发上也比较高效。在官网上可以找到它的编译方法和使用教程,Windows用户可以在此找到配置好的VS工程。

编写代码

代码总体基本不变,我们只需要再添加一个函数用于BA即可,还有一点需要注意的是,Ceres Solver默认使用双精度浮点,如果精度不够可能导致计算梯度失败、问题无法收敛等问题,因此在原来的代码中,需要将原本问Point3f型的structure,改为Point3d类型。
首先定义一个代价函数

struct ReprojectCost
{
    
	cv::Point2d observation;

	ReprojectCost(cv::Point2d& observation)
		: observation(observation)
	{
    
	}

	template <typename T>
	bool operator()(const T* const intrinsic, const T* const extrinsic, const T* const pos3d, T* residuals) const
	{
    
		const T* r = extrinsic;
		const T* t = &extrinsic[3];

		T pos_proj[3];
		ceres::AngleAxisRotatePoint(r, pos3d, pos_proj);

		// Apply the camera translation
		pos_proj[0] += t[0];
		pos_proj[1] += t[1];
		pos_proj[2] += t[2];

		const T x = pos_proj[0] / pos_proj[2];
		const T y = pos_proj[1] / pos_proj[2];

		const T fx = intrinsic[0];
		const T fy = intrinsic[1];
		const T cx = intrinsic[2];
		const T cy = intrinsic[3];

		// Apply intrinsic
		const T u = fx * x + cx;
		const T v = fy * y + cy;

		residuals[0] = u - T(observation.x);
		residuals[1] = v - T(observation.y);

		return true;
	}
};

该代价函数就是之前所说的反向投影误差,参数分别为内参、外参还有点在空间中的坐标,最后一个参数用于输出,为反向投影误差。注意,为了使BA更高效可靠,外参当中的旋转部分使用的是旋转向量而不是旋转矩阵,这样不仅使优化参数从9个变为3个,还能保证参数始终代表一个合法的旋转(如果用矩阵,可能在优化过程中,正交性不再满足)。

接下来直接使用Ceres Solver求解BA,其中使用了Ceres提供的Huber函数作为损失函数

void bundle_adjustment(
	Mat& intrinsic,
	vector<Mat>& extrinsics, 
	vector<vector<int>>& correspond_struct_idx,
	vector<vector<KeyPoint>>& key_points_for_all,
	vector<Point3d>& structure
)
{
    
	ceres::Problem problem;

	// load extrinsics (rotations and motions)
	for (size_t i = 0; i < extrinsics.size(); ++i)
	{
    
		problem.AddParameterBlock(extrinsics[i].ptr<double>(), 6);
	}
	// fix the first camera.
	problem.SetParameterBlockConstant(extrinsics[0].ptr<double>());

	// load intrinsic
	problem.AddParameterBlock(intrinsic.ptr<double>(), 4); // fx, fy, cx, cy

	// load points
	ceres::LossFunction* loss_function = new ceres::HuberLoss(4);	// loss function make bundle adjustment robuster.
	for (size_t img_idx = 0; img_idx < correspond_struct_idx.size(); ++img_idx)
	{
    
		vector<int>& point3d_ids = correspond_struct_idx[img_idx];
		vector<KeyPoint>& key_points = key_points_for_all[img_idx];
		for (size_t point_idx = 0; point_idx < point3d_ids.size(); ++point_idx)
		{
    
			int point3d_id = point3d_ids[point_idx];
			if (point3d_id < 0)
				continue;

			Point2d observed = key_points[point_idx].pt;
			// 模板参数中,第一个为代价函数的类型,第二个为代价的维度,剩下三个分别为代价函数第一第二还有第三个参数的维度
			ceres::CostFunction* cost_function = new ceres::AutoDiffCostFunction<ReprojectCost, 2, 4, 6, 3>(new ReprojectCost(observed));

			problem.AddResidualBlock(
				cost_function,
				loss_function,
				intrinsic.ptr<double>(),			// Intrinsic
				extrinsics[img_idx].ptr<double>(),	// View Rotation and Translation
				&(structure[point3d_id].x)			// Point in 3D space
			);
		}
	}

	// Solve BA
	ceres::Solver::Options ceres_config_options;
	ceres_config_options.minimizer_progress_to_stdout = false;
	ceres_config_options.logging_type = ceres::SILENT;
	ceres_config_options.num_threads = 1;
	ceres_config_options.preconditioner_type = ceres::JACOBI;
	ceres_config_options.linear_solver_type = ceres::SPARSE_SCHUR;
	ceres_config_options.sparse_linear_algebra_library_type = ceres::EIGEN_SPARSE;

	ceres::Solver::Summary summary;
	ceres::Solve(ceres_config_options, &problem, &summary);

	if (!summary.IsSolutionUsable())
	{
    
		std::cout << "Bundle Adjustment failed." << std::endl;
	}
	else
	{
    
		// Display statistics about the minimization
		std::cout << std::endl
			<< "Bundle Adjustment statistics (approximated RMSE):\n"
			<< " #views: " << extrinsics.size() << "\n"
			<< " #residuals: " << summary.num_residuals << "\n"
			<< " Initial RMSE: " << std::sqrt(summary.initial_cost / summary.num_residuals) << "\n"
			<< " Final RMSE: " << std::sqrt(summary.final_cost / summary.num_residuals) << "\n"
			<< " Time (s): " << summary.total_time_in_seconds << "\n"
			<< std::endl;
	}
}

如果求解成功,会输出一些统计信息,其中最重要的两项分别是优化前的平均反向投影误差(Initial RMSE)和优化后的该值(Final RMSE)。

main函数完全不用变,只需在最后加入如下代码,对BA进行调用即可

Mat intrinsic(Matx41d(K.at<double>(0, 0), K.at<double>(1, 1), K.at<double>(0, 2), K.at<double>(1, 2)));
vector<Mat> extrinsics;
for (size_t i = 0; i < rotations.size(); ++i)
{
    
	Mat extrinsic(6, 1, CV_64FC1);
	Mat r;
	Rodrigues(rotations[i], r);

	r.copyTo(extrinsic.rowRange(0, 3));
	motions[i].copyTo(extrinsic.rowRange(3, 6));

	extrinsics.push_back(extrinsic);
}

bundle_adjustment(intrinsic, extrinsics, correspond_struct_idx, key_points_for_all, structure);

优化结果对比

Before

Before

After

After

Before

Before

After

After

Statistics

Statistics

从统计信息可以看出,最初的重建结果,反向投影误差约为3.6个像素,BA之后,反向投影误差降为1.4个像素,如果删除一些误差过大的点,再进行一次BA,反向投影误差往往能小于0.5个像素!

这次的代码与上一篇文章的几乎一样,唯独多出来的几个函数和修改也已在上文中列出,这次就不再单独提供代码了。要运行本文的代码,需要编译Ceres Solver,还需要依赖Eigen(一个线性代数库),详细过程在Ceres的官网上均有提及。

结语

由于个人时间有限,马上就要面临毕业论文和工作等问题,这可能是该系列最后一篇文章。在此感谢各位朋友的支持,谢谢。

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/AIchipmunk/article/details/52433884

智能推荐

MyBatis源码的学习(15)---sqlSession.selectList方法_itw_zhangzx02的博客-程序员资料_etsqlsession() .selectlist

sqlSession一共俩个实现类,我们这里分析默认的DefaultSqlSession类public &lt;E&gt; List&lt;E&gt; selectList(String statement, Object parameter, RowBounds rowBounds) { try { //ms对象代表我们的xml中的一条sql。例如:&lt;select...

最强Android教程!掌握这些Android开发热门前沿知识,赶紧收藏!_左夜天的博客-程序员资料_android 实用前沿技术

前言这些题目是网友去百度、小米、乐视、美团、58、猎豹、360、新浪、搜狐等一线互联网公司面试被问到的题目。熟悉本文中列出的知识点会大大增加通过前两轮技术面试的几率。欢迎一线公司员工以及网友提交面试题库,欢迎留言。网上的都是按照公司划分的,想找具体某一方面的知识点有点不好找,我这里就根据知识点分门别类的整理了一下,想看哪一块可以快速找到。阿里巴巴面试问题还记得一些,一部分已经忘记了,为了防止再忘记,所以写出来。1:你是如何理解Android操作系统的。2:是否熟悉framework层,如果熟

revit 转换ifc_导出 IFC 文件以使用 BIM 软件进行编辑_瑞士鲁迅的博客-程序员资料

内容在编辑 BIM IFC 文件之前,了解 Revit 导出设置。用户建筑师、工程师、设计师、项目经理、BIM 经理。方法1. 定义要导出的图元在 Revit 项目中,定义专用于导出 IFC 文件的三维视图。在此视图中,您可以控制对编辑文件的团队成员可见的图元。请勿包含编辑软件不支持的类别,例如 HVAC 图元。最好在协调模型中共享这些类别。2. 检查 IFC 导出映射文件选择“文件”&gt;“导...

linux系统实验shell程序设计,实验三-shell脚本程序设计_李苦李的博客-程序员资料

《实验三-shell脚本程序设计》由会员分享,可在线阅读,更多相关《实验三-shell脚本程序设计(10页珍藏版)》请在人人文库网上搜索。1、实 验 报 告课程名称 Linux系统实践 实验项目 LINUX SHELL脚本程序设计 实验仪器 PC 系 别 计算机学院 专 业 网络工程 班级/学号 网1702/2017011463 学生姓名 孟启贤 实验日期 4.15 成 绩 指导教师 李艳平 实验...

Best Practices for Speeding Up Your Web Site_longware_新浪博客_龙威的博客-程序员资料

Minimize HTTP RequestsUse a Content Delivery NetworkAdd an Expires or a Cache-Control HeaderGzip ComponentsPut Stylesheets at the TopPut Scripts at the BottomAvoid CSS Expre...

【安卓开发系列 -- APP】APP 性能优化 -- 启动优化_奋斗企鹅CopperSun的博客-程序员资料_应用启动性能优化

【安卓开发系列 -- APP】APP 性能优化 -- 启动优化【1】APP 启动优化的必要性原因 : 用户希望应用能够及时响应并快速加载,启动时间长的 APP 不能满足该期望;启动太慢的结果 : 体验效果差,用户流失,产品失败;【2】启动流程以及分类【2.1】开机启动流程注意,该流程图所示的启动过程为系统创建并启动应用的过程,一般不需要优化;【2.2】冷启动冷启动指应用从头开始启动,系统进程在冷启动后才创建应用进程,发生冷启动的情况包括应用自设备启动后或系统终止应用后.

随便推点

PostMessage()_丫丫afc的博客-程序员资料

PostMessage函数PostMessage是WindowsAPI(应用程序接口)中的一个常用函数,用于将一条消息放入到消息队列中。消息队列里的消息通过调用GetMessage和PeekMessage取得。函数功能该函数将一个消息放入(寄送)到与指定窗口创建的线程相联系消息队列里,不等待线程处理消息就返回,是异步消息模式。消息队列里的消息通过调用GetMessage和PeekMess...

如何在Ubuntu中安装rpm包_This is bill的博客-程序员资料

原本是为了装flash player但是下载下来的是rpm的包, 所以有这个问题 解决方案Normally you install software or deb packages on Ubuntu/Mint linux via Synaptic, Ubuntu Software Center/ppa, or an apt-get command from the termina...

解决raise child_exception OSError: [Errno 2] No such file or directory的方法_DCGKCUF的博客-程序员资料_exception oserror: (2, 'no such file or directory'

如果你在执行python相关的程序时,出现如下的错误.....python2.5/subprocess.py", line 1079, in _execute_child     raise child_exception OSError: [Errno 2] No such file or directory那你就一定得检查与此相关的命令是否在你的$PATH路径下,如果不在,请一

一网打尽:贝佐斯与亚马逊时代_历史可以倒流的博客-程序员资料

全书主要描述了亚马逊从1995年刚开始成立时的一家网上书店到发展到现在成为集网上购物,云计算,Kindle业务的一家巨无霸的看似曲折却又似乎是必然的过程。说这个过程曲折,那是因为Amazon在每次想要逃脱外界定义的只是一家在网上卖书的店这个定义时,都不是一帆风顺的,无论是来自外界的阻力,还是内部的压力。而说他必然是因为贝佐斯这个整个亚马逊公司的大脑,一直在追求客户至上,扩张优先,把亚马逊是一家技术

opencart是PHP框架吗,OpenCart框架运行流程介绍_味离的博客-程序员资料

框架运行流程介绍这样的一个get请求http://hostname/index.php?route=common/home 发生了什么?1. 开始执行入口文件index.php。2. require_once(DIR_SYSTEM . 'startup.php');做一些php的配置和加载一些类声明,包括系统主框架文件(system/engine下的文件)、一些必用到的helper和library...

Highcharts 前端图表插件_阿远个人博客的博客-程序员资料

Highcharts 支持将图片下载成各种格式 Highcharts官网:https://www.hcharts.cn/download<div id="container" style="min-width:400px;height:400px"></div><script src="https://img.hcharts.cn/jquery/jquery-1.8.3.min.js"></

推荐文章

热门文章

相关标签