Bundle Adjustment

Google 做 Ceres 的一个主要原因就是用来解决大规模的 Boundle Adjustment 问题 HartleyZisserman], Triggs].

给定一组测量到的图像特征点的位置和对应关系， Boundle Adjustment 的目标是找到能使重投影误差最小的三维点位置和相机参数。这个优化问题通常被表述为一个非线性最小二乘法问题，其中误差是观测到的特征位置与相应三维点在相机图像平面上的投影之间差值的 square $L_2$ norm。

Ceres 给出的例子是解决 BAL 数据集中的问题。和前面的例子一样，第一步是定义一个 template functor 来计算重投影误差/残差。

BAL 问题中的每个残差都取决于一个三维点和一个相机，其中定义相机的九个参数是：三个旋转参数（罗德里格斯轴角向量）、三个平移参数、一个焦距参数和两个径向畸变参数。有关该摄像机模型的详细信息，请参阅 Bundler homepage and the BAL homepage。

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 - T(observed_x);
    residuals[1] = predicted_y - T(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>
       (observed_x, observed_y);
   }

  double observed_x;
  double observed_y;
};

请注意，与前面的例子不同，计算该残差函数的解析 Jacobian 有点麻烦。而自动微分功能则使计算变得简单得多。函数 AngleAxisRotatePoint() 和其他用于操作旋转的函数可以在 include/ceres/rotation.h 中找到。

bundle adjustment 问题建模如下：

ceres::Problem problem;
for (int i = 0; i < bal_problem.num_observations(); ++i) {
  ceres::CostFunction* cost_function =
      SnavelyReprojectionError::Create(
           bal_problem.observations()[2 * i + 0],
           bal_problem.observations()[2 * i + 1]);
  problem.AddResidualBlock(cost_function,
                           nullptr /* squared loss */,
                           bal_problem.mutable_camera_for_observation(i),
                           bal_problem.mutable_point_for_observation(i));
}

bundle adjustment 问题构造与曲线拟合示例非常相似，每个观测值都向目标函数中添加了一个残差项。由于这是一个大型稀疏问题（无论如何对 DENSE_QR 来说是大型问题），解决这个问题的一种方法是将 Solver::Options::linear_solver_type 设置为 SPARSE_NORMAL_CHOLESKY 并调用 Solve() 。虽然这样做是合理的，但 bundle adjustment 问题具有特殊的稀疏性结构，利用这种结构可以更高效地解决这些问题。Ceres 提供了三种专门的求解器（统称为 Schur-based 的求解器）来完成这项任务。示例代码使用了其中最简单的 "DENSE_SCHUR"。

ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_SCHUR;
options.minimizer_progress_to_stdout = true;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
std::cout << summary.FullReport() << "\n";

有关更复杂的 BA 示例，请参阅 examples/bundle_adjuster.cc ，该示例展示了 Ceres 更高级功能的使用，包括各种线性求解器、鲁棒核函数和流形。

Footnotes

• examples/simple_bundle_adjuster.cc