Ceres Solver实现简单的光束法平差

为什么使用Ceres Solver

Ceres Solver是谷歌开发的一款用于非线性优化的库，在SFM（Structure From Motion）中的最后一步，对定向和三角化后的像片位姿、相机内参以及三维空间点坐标进行Bundle Adjustment，优化像片位姿（旋转矩阵、平移矩阵，在COLMAP中，旋转矩阵以四元数的形式表示），相机内参矩阵（包括焦距、像主点、径向畸变、切向畸变），三维空间点位。

光束法平差原理简单理解

由于同名像点前方交会可能并不相交，所以取平均作为最终重建三维点，由该点及像片位姿后方交会生成的像点坐标与一开始的像点坐标不重合，因此存在一个反向投影误差，以最小化反向投影误差来迭代优化像片位姿以及空间三维坐标点。这便是光束法平差的最基本思想。

下图应该是网络上出现频率较高的：但这只是理想的前方交会结果，即两条直线刚好交会到同一个点，比如射线P1X1和P2X1交会于X1，实际并不交会，而是取平均后的结果，假设未X1附近的一个点X1'，这样反向计算出来的像点不再是原始的像点，而是像点附近的一个点，就存在了反向投影误差，即两像点之间的差值，如第二幅图所示。

（图片来源：https://blog.csdn.net/OptSolution/article/details/64442962）

光束法平差单元

按照自己的理解，Ceres Solver进行光束法平差的函数模板中存在一个平差单元，其中cost_function的定义会涉及到这个平差单元。根据上一部分原理的简单介绍，把每张像片上观测到的像点（即SFM流程中通过SIFT匹配到的可以前方交会出空间三维点的像点），定义为Observed；其所在的像片的位姿，定义为qvec和tvec；像片对应的相机内参，定义为camera_param；该像点对应的空间三维点，定义为point_3D。上述四部分便构成了一个最基本的光束法平差单元。因为同一个空间三维点可以在多张像片上找到对应的像点，以及整个平差涉及到多张像片，因此存在多个这样的基本平差单元。

代码分析

原始代码分析：

struct SnavelyReprojectionError {SnavelyReprojectionError(double observed_x, double observed_y): observed_x(observed_x), observed_y(observed_y) {}template <typename T>bool operator()(const T* const camera,const T* const point,T* residuals) const {// camera[0,1,2] are the angle-axis rotation.T p[3];ceres::AngleAxisRotatePoint(camera, point, p);// camera[3,4,5] are the translation.p[0] += camera[3];p[1] += camera[4];p[2] += camera[5];// Compute the center of distortion. The sign change comes from// the camera model that Noah Snavely's Bundler assumes, whereby// the camera coordinate system has a negative z axis.T xp = -p[0] / p[2];T yp = -p[1] / p[2];// Apply second and fourth order radial distortion.const T& l1 = camera[7];const T& l2 = camera[8];T r2 = xp*xp + yp*yp;T distortion = 1.0 + r2  * (l1 + l2  * r2);// Compute final projected point position.const T& focal = camera[6];T predicted_x = focal * distortion * xp;T predicted_y = focal * distortion * yp;// The error is the difference between the predicted and observed position residuals[0] = predicted_x - observed_x;residuals[1] = predicted_y - observed_y;return true;}// Factory to hide the construction of the CostFunction object from// the client code.static ceres::CostFunction* Create(const double observed_x,const double observed_y) {return (new ceres::AutoDiffCostFunction<SnavelyReprojectionError, 2, 9,  3>(new SnavelyReprojectionError(observed_x, observed_y)));}double observed_x;double observed_y;
};

以上部分定义了反向投影误差是怎么计算的，即residuals的计算，以及损失函数是怎么定义的，损失函数中的几个参数，第一个为代价函数的类型，第二个为代价的维度，第三个为每张像片外方位元素和相机内参，demo使用的相机内参为一个焦距加两个畸变参数，再加上3个旋转矩阵参数3个平移矩阵参数，正好9个，第四个为每一个空间三维点的坐标为3维。

if (!bal_problem.LoadFile("C:\\Users\\sggzg\\Desktop\\problem-49-7776-pre.txt")) {std::cerr << "ERROR: unable to open file " << "C:\\Users\\sggzg\\Desktop\\problem - 49 - 7776 - pre.txt" << "\n";}
const double* observations = bal_problem.observations();// Create residuals for each observation in the bundle adjustment problem. The
// parameters for cameras and points are added automatically.ceres::Problem problem;for (int i = 0; i < bal_problem.num_observations(); ++i) {
// Each Residual block takes a point and a camera as input and outputs a 2
// dimensional residual. Internally, the cost function stores the observed
// image location and compares the reprojection against the observation.ceres::CostFunction* cost_function =SnavelyReprojectionError::Create(observations[2 * i + 0],observations[2 * i + 1]);
problem.AddResidualBlock(cost_function,NULL /* squared loss */,bal_problem.mutable_camera_for_observation(i),bal_problem.mutable_point_for_observation(i));}// Make Ceres automatically detect the bundle structure. Note that the
// standard solver, SPARSE_NORMAL_CHOLESKY, also works fine but it is slower
// for standard bundle adjustment problems.ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_SCHUR;
options.max_num_iterations = 3;
options.minimizer_progress_to_stdout = true;ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
std::cout << summary.FullReport() << "\n";

这一部分中for循环主要就是添加每一个前面所说的光束法平差单元，AddResidualBlock函数即添加每一个平差单元的需要优化的参数，Create函数通过具体的像点观测值实例化前面定义的损失函数。平差的一些选项设置，平差报告的打印。

结合自己的实验数据进行调整：

主要是对反向投影误差计算的修改，根据自己的误差定义修改计算residuals的过程，以及根据相机的不同重新定义损失函数，如下：

ceres::CostFunction* cost_function = new ceres::AutoDiffCostFunction<BundleAdjustment_CostFunction, 2, 4, 3, 3, 3>(new BundleAdjustment_CostFunction(point2D.X(), point2D.Y()));problem.AddResidualBlock(cost_function, loss_function, psfmreconstruction->ImageDatas[img_idx].image.Qvec().data(),psfmreconstruction->ImageDatas[img_idx].image.Tvec().data(), psfmreconstruction->params_.data(), psfmreconstruction->points3Ds_[point2D.Point3DId()].XYZ().data());

这时可以看出损失函数中的几个参数，第一个为代价函数的类型，第二个为代价的维度，第三个为每张像片外方位元素和相机内参，我使用的相机内参为一个焦距加两个像主点参数，再加上4元数表示的qvec4个参数以及平移向量3个参数，第四个为每一个空间三维点的坐标为3维。

实验数据

原始代码针对的数据集下载地址：http://grail.cs.washington.edu/projects/bal/

该数据集为txt格式，内容截图如下：

第一行列出像片数num_image、空间三维点数num_point、所有的像点观测值数num_observed；

第二行至第num_observed行，依次为该像点对应的相机号，像点对应的空间三维点号、像点坐标值x、像点坐标值y；

第num_observed行至第（image_param_num*num_image+num_point*3）行，（此时image_param_num已经包括了每张像片的位姿以及对应的相机内参）为像片位姿、对应相机内参以及空间所有的三维点坐标值。

加载txt数据后，demo实验结果为：

自己的数据实验：

自己实现SFM过程中，初始化、增量式定向以及三角化后求解的像片位姿、空间三维坐标点、给定的相机内参并不会按照官网给的demo那样按照指定的格式输出到文本文件中，而是直接存储到对应的变量中。只需要将这些变量作为实参传递给前面分析的对应函数中。

以及优化后的几张像片位姿（qvec和tvec，qvec以四元数形式表示）

Ceres Solver实现简单的光束法平差相关推荐

Bundle Adjustment 光束法平差详解
转自http://blog.csdn.net/junshen1314/article/details/48860951 首先引述来自维基百科的定义:假设我们有一个3D空间中的点,他被位于不同位置的多个 ...
sba(sparse bundle adjustment):一个基于Levenberg-Marquardt(LM)算法的通用稀疏光束法平差C/C++软件包
[转]sba(sparse bundle adjustment):一个基于Levenberg-Marquardt(LM)算法的通用稀疏光束法平差C/C++软件包 2011-01-04 10:44 转载 ...
win7下用VS编译SBA（摄影测量光束法平差程序库）
我在学习区域网平差过程,自己编写平差程序,在解算方程过程中,变量的改正量总是达不到规定的阈值,偶然间发现sba,简直不能再爽了,sba是稀疏光束网平差库,可以根据自己的情况去修改代码,得到自己的BA程 ...
“光束法”和“空中三角测量”的辨析
"光束法"和"空中三角测量"的辨析光束法光束法(Bundle Adjustment)全称:光束法平差:又因广泛用于空中三角测量而称"光束法空中三角 ...
Ceres Solver: 高效的非线性优化库（二）实战篇
Ceres Solver: 高效的非线性优化库(二)实战篇接上篇: Ceres Solver: 高效的非线性优化库(一) 如何求导 Ceres Solver提供了一种自动求导的方案,上一篇我们已经看 ...
Ceres Solver Document学习笔记
Ceres Solver Document学习笔记 Ceres Solver Document学习笔记 1. 基本概念 2. 基本方法 2.1 CostFunction 2.2 AutoDiffCos ...
Ceres Solver 非线性优化库
Ceres Solver 非线性优化库 1. Ceres Solver 2. 下载安装 3. 简易例程 4. 环境运行 5. 非线性拟合 1. Ceres Solver Ceres solver 是谷 ...
Ceres Solver从零开始手把手教学使用
目录一 .简介二.安装三.介绍四.Hello Word! 五.导数 1 数值导数 2解析求导六.实践 Powell函数一 .简介笔者已经半年没有更新新的内容了,最近学习视觉SLAM的过程 ...
【环境配置】ceres solver安装
1. 安装 github地址 # CMake sudo apt-get install cmake # google-glog + gflags sudo apt-get install libgoo ...
R语言时间序列（time series）分析实战：简单指数平滑法预测
R语言时间序列(time series)分析实战:简单指数平滑法预测目录

Ceres Solver实现简单的光束法平差

Ceres Solver实现简单的光束法平差相关推荐

最新文章

热门文章