COLMAP: Structure-from-Motion Revisited 论文阅读

Abstract

从无序的图像集中进行三维重建时经常使用的是增量式的SFM。尽管增量式的重建系统在所有方面都取得了巨大进步，但鲁棒性、准确性、完整性和可扩展性仍然是构建真正通用pipeline的关键问题。本文提出了一种新的SfM技术，它改进了现有技术，朝着最终目标迈出了一步。

colmap: COLMAP

1. Introduction

多年来，来自无序图像的SFM经历了巨大的发展。早期的自校准度量重建系统为无序互联网图像集和城市场景的第一个系统奠定了基础。受这些作品的启发，越来越多的大规模重建系统已经被开发，从数十万到数百万到最近的一亿张互联网照片。目前已经提出了多种SfM策略，包括增量式、层次式和全局式。可以说，增量SfM是重建无序照片集最流行的策略。尽管它被广泛使用，我们仍然没有设计出真正通用的SfM系统。虽然现有系统极大地提高了现有技术的水平，但鲁棒性、准确性、完整性和可扩展性仍然是增量SfM中的关键问题，阻碍了其作为通用方法的使用。在本文中，我们提出了一种新的SfM算法来实现这一最终目标。新方法在各种具有挑战性的数据集上进行了评估，代码作为名为COLMAP的开源实现贡献给研究社区，可在https://github.com/colmap/colmap.

2. Review of Structure-from-Motion

SFM是将3D结构从其投影重建为从不同视点拍摄的一系列图像的过程。增量SFM（在本文中表示为SFM）是具有迭代重建组件的顺序处理流水线（图2）。它通常从特征提取和匹配开始，然后是几何验证。生成的场景图用作重建阶段的基础，该阶段在增量配准新图像、三角剖分场景点、过滤异常值以及使用束调整（BA）细化重建之前，通过仔细选择的双视图重建来为模型播种。以下章节详细阐述了这一过程，定义了本文中使用的符号，并介绍了相关工作。

2.1. Correspondence Search

输入：一系列输入图片；输出：经过几何校验后的图像匹配关系。

为了得到尽可能准确的匹配关系，该步骤中涉及特征提取，匹配以及几何校验。

特征提取：可以是任何一种特异性较强的特征，如SIFT（COLMAP默认），主要为后续的特征匹配服务；
匹配阶段：将输入的图像两两之间进行匹配（可以发现，这一步的时间复杂度非常大），得到潜在的场景重合部分；
几何校验：初始匹配的外点势必很多，此时需要滤出外点。因此，SfM通过尝试估计使用投影几何在图像之间映射特征点的变换来验证匹配。单应性矩阵H描述了捕捉平面场景的纯旋转或移动摄像机的变换。极线几何描述了移动摄像机通过基本矩阵E（校准）或基本矩阵F（未校准）的关系，并可以使用三焦点张量扩展到三个视图。通常情况下会用到基础矩阵（未标定）/本质矩阵（已标定）以及单应矩阵（纯旋转/共面）。图像经过上述三个步骤之后的输出为scene graph，即图像是节点，几何校验后的匹配对是边。

2.2. Incremental Reconstruction

输入：「scene graph」；输出：相机位姿以及3D路标点；

为了完成增量重建，需要完成初始化，图像注册，三角化以及BA优化这些步骤。

初始化：SfM在初始化时需要非常仔细的选择两帧进行重建；此时需要尽量选择scene graph中相机间可视区域多的两视角进行初始化，文中称这种选择增加了“redundancy”进而增加了重建的鲁棒性与精确性。
图像注册：即根据已有的SfM模型估计图片的位姿（更准确的是估计拍摄该图像相机的位姿）；根据未知的图像与地图中图像的2D-2D匹配，进而得到二者之间的2D-3D匹配，然后利用RANSAC-PNP解算出该相机位姿，完成图像注册。为了提高位姿结算精度以及可靠的三角化，本文设计了新颖的鲁棒后续帧选择策略，后续章节进行介绍。
三角化：新注册的图像需要对已有的场景点有足够多的观测，同时也可以通过三角化扩展场景点。现有的三角化方法鲁棒性不足且耗时严重。本文设计了一种全新的三角化策略，后续章节进行介绍。
BA优化：由于图像注册与三角化是分开进行的，但二者是强相关的两个过程：错误的位姿会导致错误的三角化点点产生，反之亦然。BA能够同时优化相机位姿以及地图点，使模型的redundancy更强。这里顺便提一下BA，优化目标是通过调整相机位姿与地图点位置使重投影误差最小：

EEE = ∑j\sum_j∑j ρj\rho_jρj(∣\mid∣∣\mid∣π\piπ(PcP_cPc, XkX_kXk) - xjx_jxj∣\mid∣∣22\mid_{2}^{2}∣22)

BA问题通常可以通过LM(Levenberg-Marquardt)算法进行求解。作者提到，当SfM面对网络图像时，特别对于那些几乎一样的图像时，上述优化过程会占用极长的时间。本文尝试对该问题进行解决，后续章节进行介绍。

3. Challenges

虽然目前最先进的SfM算法可以处理大规模互联网照片集合中多样化和复杂的图像分布，但它们在完整性和鲁棒性方面往往不能产生完全令人满意的结果。通常情况下，系统无法对大量经验上应该可配准的图像进行配准，或者由于错配准或漂移，系统产生了破碎的模型。首先，这可能是由于对应搜索产生了一个不完整的场景图，例如，由于匹配近似，因此既没有为完整的模型提供必要的连通性，也没有为可靠的估计提供足够的冗余。第二，这可能是由于场景结构缺失或不准确，导致重建阶段无法配准图像-图像配准和三角剖分是一种共生关系，图像只能注册到已有的场景结构中，而场景结构只能从注册的图像中三角剖分。在增量重建过程中，在每一步中最大化两者的准确性和完整性是SfM的一个关键挑战。在本文中，我们解决了这一挑战，并在提高效率的同时，在完整性、鲁棒性和准确性方面显著提高了当前技术状态(第5节)的结果。

4. Contributions

本节介绍了一种新算法，该算法改进了SfM中的主要挑战。首先，我们引入了一种几何验证策略，该策略使用信息增强场景图，从而提高初始化和三角剖分组件的鲁棒性。第二，最大化增量重建过程的鲁棒性和准确性的次最佳视图选择。第三，一种稳健的三角剖分方法，以降低的计算成本产生比现有技术显著更完整的场景结构。第四，迭代BA、retriangulation和异常值过滤策略，通过减少漂移效应显著提高完整性和准确性。最后，通过冗余视图挖掘为密集照片集提供更有效的BA参数化。这导致系统在鲁棒性和完整性方面明显优于现有技术，同时保持其效率。我们对比了对当前最先进的系统Bundler （开源和VisualSFM（闭源代码）的贡献。本系统作为开源实现发布。

4.1. Scene Graph Augmentation

提出了一种多模型几何验证策略，用适当的几何关系对场景图进行增强。对于输入的图像对，具体过程如下：

首先估计基础矩阵FFF(fundamental matrix)，如果此时至少有NFN_FNF个内点数(inliers)，那么就认为该图像对通过了几何验证。
然后估计单应矩阵HHH，通过其内点数(inliers)NHN_HNH进行分类；如果NHN_HNH/NFN_FNF < ϵHF\epsilon_{HF}ϵHF，就认为是一般场景(general scene)
对于已经标定的相机，估计本质矩阵EEE (essential matrix)，记内点数为NEN_ENE ；如果NEN_ENE/NFN_FNF > ϵEF\epsilon_{EF}ϵEF，就认为相机的内参是正确的(correct calibration)。
如果此时场景标定是正确的且是一般场景(correct calibration and NHN_HNH/NFN_FNF < ϵHF\epsilon_{HF}ϵHF)，通过分解本质矩阵(essential matrix EEE )得到相机位姿，然后对点三角化(triangulate)，计算三角化点点平均角度 αm\alpha_mαm(median triangulation angle )；随后利用αm\alpha_mαm判断场景是纯旋转(pure rotation, panoramic)还是平面(planar scenes, planar）。
网络图像可能存在的问题有水印(watermarks)/时间戳(timestamps)/镶边(frames), 简称WTFs问题。此时计算图像对的相似矩阵 SSS，记录图像对在图像边缘的内点数 NSN_SNS；若此时NSN_SNS/NFN_FNF > ϵSF\epsilon_{SF}ϵSF或者NEN_ENE/NEN_ENE > ϵSE\epsilon_{SE}ϵSE，则认为该场景为WTFs，此时该图像不加入场景图**「scene graph」**。

综上，对**「scene graph」**中的图像对分为：常规，全景，平面(general, panoramic, planar)。仅从非全景和已标定图像中选择参与重建。

注意：不对全景图像进行三角化，以避免退化点，可以提高三角化和后续图像注册的鲁棒性。

4.2. Next Best View Selection

参与重建的图像的选择对重建的效果的影响是非常大，若选择了不好的图像，可能会导致相机无法注册或者三角化失败。

通常的做法是选择能看到三角化点最多的图像进行注册。但仅数量多并不意味着好，已经通过实验验证**「PNP位姿估计的准确率不仅与观测点数量相关而且与观测点在图像上的分布有关」**：若点分布越均匀，位姿估计的准确率越高。

选择方法：

找到候选图片中至少看到NtN_tNt>0个三维点的图片(用一个feature tracks graph可以有效的统计这个信息），但是网上的数据集中看到同一个场景的图片非常多，从而导致feature tracks graph会非常密集，本文提出用一个多分辨率分析的方法高效的解决这个问题。

作者认为：可见点的数量越多且分布越均匀得分SSS越高，得分越高点的图像越会被先注册。得分的计算方法：
- 在x,y轴上分别划分出KlK_lKl个固定大小的网格,每个网格的可以有两个状态—empty || full.
- 在重建过程中,当一个网格中的点visible了,这个网格就从empty转变成full,并且得分SSS的值根据权重wlw_lwl增加.
当可见点的数量NtN_tNt<<Kl2K_l^2Kl2时，这个方法可能不能很好的描述分布情况，因为点可能都落在一个网格中。为了解决这个问题,提出一种多分辨率的方法（l = 1, 2, …, L），其中KlK_lKl=2l2^l2l，权重wlw_lwl=Kl2K_l^2Kl2(原文这里应该出错了)，按照图三的计算权重应该是wlw_lwl=KlK_lKl，最后把每个分辨率的结果加起来得到分数SSS。

Score左上Score_{左上}Score左上 = 1 ×\times× 2 + 4 ×\times× 4 + 6 ×\times× 8 = 66

Score右上Score_{右上}Score右上 = 1 ×\times× 2 + 4 ×\times× 4 + 16 ×\times× 8 = 146

Score左上Score_{左上}Score左上 = 4 ×\times× 2 + 6 ×\times× 4 + 6 ×\times× 8 = 80

Score左上Score_{左上}Score左上 = 4 ×\times× 2 + 16 ×\times× 4 + 16 ×\times× 8 = 200

4.3. Robust and Efficient Triangulation

transitive corrcorrespondence可以增加三角测量的completeness，也对后序图片的注册有帮助，基于appearance的近似匹配产生的图片对通常是baseline较小的(相机移动不大的)，因此利用transitive(传递性)可以使三角测量更加准确。

现存的一些多视角三角测量方法对于一定程度的外点具有鲁棒性，但是还没有方法能对含有高比例外点(outlier)的特征轨迹(feature tracks)进行处理，本文提出一种高效的、基于采样(sampling-base)的方法，从含有外点的特征轨迹(feature tracks)中稳健的估计所有点。

**Bundler：**在track中取样所有的点对儿，进行two-view三角测量，一旦找到一对儿点满足三角测量angle，就在这个track上进行multi-view triangulation，一旦track中的所有点都满足cheirality constrain，就接受这个track，这种方法有两个limitation：1. 对于外点不robust，因为它没法recover融合到track中的single point（因为是进行two-view）；2。对于所有pairwise点进行two-view 三角测量计算代价大。

本文使用了RANSAC对多帧观测进行三角化。记特征追踪 τ\tauτ = {TnT_nTn | nnn = 1…NTN_TNT}，它是由一系列的内点率(inliers) ϵ\epsilonϵ 未知的观测构成的集合。观测TnT_nTn包括归一化坐标 xnˉ\bar{x_n}xnˉ ∈\in∈ R2\Bbb R^2R2以及它对应的投影矩阵PnP_nPn ∈\in∈ SE(3)SE(3)SE(3) ，其中P=[RT−RTt]P = [R^T - R^Tt]P=[RT−RTt] ，旋转矩阵R∈SO(3)R \in SO(3)R∈SO(3)，平移矩阵t∈R3t \in \Bbb R^3t∈R3。最终目标是最大化三角化模型对观测的复合程度，于是三角化的点XabX_{ab}Xab为：

XabX_{ab}Xab ~ τ(xaˉ,xbˉ,Pa,Pb)\tau(\bar{x_a}, \bar{x_b}, P_a, P_b)τ(xaˉ,xbˉ,Pa,Pb)

其中 τ\tauτ 是任何一种三角化方法(本文使用了DLT方法)。一个比较好的三角化点需要满足两个条件：

足够大的三角化角度α\alphaα；

cos⁡α\cos{\alpha}cosα = $\frac{t_a - X_{ab}}{\mid\mid t_a - X_{ab}\mid\mid_2} $ ⋅\cdot⋅ $\frac{t_b - X_{ab}}{\mid\mid t_b - X_{ab}\mid\mid_2} $

三角化点深度$d_a 和和和d_b为正，且该点的重投影误差为正，且该点的重投影误差为正，且该点的重投影误差e_n小于阈值小于阈值小于阈值t$;

ddd = [p31[p_{31}[p31 p32p_{32}p32 p33p_{33}p33 p34]p_{34}]p34] [XabT[X_{ab}^T[XabT 1]T1]^T1]T ;

en=∣∣xnˉ−[xˊ/zˊyˊ/zˊ]∣∣2with[xˊyˊzˊ]=Pn[Xab1]e_n = \left|\left|\bar{x_n} - \left[ \begin{matrix} \acute{x} / \acute{z}\\ \acute{y} / \acute{z}\\ \end{matrix} \right] \right|\right|_2 with \left[ \begin{matrix} \acute{x} \\ \acute{y} \\ \acute{z} \\ \end{matrix} \right] = P_n \left[ \begin{matrix} X_{ab} \\ 1 \\ \end{matrix} \right] en=∣∣∣∣xnˉ−[xˊ/zˊyˊ/zˊ]∣∣∣∣2with⎣⎡xˊyˊzˊ⎦⎤=Pn[Xab1]
值得注意的是，三角化的过程中使用了RANSAC，即从上述特征追踪中随机选择2个点(一对点)进行三角化，统计内点数(即，重投影误差小于阈值记为内点)，选择内点数最多的那个作为三角化点。这个过程中存在一个问题：假如该点被追踪到了比较少的次数，此时随机采样会重复选择相同的一对点进行三角化，这样会造成不必要的资源消耗。本文设计了一种可以产生不重复随机组的random sampler；另外，为了能够保证以置信度为的概率至少选择一对无外点的最小解集，该过程至少重复进行KKK次。因为一个特征轨迹可能包含多个独立的点，我们递归地运行这个过程(这个过程应该是RANSAC)，从剩下的测量中删除共识集(consensus set)。循环的终止条件是最后一次的共识集大小小于3。

4.4. Bundle Adjustment

一般在图像注册和三角测量之后进行BA，增量式三维重建中，增加的点只会对局部产生影响，因此只对most-connected image 注册之后进行局部BA即可。同VisualSfM类似，只在模型增长达到某个比率时才进行全局BA。

Parameterization：用柯西核函数(Cauchy function)作为局部BA的损失函数；当规模只有几百个相机的时候，使用稀疏直接求解器(sparse direct solver)；当问题规模很大的时候使用PCG；文章中使用Ceres Solver(谷歌开发的非线性优化的库)，并提供在任何图像组合中共享任意复杂的相机模型的选项。对于无序的网络照片，由于估计依赖于纯自校正，文章使用了一个只有径向畸变参数的简单相机模型。
Filtering ：在使用了BA之后，要过滤掉一些不符合模型的观测(observations)。首先滤除掉重投影误差大的点，然后滤除triangulation angle小于minimum的点；在全局BA之后，我们还检查退化相机，例如全景或人工增强图像引起的退化相机(e.g. those caused by panoramas or artificially enhanced images)，焦距和畸变参数在BA中自由优化，对于未标定的相机固定在图像中心点，视场异常或畸变系数较大的摄像机被认为是不正确的估计，在全局BA之后将会过滤。
Re-Triangulation ：类似于VisualSfM，文章中执行重新三角测量(RT)来解释全局BA (pre-BA RT)之前的漂移效应。然而，BA常常显著改善相机和点的参数。因此，我们建议用一个额外的post-BA RT步骤来扩展非常有效的pre-BA RT。由于BA对位姿进行了优化，这使得之前由于不准确/错误位姿导致无法三角化的点有了能够被重新三角化的可能性，所以在BA之后执行一步重三角化(类比于BA之前的三角化)；另外也尝试merge tracks。
Iterative Refinement：受到外点(outliers)的影响，一次BA通常会导致很大一部分点被滤除，然后经过post-BA RT对地图点进行了扩充；因此，文章中提出循环运行BA、RT、fliting，直到the number of filtered observations and post-BA RT points diminishes.，一般是做完第二遍之后模型的完整性就会得到非常不错的提升。

4.5. Redundant View Mining

BA是SfM的瓶颈，本文针对“稠密数据采集”带来的问题设计了一种高场景覆盖度的图像聚类的BA参数化策略。

Ni et al.[43]将划分作为连接的摄像机和点参数图上的图切割问题，将摄像机和点划分为通过分隔变量连接的子图(submaps)。文章中提出的方法是受这些前人工作的启发。与Ni等人类似，我们将问题划分为子映射(submaps)，其内部参数被去除(whose internal parameters are factored out)。这个方法共有三个贡献：

利用SfM固有属性提出一个高效的相机分组方案
不像之前的方法将许多摄像机聚集到一个子图中，本文中的方法是将场景划分为许多小的、高度重叠的摄像机组。每组中的摄像机被分解为单个摄像机，从而降低了解决减少的摄像机系统的成本。
由于相机组非常小，高度重叠，，省去了之前方法(Ni et al.)中分隔符变量的优化

图像与地图点根据重建过程中是否受到最新帧点影响分为两类。对于比较大的重建问题，最新的重建帧仅会影响其邻近的图像以及地图点，很大一部分图像/点并不会发生改变。BA自然会对新扩展的部分进行更多优化，而其他部分只在漂移的情况下进行改进。

本文将不受影响的那部分定义很多组的集合 G={GrG_rGr∣r=1…NGN_GNG} ，其中 GrG_rGr 表示一组高重合度的相机，它被参数化为一个相机。受到影响的那些独立分组，以便对他们参数进行优化（受到影响有新加入的图像以及一组图像中有超过ϵr\epsilon_rϵr比例的图像的重投影误差(reprojection error )都超过了r像素）。

组中的图像应该尽可能冗余，图像之间共可见点的数量是描述它们相互作用程度的度量。另一个场景中有NXN_XNX个点，那么每张图像都可以表示为一个向量vi∈{0,1}NXv_i \in \lbrace0,1\rbrace^{N_X}vi∈{0,1}NX，其中图像iii可以看见点XnX_nXn时，viv_ivi中的第n位就为1，否则为0。图像a和b之间的相互作用程度是通过对它们的向量viv_ivi的位运算来计算的:
Vab=∣∣va∧vb∣∣/∣∣va∨vb∣∣V_{ab} = \mid\mid v_a \wedge v_b \mid\mid / \mid\mid v_a \vee v_b \mid\mid Vab=∣∣va∧vb∣∣/∣∣va∨vb∣∣
建组方式(To build groups):

将图像按照viv_ivi的大小(∣∣vi∣∣≥∣∣vi+1∣∣\mid\mid v_i\mid\mid \geq \mid\mid v_{i+1}\mid\mid∣∣vi∣∣≥∣∣vi+1∣∣)排序。
初始化GrG_rGr时，首先把排好序的图像中的第一张图像IaI_aIa移除，然后找到IbI_bIb，使得VabV_{ab}Vab最大。如果VabV_{ab}Vab > VVV并且∣Gr∣\mid G_r \mid∣Gr∣ < SSS，就从排好序的图像中移除IbI_bIb并添加到GrG_rGr中,否则初始化新的组。为了减少搜索IbI_bIb的时间，将搜索限制在±β\pm \beta±β范围内的KτK_{\tau}Kτ最近邻。

分组图像的BA的损失函数(cost function)为:
Eg=∑jρj(∣∣πg(Gr,Pc,Xk)−xj∣∣22)E_g = \sum_j \rho_j ( \mid\mid \pi_g (G_r, P_c, X_k) - x_j\mid\mid_2^2) Eg=j∑ρj(∣∣πg(Gr,Pc,Xk)−xj∣∣22)
其中相机外参Gr∈SE(3)G_r \in SE(3)Gr∈SE(3)，PcP_cPc是图像的位姿，在r组中的图像的投影矩阵(projection matrix)就是拼接组的和图像的位姿：Pcr=PcGrP_{cr} = P_cG_rPcr=PcGr。总的cost就是组的和非组的cost之和。PcP_cPc和GrG_rGr中的旋转变量使用四元数(quaternions)。

组的大小越大，性能优势越大。请注意，即使是两个群体的情况，文章中说也观察到计算上的好处。此外，性能收益取决于问题的大小，因为相机数量的减少对直接方法的立方计算复杂度的影响大于对间接方法的线性复杂度的影响。

5. Experiments

作者在各种各样的数据集上(17个数据集，包含144953张无序网络图)与增量式SFM(Bundler、VisualSFM)和全局式SFM系统(DISCO、Theia)进行对比。QUAD拥有真实的(ground-truth)相机位置。在所有的实验中，我们使用RootSIFT使用在不相关的数据集上训练的词汇树，将每张图像与它的100个最近的邻居进行匹配。为了确保不同方法之间的可比性，不计算他们的关联搜索时间(correspondence search)，机器的配置是2.7GHz、256GB RAM。

Next Best View Selection：作者使用了L=6的空间金字塔以及生成高斯分布的点来模拟。σ越大，µ越中心，对应的分布越均匀，正确地产生较高的分数。当分数很低时，分数主要由点的数量决定，否则则由它们在图像中的分布决定。
Robust and Efficient Triangulation：在Dubrovnik数据集上的一个实验(在由47M个验证匹配组成的2.9M特征轨迹上)评估我们的方法。与Bundler进行比较。作者把α\alphaα(4.3节中三角化时要求的最小的角度)设成2°2^°2°，把t(4.3节中重投影误差的阈值)设成8px，把ϵ0\epsilon_0ϵ0(RANSAC时inlier的初始的比例)设为0.03。为了避免组合爆炸，我们将RANSAC的迭代次数限制为10K次，即ϵmin\epsilon_{min}ϵmin约为0.02，η\etaη为0.999(η\etaη是RANSAC的成功率)。递归方法明显比非递归方法恢复更长的轨迹和更多的轨迹元素。注意，基于ransac的递归方法的点数越高，对应的跟踪长度就越小。基于ransac的方法产生稍差的轨迹(tracks)，但速度要快得多(10-40倍)。通过改变η\etaη值，很容易平衡速度和完整性。
Redundant View Mining：通过VVV(看4.5节，对照片分组时要求Vab>VV_{ab}>VVab>V)这个参数可以控制重建的时间和质量，随着VVV的减小可以减少运行时间，同时会增加重投影误差降低重建的质量。
System：系统整体上的运行时间Theia是最快的方法，而本文的方法比VisualSFM略差，比Bundler快50多倍。但是本文的方法结果最好。

6. Conclusion

本文提出了一种克服关键挑战的SfM算法，向通用的SfM算法系统又迈进了一步。提出的算法组件在完整性、鲁棒性、准确性和效率方面改进了现有的算法。在具有挑战性的大规模数据集上进行的综合实验证明了单个组件和整个系统的性能。整个算法作为开源实现向公众发布。

参考

[1] Schönberger, J. L., & Frahm, J. M. (2016). Structure-from-Motion Revisited. In IEEE Conference on Computer Vision & Pattern Recognition.

[2] 博客：论文笔记]Structure-from-Motion Revisited - 知乎 (zhihu.com)

[3] 博客：三维重建系列之COLMAP: STRUCTURE-FROM-MOTION REVISITED - 知乎 (zhihu.com)

[4] 博客：COLMAP: Structure-from-Motion Revisited 论文解读_枫仙森的博客-CSDN博客

[5] 博客：三维重建系列之COLMAP: Structure-from-Motion Revisited_Tom Hardy的博客-CSDN博客

[6] 博客： structure-from-motion revisited论文笔记_你吃过卤汁牛肉吗的博客-CSDN博客