Colmap中重三角化方法总结

这篇文章回答以下问题：

**Colmap中两种重三角化是如何实现的？**

Colmap作为现在最先进，最完整的一个SFM+MVS系统，其中有很多地方值得学习和借鉴。对于其中的SFM模块，采用了大量的tricks来提高SFM重建的准确性和完整性，本文试图解释是其中之一：重三角化。

本文只讨论重三角化，其在全局BA前后完成，目的是为了提高重建的完整性和准确性。主要参考资料，来自Colmap工程代码和论文《Structure from motion Revisited》。

目录

重三角化

在原论文中，有一段文字提到重三角化，其大致翻译如下：

重三角化 和VisualSFM类似，本文也使用了在全局BA前重三角化来减小漂移效应。但是，BA总是可以很好的优化相机和点的参数。因此，本文提出拓展有效的pre-BA 重三角化，增加post-BA 重三角化。目的在于提高重建的完整性，通过继续先前由于不准确的相机位姿导致三角化失败的点的追踪。我们只继续追踪哪些误差在滤波阈值下的观察点。此外，我们尝试合并追踪，因此提高下一步BA的完整性

论文中提到，重三角化可以提高重建三维点的完整性，并且减小累积误差带来的漂移。Colmap扩展了重三角化到全局BA之后，所以使用了两次重三角化。

这里听起来似乎，两种BA似乎使用同样的计算流程去计算，但是事实上并非如此。那么两种重三角化是如何实现的呢？这就需要看代码细节了。

pre-BA RT: 全局BA前的重三角化

在colmap源代码中sfm/incremental_triangulation.cc，很容易找到一个ReTriangulate函数，其在IterrativeGlobaBA（见controllers/incremental_mapper.cc）中使用，如下：

void IterativeGlobalRefinement(
const IncrementalMapperOptions& options,
IncrementalMapper* mapper) {
  PrintHeading1("Retriangulation");
  CompleteAndMergeTracks(options, mapper);
  std::cout << "  => Retriangulated observations: "
        << mapper->Retriangulate(options.Triangulation()) << std::endl;

/*..全局BA优化代码..*/ 
}

可以总结这个重三角化的大致计算流程如下：

1. 设置最大三角化误差为5.0

2. 对于每个图像对，只在under-reconstructed 图像对之间使用重三角化：

   已有三维点/匹配特征点匹配数量<0.2

   2.1 检测是否是under-reconstructed 图像？否则跳过

   2.2 图像是否在已经在重建中使用？否则跳过

   2.3 图像相机参数是否有问题？否则跳过

   2.4 对于每个图像对中的每个特征点匹配对

   ​ 2.4.1 **都有对应三维点？**是则跳过（无论是同一三维点/不同三维点，这里不处理后一种情况）。

   ​ 2.4.2 如果两个特征点都没有三维点，使用EstimateTriangulation (estimators/triangulate.cc)，这是在论文中介绍过的采样三角化方法，在这里不解释，理解为普通三角化就好。如果误差满足（注意此处阈值为5.0），则添加该3D点和观测关系到重建中。

   ​ 2.4.3 如果其中一个特征点有对应重建三维点，设置另外一个特征点的3D点为该点，添加观测关系到重建中。

post-BA RT：全局BA后的重三角化

BA后的重三角化不是很好理解，并且在代码中也没有直接找到这样的实现。基于colmap.cc的提示，和Github上Issue#414作者本人回答，post-BA重三角化在这个函数中实现:

/////////////////////////////////////////////////////////////////////////
// Retriangulation
/////////////////////////////////////////////////////////////////////////
PrintHeading1("Retriangulation");
CompleteAndMergeTracks(mapper_options, &mapper);

回去看IterativeGlobalBA函数：

void IterativeGlobalRefinement(const IncrementalMapperOptions& options,
                               IncrementalMapper* mapper) {
  PrintHeading1("Retriangulation");
  //重三角化
  CompleteAndMergeTracks(options, mapper);
  std::cout << "  => Retriangulated observations: "
            << mapper->Retriangulate(options.Triangulation()) 
            << std::endl;
    
  for (int i = 0; i < options.ba_global_max_refinements; ++i) {//每次BA迭代
    const size_t num_observations =
        mapper->GetReconstruction().ComputeNumObservations();
    size_t num_changed_observations = 0;
    AdjustGlobalBundle(options, mapper);
    //重三角化
    num_changed_observations += CompleteAndMergeTracks(options, mapper);
    //滤除部分重投影误差较大的点
    num_changed_observations += FilterPoints(options, mapper);
    const double changed =
        static_cast<double>(num_changed_observations) /num_observations;
    std::cout << StringPrintf("  => Changed observations: %.6f",changed)
              << std::endl;
    if (changed < options.ba_global_max_refinement_change) {
      break;
    }
  }

  FilterImages(options, mapper);
}

如果说CompleteAndMergeTracks(options, mapper)中有post-BA RT的实现，就很让人疑惑，在pre-BA RT前居然也来了一次post-BA RT？ 可能这样效果比较好，这个问题暂时先搁置。

可以看到，这个函数中BA使用前后确实各使用了RT，并且每次迭代前后都使用post-BA RT，滤去外点，个人认为这是很工程化的做法，有助于提高重建模型的精度和完整性。

下面来看一下post-BA RT的具体实现方法：

CompleteTracks：

对于一些追踪成功的特征点，在某些图像中没有被三角化，如果误差满足要求，则这些特征点添加为对应三维的观察值，具体程序步骤如下：

1. 设置最大三角化误差为4.0，pre-BA RT为5.0（issue中作者回答：RT with tighter threshold）

2. 对于3D点的追踪元素（理解为3D点对应某图像某个特征点在其它图像中特征点的对应关系）

   2.1 对于每个图像对中的每个特征点匹配对

   ​ 2.1.1 图像是否在已经在重建中使用？图像相机参数是否有问题？否则跳过

   ​ 2.1.2 有对应三维点？是则跳过

   ​ 2.1.3 计算重投影误差是否在阈值范围内？否则跳过

   ​ 是则继续，添加该图像特征点为3D点的观测（增加约束，增加完整性）

   ​

MergeTracks：

合并不同图像中同一特征点的不同观察，用一个基于追踪长度权重的公式计算，如果所有追踪点误差都满足，则合并。具体程序步骤如下：

1. 设置最大三角化误差为4.0，pre-BA RT为5.0

2. 对于3D点的追踪元素（理解为3D点对应某图像某个特征点在其它图像中特征点的对应关系）

   2.1 对于每个图像对中的每个特征点匹配对

   ​ 2.1.1 图像是否在已经在重建中使用？图像相机参数是否有问题？否则跳过

   ​ 2.1.2 有对应三维点？是则跳过

   ​ 2.1.3 2D特征点是否存在对应三维点？否则跳过该三维点是否已经合并过？是则跳过

   ​ 如果继续：

   ​ a. 基于追踪长度合并三维点，即两个三维点各自对应观测点个数，合并公式为：
   $mergedXYZ=(point3D*length1+Corrpoint3D*length2)/(length1+length2)$mergedXYZ=(point3D∗length1+Corrpoint3D∗length2)/(length1+length2)
   ​ b. 计算合并值，和每个追踪2D点在图像上重投影误差。

   ​ 全部满足阈值要求，则合并成功

   ​ c. 合并成功，删除初始点。

总结

回到两种重三角化方法，全局BA前的重三角化重构了一些地图点；BA后的RT没有增加新的地图点，而是给地图点增加了更多的观测（2D图像点），甚至更高的阈值可能还会减少地图点。Colmap中两种方法的结合，我认为其实就是一个Coarse to Fine的过程，先增加地图点，然后增加约束数量，滤去外点，再做BA优化，从而获得更好的重建结果。这些想法都是作者不断调试所获得方法，具体应用重三角化到自己的场景中，也需要更多的调试来获得更适合的处理方法。

以上colmap中两种重三角化方法，是VisualSFM中RT的延伸，具体VisualSFM如何实现，我们无从得知。但是重三角化作为一种减小增量式sfm减小误差的方法，还是很值得借鉴，潜在的用法是，我们可以将其用于特征点法SLAM全局地图的构建时，减小累积误差。