自动化学报  2018, Vol. 44 Issue (3): 506-516   PDF    
基于凸优化改进的相机全局位置估计方法
谢理想1, 万刚1, 曹雪峰1, 王庆贺1, 王龙2     
1. 解放军信息工程大学 郑州 450000;
2. 国防信息学院 武汉 430000
摘要: 相机全局位置估计作为运动恢复结构算法(Structure from motion,SfM)中的核心内容一直以来都是计算机视觉领域的研究热点.现有相机全局位置估计方法大多对外点敏感,在处理大规模、无序图像集时表现的尤为明显.增量式SfM中的迭代优化步骤可以剔除大部分的误匹配从而降低外点对估计结果的影响,而全局式SfM中没有有效地剔除误匹配的策略,估计结果受外点影响较大.针对这种情况,本文提出一种改进的相机全局位置估计方法:首先,结合极线约束提出一种新的对误匹配鲁棒的相对平移方向估计算法,减少相对平移方向估计结果中存在的外点;然后,引入平行刚体理论提出一种新的预处理方法将相机全局位置估计转化为一个适定性问题;最后,在此基础上构造了一个对外点鲁棒的凸优化线性估计模型,对模型解算获取相机位置估计全局最优解.本文方法可以很好地融合到当下的全局式SfM流程中.与现有典型方法的对照实验结果表明:在处理大规模、无序图像时,本文方法能显著提高相机全局位置估计的鲁棒性,并保证估计过程的高效性和估计结果的普遍精度.
关键词: 运动恢复结构     全局位置估计     平行刚体     凸优化     极线几何    
An Improved Method for Camera Location Estimation Through Convex Optimization
XIE Li-Xiang1, WAN Gang1, CAO Xue-Feng1, WANG Qing-He1, WANG Long2     
1. PLA Information Engineering University, Zhengzhou 450000;
2. National Defense Information College, Wuhan 430000
Manuscript received : September 7, 2016, accepted: January 16, 2017.
Foundation Item: Supported by National Natural Science Foundation of China (41301428, 4110143), Open Research Fund of State Key Laboratory of Geographical Information Engineering (SKLGIE2016-M-3-4)
Author brief: XIE Li-Xiang Master student at PLA Information Engineering University. His research interest covers computer vision and 3D reconstruction;
CAO Xue-Feng Lecturer at PLA Information Engineering University. He received his Ph. D. degree in cartography and GIS from PLA Information Engineering University in 2012. His research interest covers point cloud 3D reconstruction, simultaneous localization, and mapping;
WANG Qing-He Master student at PLA Information Engineering University. His research interest covers intelligent control of UAVs and cooperative for multiple UAVs;
WANG Long Assistant lecturer at National Defense Information College. He received his bachelor degree and master degree from PLA Information Engineering University in 2012 and 2015, respectively. His research interest covers data mining, large data analysis, and data visualization
Corresponding author. WAN Gang Professor at PLA Information Engineering University. He received his Ph. D. degree in cartography and GIS from PLA Information Engineering University in 2006. His research interest covers UAV surveying and mapping, virtual geographic environment. Corresponding author of this paper
Recommended by Associate Editor JIA Yun-De
Abstract: As a core module of structure from motion (SfM), location estimation of cameras in a global framework has been a research hotspot of computer version. State-of-the-art methods for location estimation are sensitive to outliers, especially for large scale, unordered images. The incremental SfM reduces the influence of outliers through an iterative optimization. The global SfM does not have an efficient strategy to remove mismatch, so the result of estimation is influenced deeply by outliers. Therefore, we introduce an improved method for location estimation. First, combined with the epipolar constraint we propose a new pairwise direction estimation algorithm. Then, we make the problem well-posed by introducing a new preprocessing method based on parallel rigidity. Finally, we propose a robust linear estimation model based on convex programing. We can get a global optimum solution by resolving this model. The method can integrate well with state-of-art global SfM pipeline. Multiple group experiments have proved the robustness of our methods without any loss of efficiency and common precision.
Key words: Structure from motion (SfM)     global position estimating     parallel rigidity     convex optimization     epipolar geometry    

运动恢复结构算法(Structure from motion, SfM)是多视图三维重建技术的关键内容, 一直都是计算机视觉领域的基础问题. SfM在虚拟现实、增强现实、运动追踪、逆向工程、城市场景建模等方面都有广泛应用. SfM算法解决的是一个从二维图像恢复拍摄时的相机运动和场景三维稀疏结构的问题. SfM算法通常由以下三个步骤组成: 1)图像特征匹配及相机相对位姿估计(相对平移和相对旋转); 2)相机运动估计, 即根据相机相对位姿估计相机的全局旋转和全局位置. 3)利用经过最小化重投影误差优化的相机参数进行场景稀疏结构恢复, 常用光束法平差[1] (Bundle adjust, BA)作为优化算法.在目前SfM算法中对第1)和第3)步骤的研究的较为成熟, 出现了较为完备的理论和算法.相对而言, 对第2)步骤的研究尚不够成熟, 尤其是相机全局位置估计部分, 现有的位置估计方法多对外点敏感.

现有的SfM算法可大致分为两类.一类是增量式方法[2-8] (Incremental method), 选取两张图像作为初始像对进行初始重建, 然后不断加入新的图像以扩大重建的范围, 直到所有的图像全部重建完毕.在这个迭代过程中, 通常会进行多次BA处理以优化相机运动参数及重建的场景三维稀疏结构.另一类是全局式方法[9-11] (Global method), 相对于增量式方法来说, 全局式方法不是一个逐渐累加的过程, 而是针对所有图像同时进行估计, 一次性解算所有的相机运动参数, 参数解算完毕后再进行场景三维稀疏结构的估计, 最后只进行一次全局性的BA优化, 最终得到相机参数和稀疏场景结构.

增量式方法首先利用初始像对估计相机运动和局部场景三维稀疏结构, 后续的迭代重建是基于初始像对的重建结果进行的, 最终的计算结果对初始像对选取的依赖性较高.此外, 在不断的迭代BA过程中, 误差会不断累积, 导致后续图像重建的计算误差大于前者, 最终可能产生场景漂移现象.全局式方法对所有的图像同时进行处理, 只进行一次的BA优化.这不仅能有效地避免误差累积, 将误差平均分配, 并且省去了耗时反复的迭代BA过程.另外全局式方法中的相机运动估计和场景三维稀疏结构恢复是分开进行的, 在所有相机运动全部估计完毕后再进行场景三维稀疏结构恢复, 大幅减少了每次需要估计的参数个数, 这也从另一方面提升了算法效率.

虽然全局式SfM算法特点突出, 但现有的SfM算法中, 增量式方法明显多于全局式方法.这是由于现有的相机运动估计算法大多对外点敏感, 这里的外点主要包括前文步骤1)中图像特征匹配过程中产生的误匹配值以及在相机运动估计过程中由于误匹配存在产生的外点(例如相机相对平移方向估计结果中的外点).增量式方法通过反复迭代BA过程可以有效地剔除大部分的误匹配, 降低外点对估计结果的影响.而全局式方法并没有这一剔除过程, 当外点较多时估计的精度无法保证.这种情况在处理大规模、无序图像时表现的尤为明显.因此对外点鲁棒的相机运动估计算法具有很高的研究价值, 可以充分发挥全局式SfM方法的优势, 提高SfM处理效率和精度.

在全局式SfM过程中相对运动估计分为相机全局旋转和相机全局位置估计两部分, 通常这两部分开进行.近20年, 相机运动估计的研究取得了明显的进展, 其中相机全局旋转估计方面已经出现不少成熟稳定的算法[12-15].然而, 在相机全局位置估计方面, 尤其是在处理大规模、无序图像集时, 现有的相机全局位置估计算法存在几个问题:估计结果不稳定[7, 16-17]、对相对平移方向中的外点的敏感[18-19]以及非凸方程易陷入局部最优解[11].因此, 鲁棒、高效并且能保证收敛于全局最优解的全局位置估计方法就显得很重要.

早期的全局位置估计算法[16, 20]是基于相机间相对平移方向构建线性方程并求取方程的最小二乘解来进行全局位置的估计.这类线性方法具有很高的计算效率, 但后续研究发现通过这类方法获取的解并不稳定, 在少数位置可能会产生偏离真值的伪解(Spectral solution).文献[7]试图通过估计相机之间相对平移的比例来消除这种伪解, 但并未有效提升估计结果的精度.文献[11]提出Lie代数平均法, 可有效解决伪解的问题, 但因为采用非凸的优化方程, 导致估计结果收敛到局部最优解而无法保证全局最优.文献[18]提出了一种基于范数的拟凸优化方法.这种方法可以很好地解决局部最优解的问题, 但范数通常易受外点影响, 估计结果难以保证.同本文方法最相似的是文献[10]提出的方法, 文献[10]构建了一个基于几何距离约束的线性方程, 但采用的是范数约束, 导致对外点敏感.另一个和本文相关的是文献[21]提出的基于相对平移方向构成线性方程求取最小化范数解的方法, 文献[21]通过加入约束来有效消除子伪解问题, 但因相对平移方向观测值中外点的影响, 估计结果的精度会降低.另外, 文献[22]中提出了一种基于经典刚体平移理论的全局位置估计算法.同文献[22]相比, 本文采用的是平行刚体理论, 使得估计方程截然不同.文献[23]提出了一种称为1DSfM的方法, 首先通过预处理剔除相对平移方向测量中的外点, 再设计一个非凸的优化方程对全局位置进行估计.

针对现有全局式SfM算法中全局位置估计存在的不足, 结合平行刚体理论、凸优化方法以及极线几何理论, 本文提出一种新的高效鲁棒的全局位置估计方法.并且本文方法可以很好地融合于现有的全局式SfM算法流程中.本文内容结构如下:第1节简单介绍了全局式SfM算法.第2节和第3节是对相机全局位置估计方法的改进.第2节阐述了一种结合极线约束的相对平移方向估计算法.第3节引入平行刚体理论, 阐述一种基于凸优化的鲁棒性相机全局位置估计算法.第4节是实验部分, 使用8组公开数据集对本文算法测试并和相关典型算法进行比较.

1 全局式SfM算法简介

由于本文工作是在全局式SfM方法基础上对其中的相机全局位置估计部分进行了改进, 所以首先对全局式SfM算法进行简要介绍.

全局式SfM一般包含以下5个步骤:

步骤1. 图像特征提取与匹配.

步骤2. 相机相对位姿估计.

步骤3. 相机全局旋转估计

步骤4. 相机全局位置估计.

步骤5. 光束法平差(BA).

其中相机全局位置估计步骤包含相机相对平移方向估计和基于相对平移方向的全局位置解算两部分.由于在全局式SfM处理过程中, 只进行了一次BA优化, 中间过程中没有类似增量式方法的迭代BA处理, 不能有效地剔除误匹配的影响.误匹配的存在会直接导致相机相对平移方向估计结果中的外点产生.而相机相对平移方向估计结果中存在的外点会直接影响相机全局位置估计的最终精度, 甚至导致计算结果无法收敛.

因此增加相机全局位置估计方法对外点的鲁棒性, 提高估计结果的精度和稳定性可以从两个方面考虑:减少误匹配对相机相对平移估计过程的影响, 削减相机相对平移方向估计结果中的外点, 以及在相对平移方向估计结果存在外点的情况下进行相机全局位置的鲁棒性估计.针对这两方面, 本文对全局式SfM中的相机全局位置估计方法做出以下两点改进:首先提出一种新的基于极线约束的相机相对平移估计方法.然后引入平行刚体理论将相机全局位置估计问题转换为一个适定性问题并提出一种基于凸优化的能收敛到全局最优的鲁棒性相机位置估计方程.并且改进后的相机全局位置估计方法可以很好地融合到传统的全局式SfM流程中, 如图 1所示.图 1中的虚线框部分为本文所做的改进.

图 1 融合本文改进的全局式SfM算法流程图 Figure 1 The processing pipeline of global SfM fusion our modifying
2 基于极线约束的相对平移方向估计

相对平移方向估计指的是对全局坐标系下相机之间相对平移方向的估计.全局式SfM算法中, 相机相对平移方向的估计是相机全局位置估计的基础, 相对平移方向估计的结果直接影响到了后续的全局位置估计的精度.

通常, 估计相对平移方向是通过对本质矩阵 $E_{ij}$ 分解获得相对旋转参数 $R_{ij}$ 和平移参数 $\pmb t_{ij}$ (这里的相对平移和旋转是在相机坐标系下的).然后相对旋转参数 $R_{ij}$ 估算出全局旋转 $R_{i}$ , 继而据此解算出相对平移方向估计值$ ({\pmb \gamma_{ij}}={R_i {\pmb t_{ij}}/||\pmb t_{ij}||}) $.在误匹配存在或者匹配点过少的情况下, 利用这种方法估计出的相对平移方向通常和真值存在较大偏差.全局式SfM算法中没有同增量式方法中一样的迭代优化过程, 所以误匹配点存在并参与运算的问题无法避免, 因此使用这种估计方法对全局式SfM过程中相对平移方向进行估计通常不能获得一个很精确的结果, 当处理大规模、无序图像集时这种情况表现得尤为明显.

对此, 本文提出一种基于极线约束的相对平移方向估计方法, 旨在当存在图像误匹配情况时仍能保持估计结果的鲁棒性(减少估计结果中存在的外点, 提高第3节中的全局位置估计精度).算法流程如图 2所示.

图 2 本文改进的相对平移方向估计算法流程 Figure 2 The processing pipeline of relative translation estimation based on our modifying

全局旋转的估计会在很大程度上影响到相对平移方向的估计精度.文献[14]提出了一种高效、对外点鲁棒的全局旋转估计算法, 其主要思想是首先利用一个速度相对较慢但是高鲁棒的 $L_1$ 范数优化方程最小化相对旋转误差, 然后使用一个速度较快的基于 $L_2$ 范数的优化方程对结果进行进一步的优化.因此, 为了提高相对平移方向估计的精度, 本文首先用文献[14]方法估计出全局旋转 $R_i$ , 然后结合鲁棒性的估计模型解算相对平移方向值.

为了表述方便, 用$(\{I_i\}_{i=1}^n) $表示图像集合. $\pmb P$ 表示三维场景点. $ (\{R_i\}_{i=1}^n) $$ (\{\pmb t_i\}_{i=1}^n)$$ (\{f_i\}_{i=1}^n) $分别表示各个相机的全局旋转、位置以及相机焦距. $\pmb P_i$ 表示场景点在以第 $i$ 个相机为中心的坐标系下的位置, $\pmb p_i$ 表示场景点投影到第 $i$ 个相机对应的图像平面上的位置. $\pmb P$ $\pmb P_i$ $\pmb p_i$ 之间可通过空间相似变化和投影变换相互转换, 具体转换公式如下:

$ \pmb P_i=R_i({\pmb P}-{\pmb t_i})=(X_i, Y_i, Z_i)^{\rm T} $ (1)
$ \pmb p_i=(f_i/Z_i){\pmb P_i}=(x_i, y_i, f_i)^{\rm T} $ (2)

图 3所示的图像 $I$ $J$ 之间的本质矩阵为$(E_{ij}=[\pmb t]_{\times} R_{ij})~(其中(R_{ij}=R_i^{\rm T} R_j, \pmb t_{ij}=R_i^{\rm T}(\pmb t_j-{\pmb t_i})) $分别表示图像 ${i, j}$ 之间的相对旋转和相对平移, $ ([\pmb t]_{\times}) $为平移向量的叉乘).

图 3 图像之间的极线关系 Figure 3 Epipolar relationship between image pairs

图像 $I$ , $J$ 的极线约束可表示为

$ \begin{equation} {\pmb p}_i^{\rm T} E_{ij} {\pmb p_j} = 0 \end{equation} $ (3a)

$ (p_i'=(x_i, y_i)^{\rm T})$表示三维点投影到第 $i$ 个像平面上的平面坐标, 则(3a)可改写为

$ \begin{equation} \left[ \begin{array}{c} \dfrac{x_i}{f_i} \\[4mm] \dfrac{y_i}{f_i} \\ 1 \\ \end{array} \right] E_{ij} \left[ \begin{array}{c} \dfrac{x_j}{f_j} \\[4mm] \dfrac{y_j}{f_j} \\ 1 \\ \end{array} \right] = 0 \end{equation} $ (3b)
$ \begin{equation} \left[ \begin{array}{c} \dfrac{p_i'}{f_i} \\ 1 \\ \end{array} \right] E_{ij} \left[ \begin{array}{c} \dfrac{p_j'}{f_j} \\ 1 \\ \end{array} \right] = 0 \end{equation} $ (3c)

为了获取由相对平移方向的线性估计模型, 首先对原有的极线约束方程(3)进行变形, 获得表达式(4):

$ \begin{align} {\pmb p}_i^{\rm T} E_{ij} {\pmb p_j} &= {\pmb p}_i^{\rm T} [R_i^{\rm T}({\pmb t_i}-{\pmb t_j})]_{\times}R_i^{\rm T} R_j {\pmb p_j} =\nonumber \\ & {\pmb p}_i^{\rm T} R_i^{\rm T}(({\pmb t_i}-{\pmb t_j})_ {\times}R_j {\pmb p_j})= \nonumber \\ &(R_i{\pmb p_i}{\times}R_j {\pmb p_j})^{\rm T}({\pmb t_i}- {\pmb t_j})=0 \nonumber \\ \Longleftrightarrow \quad& \pmb v_{ij}^{\rm T}(\pmb t_i-{\pmb t_j}) = 0 \nonumber \\ \Longleftrightarrow \quad& \pmb v_{ij}^{\rm T} {\pmb \gamma}_{ij} = 0 \nonumber \\ {\pmb v}_{ij} &= (R_i {\pmb p_i}{\times}R_j {\pmb p_j}) =\nonumber \\ & \left[ \left( R_i \left[ \begin{array}{c} \dfrac{p_i'}{f_i} \\[2mm] 1 \end{array} \right] \right) \right] {\times} \left[ \left( R_j \left[ \begin{array}{c} \dfrac{p_j'}{f_j} \\[2mm] 1 \end{array} \right] \right) \right] \end{align} $ (4)

将式(4)推广到 $k$ 对匹配点, 得到下面的基于 $L_1$ 范数的线性估计模型:

$ \begin{align} \mathop {\rm minimize}_{\{\pmb \gamma_{ij}'\}_{ij}}&\sum\limits_{k=1}^{m_{ij}}\delta \nonumber \\ \mbox{subject to} \quad& |(\pmb \gamma_{ij}')^{\rm T} {\pmb v}_{ij}^k| \le \delta, \quad \forall i, j \nonumber \\ &||\pmb \gamma_{ij}'||=1, \quad \forall i, j \end{align} $ (5)

对式(5)进行迭代加权最小二乘法(Iteratively reweighted least squares, IRLS) [24]处理, 可以获取符号不确定的相对平移方向估计值 $\pmb \gamma_{ij}'$ .

具体步骤为:

步骤1. 设置初始权值和最大迭代次数;

步骤2. 计算约束矩阵并对矩阵分解获取相对平移方向估计值;

步骤3. 计算残差并更新权值;

步骤4. 判断估计结果是否收敛, 若不收敛, 重复步骤2和3, 否则结束运算.

在迭代过程中参与计算的每对匹配点的权值是根据每次迭代结果的残差分配的, 对残差贡献大的匹配点分配更小的权值, 反之分配更大的权值.

得到 $\pmb \gamma_{ij}'$ 后, 离最终估计 $\pmb \gamma_{ij}$ 还差一步符号转换:

$ \begin{equation} \pmb \gamma_{ij} = s{\pmb \gamma'_{ij}}, \quad s \in \{-1, +1\} \end{equation} $ (6)

式中, $s$ 表示平移方向的符号.重建出来的稀疏场景点必定位于相机平面前方这一既定事实是确定符号 $s$ 正负的主要依据[25].具体方法是:首先令 $s=+1$ , 然后计算出位于相机平面前方的场景点个数 $n_f$ , 如果大于所有场景点总数的 $1/k$ , 则保持符号不变(一般 $k$ 取2, 因为场景点有一半位于相机平面前面足以说明 $s$ 符号为正); 否则, 令 $s=-1$ .最终获得相对平移方向估计值 $\pmb \gamma_{ij}$ , 从而实现相对平移方向的鲁棒性估计.

相对于传统的平移方向估计算法, 本文不是直接利用本质矩阵分解得到 $R_{ij}, \pmb t_{ij}$ 进而获得相对平移方向值, 而是通过由极线约束变形构造一个线性的优化方程, 直接确立相对平移方向同匹配点对 $\pmb p_i, \pmb p_j$ 及绝对旋转 $R_i$ 之间的关系.这种估计模型基于 $L_1$ 范数的, 并且通过计算过程中更新匹配点对对应的权值可以有效降低误匹配对估计结果的影响, 减少相对平移方向估计结果中的外点, 提高后续相机全局位置估计的精度.

3 基于凸优化的鲁棒性全局位置估计

在全局式SfM算法流程中, 相机全局位置估计是指在全局坐标系下对参与重建的图像所对应的相机位置的估计.为了表述方便, 文中使用对极图 $G_l=(V_l, E_l)$ 来表示图像拍摄时相机之间的相对关系(相机全局位置和相对平移方向), 其中顶点 $V_l=\{1, 2, 3, \cdots, n\}$ 对应于 $n$ 个相机在全局坐标系空间位置 $\pmb t_i~(1\le i \le n)$ , 顶点之间连线的方向 $(i, j) \in E_l$ 对应的是全局坐标系下相机之间的相对平移方向 $\pmb \gamma_{ij}$ .

为了获取高效、鲁棒、全局最优的相机全局位置估计值, 本文构建了一个基于 $L_1$ 约束的凸优化线性估计模型(第3.2节).由于当需要估计的框架本身不唯一时(即出现如图 5(a)所示情形时, 详见第3.1节), 直接利用凸优化估计模型进行解算并不能保证获取唯一解.因此在将已知量代入模型进行未知参数解算之前还需要做一些预处理步骤.预处理步骤的目的在于将估计问题转换为一个适定性问题.文献[10, 26]中提出的获取三视图间三焦张量的方法可以有效解决这种问题, 但三焦张量的计算会消耗额外的时间并增加了后续位置估计模型的复杂度.本文基于平行刚体理论提出一种新的高效预处理方法.

图 4 相对平移方向和相机全局位置 Figure 4 Relative translation direction and cameras' global position
图 5 平行刚体示例 Figure 5 An example of parallel rigid
3.1 平行刚体

平行刚体理论在基于方位测量的(Bearing-based)传感器网络定位、机器人自主定位与导航、目标跟踪等方面有很多应用[27-32].本文首次将其引入到全局式SfM算法的相机全局位置估计研究中.

针对相机位置的全局估计, 首先考虑以下问题:已知相机之间的相对平移方向 $\pmb \gamma_{ij}$ 时, 能否唯一确定所有相机的位置 $\pmb t_i~(1 \le i \le n)$ (这里的唯一性不考虑全局的平移和缩放, 如图 2所示)?当对极图有什么样的属性时, 可唯一确定相机位置?如果不能唯一确定相机的位置, 是否可以从对极图中提取出子图, 从而唯一确定所有相机的全局位置?

上述几个问题符合平行刚体理论所研究的内容:已知图 $G$ 中所有边的方向, 当 $G$ 是平行刚体时, 能唯一确定图中所有顶点的位置; 当 $G$ 不是平行刚体时, 通过提取 $G$ 中的最大平行刚体成分可唯一确定图中所有顶点的位置.因此可将相机全局位置估计同平行刚体理论有效结合起来, 从而将相机全局位置估计转换为一个适定性问题.

以下为确定平行刚体存在性的定理:

定理1 [28].已知图 $G=(V, E)$ , 令 $(d-1)E$ 表示 $E$ $d-1$ 个备份, 当且仅当非空集合 $D\subseteq(d-1)E$ 且存在 $|D|=d|V|-(d-1)$ $G$ 为平行刚体. $D$ 的子集 $D'\subseteq D$ 满足下面的不等式:

$ \begin{equation} |D'| \le {d|V(D')|-(d+1)} \end{equation} $ (7)

式中, $V(D')$ 表示 $D'$ 所对应的顶点, $|*|$ 表示*中元素的个数.

对于定理1, 可以结合图 5进行理解.当 $d=2$ 时, 即在二维情形下, 此时 $D=E$ , 要满足定理1中的情形必须存在 $|D|=2|V|-3$ . 图 5(a) $|D|=6, |V|=5$ 显然不满足条件. 图 5(b)为从图 5(a)中提取出的子图 $|D|=3, |V|=3$ , 满足该条件. 图 5(c) $|D|=7, |V|=5$ , 满足该条件.所以当 $d=2$ 时, 图 5(b)(c)中描述的情形为平行刚体, 而图 5(a)则为非平行刚体.

平行刚体与解唯一的一致性还可以从图 5中直观的理解, 在图 5(a)中已知图 $G=(V, E)$ 中的顶点和顶点间连线的方向.对于图 5(a)中的情形可以获得多个独立的解, 例如 $\{t_1, t_2, t_3, t_4, t_5\}$ $\{t_1, t_2, t_3, t_4', t_5'\}$ .而图 5(b)(c)中的情形只会得到唯一的解.

由定理1衍生了一系列检测平行刚体的算法.文献[33]介绍了一种时间复杂度为O $(n^2)$ 的Pebble game变种算法检测平行刚体是否存在以及文献[34]介绍了一种算法用来从非平行刚体结构中提取出其中的最大平行刚体成分, 其时间复杂度为多项式时间.

当无外点的情况下(相机相对平移方向观测值均为真值), 结合平行刚体理论和基本的位置估计模型(式(8))可准确地恢复相机的全局位置.但对于真实的图像, 外点是普遍存在的, 尤其是面对大规模、无序图像时(因为误匹配难以避免).平行刚体理论应用于相机全局位置估计的意义在于将原本的估计问题转换成一个适定性问题, 消除由于待估计框架不唯一带来的估计结果歧义性, 保证相机位置估计的稳定性.

将平行刚体理论引入到全局式SfM中, 进行相机全局位置估计基本思路是首先确定对极图 $G_l=(V_l, E_l)$ 是否平行刚体, 如果是则直接进行鲁棒性的相机全局位置估计(第3.2节), 否则提取出其中最大平行刚体成分(利用文献[32]中算法), 然后再进行相机全局位置的鲁棒性估计.

3.2 鲁棒性全局位置估计

当相对平移方向估计结果中无外点时, 可以从对极图 $G_l=(V_l, E_l)$ 中提取出最大的平行刚体, 然后通过下式精确恢复出相机全局位置.

$ \begin{equation} \pmb \gamma_{ij}=\frac{\pmb t_i-{\pmb t_j}}{||\pmb t_i-{\pmb t_j}||} \end{equation} $ (8)

式中, $\pmb \gamma_{ij}$ 表示相机 $i, j$ 之间的相对平移方向(相机相对平移方向的鲁棒性估计见第2节), $\pmb t_i, \pmb t_j$ 表示相机 $i, j$ 的全局位置.

然而由于外点的普遍存在, 直接利用式(8)进行全局位置估计的结果通常与真值有很大偏差.如何在外点存在的情况下鲁棒的估计相机全局位置是个很重要的问题, 这也是本节探讨的内容.

当已知相机 $i, j$ 之间的相对平移方向观测值为 $\pmb \gamma_{ij}$ , 同真值之间的差值为 $\pmb \varepsilon_{ij}$ 时, 对于每个 $\pmb \gamma_{ij}$ 满足如下等式:

$ \begin{equation} \pmb \gamma_{ij}=\frac{\pmb t_i-{\pmb t_j}}{||\pmb t_i-{\pmb t_j}||}+\pmb \varepsilon_{ij} \end{equation} $ (9)

为了误差 $\varepsilon_{ij}$ 存在的情况下能快速鲁棒估计出相机的全局位置, 将式(9)改写为如下线性形式:

$ \begin{equation} \pmb \varepsilon_{ij}'=\pmb t_i-{\pmb t_j}-{\pmb \lambda_{ij}}{\pmb \gamma_{ij}} \end{equation} $ (10)

式中, 为了简化计算、提高计算效率, 令 $\lambda_{ij}=||\pmb t_i-{\pmb t_j}||$ .且使用相机间的距离值误差 $\pmb \varepsilon_{ij}'$ 来代替相对平移方向误差 $\pmb \varepsilon_{ij}$ , 最终组成一个由$ (\pmb \varepsilon_{ij}', \pmb t_i, \pmb t_j, \lambda_{ij}, \pmb \gamma_{ij})$构成的线性方程.相机之间的距离误差 $\pmb \varepsilon_{ij}'$ 和相对平移方向误差 $\pmb \varepsilon_{ij}$ 成线性关系( ${\pmb \varepsilon_{ij}'}=-{\pmb \lambda_{ij}}{\pmb \varepsilon_{ij}}$ ), 所以最小化相对平移方向误差 $\pmb \varepsilon_{ij}$ 的问题可以转化为最小化相机间距离误差 $\pmb \varepsilon_{ij}'$ 的问题.为了提高对外点的鲁棒性, 本文基于 $L_1$ 范数约束对距离误差 $\pmb \varepsilon_{ij}'$ 进行最小化.

由此可以获得如下的基于约束的估计模型:

$ \begin{align} &\mathop{\rm minimize}_{\{\pmb t_i\}_i, \{\lambda_{ij}\}_{ij}, \{\pmb \varepsilon_{ij}'\}_{ij}} \sum\limits_{(i, j) \in E_l}{\pmb \varepsilon_{ij}'} \nonumber \\ &\mbox{subject to}\ ||\pmb t_i-{\pmb t_j}-\lambda_{ij}{\pmb \gamma_{ij}}|| \le \pmb \varepsilon_{ij}' \nonumber \\ &\qquad\sum\limits_{i \in V_l}{\pmb t_i}=0 \nonumber \\ &\qquad\lambda_{ij}=||\pmb t_i-{\pmb t_j}|| \ge c, \forall (i, j) \in E_l \end{align} $ (11)

式(11)中估计模型是非线性且非凸的, 直接使用该模型进行相机全局位置的估计会存在解算复杂度高, 求解困难的问题.因此舍弃非凸约束 $\lambda_{ij}=||\pmb t_i-{\pmb t_j}||$ , 将 $\lambda_{ij}$ 作为独立的参数求解(此时 $\lambda_{ij}$ 表示的是对应于每个相机对的未知尺度因子).这种处理策略一方面可以将式(11)中估计模型线性化(这种线性化方法广泛应用于相机全局位置估计过程中, 文献[10, 17, 20]等均采取这种策略直接将 $\lambda_{ij}$ 作为一个独立的未知尺度因子, 其中文献[20]首次提出利用线性估计模型进行相机全局位置的估计).另一方面可以将式(11)中的估计模型转换为具有凸属性的估计模型.

通过上述舍弃非凸约束的策略将式(11)中估计模型从非凸非线性估计模型转换为具有凸属性的线性估计模型, 能够提高估计过程的效率和精度.如果使用原始的非凸非线性模型即使在数据量很小的情况下也会出现计算复杂度过高的问题.

在相机全局位置估计过程中有大量的多余观测(Redundance data), 根据式(11)可以列出多个同相机 $i$ 相关的约束参与迭代解算(约束个数等于与 $i$ 相关的相机对的总数, 例如相机对 $i, j$ $i, j+k$ $i, j+n\cdot$ , 理论上最多可列 $N-1$ 个约束, 其中 $N$ 为参与估计的相机总数).在这种情况下将 $\lambda_{ij}$ 作为一个独立的未知尺度因子可以降低参数之间的耦合性, 简化计算, 并不会降低最终估计的精度.

最终获得如下基于 $L_1$ 约束的凸优化线性估计模型:

$ \begin{align} &\mathop{\rm minimize}_{\{\pmb t_i\}_i, \{\lambda_{ij}\}_{ij}, \{\pmb \varepsilon_{ij}'\}_{ij}} \sum\limits_{(i, j) \in E_l}{\pmb \varepsilon_{ij}'} \nonumber \\ &\mbox{subject to}\ ||\pmb t_i-{\pmb t_j}-\lambda_{ij}{\pmb \gamma_{ij}}|| \le \pmb \varepsilon_{ij}' \nonumber \\ &\qquad \sum\limits_{i \in V_l}{\pmb t_i}=0 \nonumber \\ &\qquad \lambda_{ij} \ge c, \forall (i, j) \in E_l \end{align} $ (12)

式中, 约束 $\sum_{i \in V_l}{\pmb t_i}=0$ $\lambda_{ij} \ge c$ 分别是为了消除估计结果中平移和尺度的二义性(同文献[10]一致, 本文取 $c$ 值为1).另外约束 $\lambda_{ij} \ge c$ 可以避免诸如 $\lambda \equiv 0, \pmb t \equiv 0$ 之类的平凡解.

对于上述结构的凸优化估计模型, 通常使用二次规划(Quadratic programing, QP)求解.为了提高估计结果的鲁棒性, 本文结合迭代加权最小二乘(Iteratively reweighted least squares, IRLS)对模型进行求解.

IRLS主要思想是:对目标模型进行迭代求解, 在每次迭代过程中更新观测值权重.权重的大小取决于上一次迭代结果的残差, 对残差贡献小的观测值给于更高的权重, 反之更低.当应用于本节中相机全局位置估计问题时, 如下式(13)所示.

$ \begin{equation} \mathop{\rm minimize}\limits_{\{\pmb t_i\}_i, \{\lambda_{ij}\}_{ij}, \{\gamma_{ij}\}_{ij}} \sum\limits_{(i, j) \in E_l}w_{ij}||\pmb t_i-{\pmb t_j}-\lambda_{ij}{\pmb \gamma_{ij}}|| \end{equation} $ (13)

对残差贡献小的相对平移方向值观测值给予更大的权重.这种方式可有效地提高算法对外点的鲁棒性.

主要算法步骤如下:

步骤1. 设置QP和IRLS的迭代次数以及初始权重.

步骤2. 根据式(12)计算约束矩阵 $A$ .

步骤3. 使用QP求解矩阵方程.

步骤4. 计算残差并更新权重.

步骤5. 重复步骤4和5直至达到迭代次数的上限.

同文献[10]方法相比, 本文使用了基于 $L_1$ 范数而不是 $L_{\infty}$ 范数的估计模型, 对外点更具有鲁棒性.同文献[23] 1DSfM方法相比, 本文使用了线性的而不是非线性的优化方程, 在计算效率方面更高, 并且本文使用的是凸优化估计模型能够保证收敛到全局最优, 而1DSfM只能收敛到局部最优.

4 实验

实验数据为文献[23]提供的8组公开图像数据集.这8组数据都具有大规模、无序的特点.为了验证本文算法的效率和精度, 对这8组数据进行处理并同文献[20, 23]的实验结果进行比较.下面给出实验结果及结果分析

4.1 实验结果

实验结果分为两部分, 第一部分是全局位置估计精度的比较(见表 1), 第二部分是处理时间的比较(见表 2), 另外给出了利用本文算法对8组数据处理获取的场景三维稀疏结构(图 11).

表 1 本文方法同1DSfM、文献[20]处理精度比较 Table 1 Comparison of accuracy: our method、1DSfM and [20]
表 2 本文方法同1DSfM、文献[20]、Bundler处理时间比较 Table 2 Comparison of efficiency: our method、1DSfM、Bundler and [20]
图 11 利用本文算法对8组公开数据集处理获取的场景三维稀疏结构 Figure 11 The experimental result with 8 groups of datasets based on our method

为了控制其他因素的影响, 保证实验所采用的特征提取与匹配、相对运动估计等算法同文献[20, 23]保持一致, 采用文献[23]给出的相对运动估计结果作为本文算法的输入.为了保证公平性, 程序只使用了单线程进行处理.

表 1 $\overline{x}, \widetilde{x}$ 分别表示估计距离同参考距离之间差值的平均值和中位值(单位近似为米)其中1DSfM和文献[20]的实验结果引用自文献[23] (1DSfM结果中未给出BA前的平均值, 文献[20]结果中只给出了BA后的中位数).

增量式Bundler方法[3]作为处理大规模无序图像的经典算法, 在相机位置估计精度方面比较可靠.因此, 同文献[23]一样, 为了对实验结果的精度进行量化比较, 使用Bundler的处理结果作为参考(该结果从文献[23]给出的网站获取), 分别比较估计值同参考值的差值的平均值、中位数(采用了一种基于RANSAC (Random sample consensus)的绝对定位方法, 将实验结果对齐到参考框架中从而获取平均值和中位数).选择平均值和中位数作为比较量的原因在于:平均值反映的是全体数据的平均水平, 平均值的大小易受极端值得影响, 因此可以反映算法对外点的鲁棒性; 中位数反映的是全体数据的集中趋势, 不受少数极端值影响, 中位数的大小可以反映估计结果的普遍精度.

表 2中比较了本文方法同1DSFM、文献[20]、及传统的增量式Bundler方法的处理时间(单位为秒).表中 $T_R$ 表示绝对旋转所需的时间, $T_O$ 表示对相对平移方向进行优化所需时间, $T_S$ 表示全局位置估计所需的时间, $T_{BA}$ 表示进行BA所需的时间, $\sum$ 表示各个阶段时间的总和.

4.2 结果分析

数据处理精度方面, 从表 1可以看出, 本文方法在BA优化之前的平均值比1DSfM方法经过BA优化后平均值更低, 在本文方法经过BA优化后两者的差距进一步扩大.其中Vienna和Alamo两个数据的中位数更是相差了2E3、1E7倍(图 7, 从图 6中也可以看出对于数据Alamo使用本文方法只有少量点偏离参考值相对较大).这表明了本文实验结果中具有更少的偏离参考值的极端值, 充分反映了本文方法对外点的鲁棒性.并且同1DSfM相比, 在BA优化前后本文方法都具有与前者近似或者更低的中位数(图 9图 10).同文献[20]相比, 本文方法显然具有更低的中位数(图 8).这表明了本文方法精度方面的可靠性.

图 6 相机全局位置散点图(数据Vienna, 相机个数821) Figure 6 Global location of cameras represented by scatter diagram
图 7 BA优化后本文同1DSfM实验结果平均值比较图 Figure 7 Comparison result of mean between our method and 1DSfM after BA
图 8 BA优化后本文同文献[20]实验结果中位数比较图 Figure 8 Comparison result of median between our method and [20] after BA
图 9 BA优化前本文同1DSfM实验结果中位数比较图 Figure 9 Comparison result of median between our method and 1DSfM before BA
图 10 BA优化后本文同1DSfM实验结果中位数比较图 Figure 10 Comparison result of median between our method and 1DSfM after BA

数据处理的效率方面, 从表 2可以看出, 相比较传统的增量式方法, 本文方法有6 $\sim$ 10倍加速; 相比1DSfM, 具有接近于2倍的提速.对于大部分实验数据, 本文处理速度不低于文献[20]方法.因此, 本文的方法在处理效率上优于或者接近于现有的典型方法.

综上, 实验结果表明:本文方法对外点具有更好的鲁棒性, 并且处理结果的普遍精度和处理效率都不低于典型方法, 通过本文方法可以高效、鲁棒地进行相机全局位置估计.

5 结论

针对现有全局式SfM算法流程中的相机全局位置估计部分存在的不足, 本文提出一种改进的相机全局位置估计方法.本文工作主要体现在两个方面: 1)结合极线约束, 提出一种新的对误匹配鲁棒相对平移方向估计算法; 2)引入平行刚体理论将全局位置估计转换为一个适定性问题并提出一种基于约束的凸优化线性估计模型.实验结果表明本文改进的相机全局位置估计方法对外点具有更好的鲁棒性, 并且在计算效率和估计结果的普遍精度方面也具有一定优势.

在实验过程中我们发现鲁棒性BA对实验结果也有很重要的影响, 所以下一步的研究工作是尝试将本文的估计方法和鲁棒性BA结合以获取更好的估计结果.

参考文献
1
Triggs B, McLauchlan P, Hartley R I, Fitzgibbon A W. Bundle adjustment-a modern synthesis. Vision Algorithms: Theory and Practice: Lecture Notes in Computer Science. Berlin, Heidelberg: Springer, 1999. 298-372 http://rd.springer.com/book/10.1007/3-540-44480-7
2
Zhang Z Y, Shan Y. Incremental Motion Estimation Through Local Bundle Adjustment. Technical Report MSR-TR-01-54, Microsoft Research, Redmond, WA, 2001.
3
Snavely N, Seitz S M, Szeliski R. Photo tourism:exploring photo collections in 3D. ACM Transactions on Graphics, 2006, 25(3): 835-846. DOI:10.1145/1141911
4
Snavely N, Seitz S M, Szeliski R. Skeletal graphs for efficient structure from motion. In: Proceedings of the 2008 IEEE Conference on Computer Vision and Pattern Recognition. Anchorage, USA: IEEE, 2008. 1-8 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=4587678
5
Havlena M, Torii A, Knopp J, Pajdla T. Randomized structure from motion based on atomic 3D models from camera triplets. In: Proceedings of the 2009 IEEE Conference on Computer Vision and Pattern Recognition. Miami, USA: IEEE, 2009. 2874-2881 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5206677
6
Furukawa Y, Curless B, Seitz S M, Szeliski R. Towards internet-scale multi-view stereo. In: Proceedings of the 2010 IEEE Conference on Computer Vision and Pattern Recognition. San Francisco, USA: IEEE, 2010. 1434-1441 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5539802
7
Sinha S N, Steedly D, Szeliski R. A multi-stage linear approach to structure from motion. In: Proceedings of the 11th European Conference on Trends and Topics in Computer Vision. Berlin, Heidelberg: Springer, 2010. 267-281 https://link.springer.com/chapter/10.1007%2F978-3-642-35740-4_21
8
Kneip L, Chli M, Siegwart R Y. Robust real-time visual odometry with a single camera and an IMU. In: Proceedings of the 22nd British Machine Vision Conference. Scotland, UK: BMVA, 2011. http://dx.doi.org/10.3929/ethz-a-010025746
9
Tomasi C, Kanade T. Shape and motion from image streams under orthography:a factorization method. International Journal of Computer Vision, 1992, 9(2): 137-154. DOI:10.1007/BF00129684
10
Moulon P, Monasse P, Marlet R. Global fusion of relative motions for robust, accurate and scalable structure from motion. In: Proceedings of the 2013 IEEE International Conference on Computer Vision. Sydney, AU: IEEE, 2013. 3248-3255 http://ieeexplore.ieee.org/document/6751515/
11
Govindu V M. Lie-algebraic averaging for globally consistent motion estimation. In: Proceedings of the 2004 IEEE Conference on Computer Vision and Pattern Recognition. Washington, USA: IEEE, 2004. 684-691 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=1315098
12
Martinec D, Pajdla T. Robust rotation and translation estimation in multiview reconstruction. In: Proceedings of the 2007 IEEE Conference on Computer Vision and Pattern Recognition. Minnesota, USA: IEEE, 2007. 1-8 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=4270140
13
Hartley R, Aftab K, Trumpf J. L1 rotation averaging using the Weiszfeld algorithm. In: Proceedings of the 2011 IEEE Conference on Computer Vision and Pattern Recognition. Providence, USA: IEEE, 2011. 3041-3048 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5995745
14
Fredriksson J, Olsson C. Simultaneous multiple rotation averaging using lagrangian duality. In: Proceedings of the 11th Asian Conference on Computer Vision. Berlin, Heidelberg: Springer, 2012. 245-258 https://link.springer.com/chapter/10.1007%2F978-3-642-37431-9_19
15
Chatterjee A, Govindu V M. Efficient and robust large-scale rotation averaging. In: Proceedings of the 2013 IEEE International Conference on Computer Vision. Sydney, AU: IEEE, 2013. 521-528 http://ieeexplore.ieee.org/document/6751174/
16
Brand M, Antone M, Teller S. Spectral solution of large-scale extrinsic camera calibration as a graph embedding problem. In: Proceedings of the 8th European Conference on Computer Vision. Berlin, Heidelberg: Springer, 2004. 262-273 http://www.springerlink.com/content/0bx5vqf0688nxcw3
17
Arie-Nachimson M, Kovalsky S Z, Kemelmacher-Shlizerman I, Singer A, Basri R. Global motion estimation from point matches. In: Proceedings of the 2nd International Conference on 3D Imaging, Modeling, Processing, Visualization and Transmission. Zurich, Switzerland: IEEE, 2012. 81-88 http://ieeexplore.ieee.org/document/6374980/
18
Sim K, Hartley R. Recovering camera motion using L minimization. In: Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. New York, NY, USA: IEEE, 2006. 1230-1237
19
Kahl F, Hartley R. Multiple-view geometry under the L-norm. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2008, 30(9): 1603-1617. DOI:10.1109/TPAMI.2007.70824
20
Govindu V M. Combining two-view constraints for motion estimation. In: Proceedings of the 2001 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. Kauai, HI, USA: IEEE, 2001. 218-225 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=990963
21
Tron R, Vidal R. Distributed image-based 3-D localization of camera sensor networks. In: Proceedings of the 48th IEEE Conference on Decision and Control. Shanghai, China: IEEE, 2009. 901-908 http://ieeexplore.ieee.org/xpls/icp.jsp?arnumber=5400405
22
Li H D. Multi-view structure computation without explicitly estimating motion. In: Proceedings of the 2010 IEEE Conference on Computer Vision and Pattern Recognition. San Francisco, CA: IEEE, 2010. 2777-2784
23
Wilson K, Snavely N. Robust global translations with 1DSfM. In: Proceedings of the 13th European Conference on Computer Vision. Berlin, Heidelberg: Springer, 2014. 61-75 https://link.springer.com/chapter/10.1007%2F978-3-319-10578-9_5
24
Daubechies I, Devore R, Fornasier M, Güntürk C S. Iteratively reweighted least squares minimization for sparse recovery. Communications on Pure and Applied Mathematics, 2010, 63(1): 1-38. DOI:10.1002/cpa.v63:1
25
Andrew A M. Multiple view geometry in computer vision. Kybernetes, 2001, 30(9-10): 1333-1341.
26
Jiang N J, Cui Z P, Tan P. A global linear method for camera pose registration. In: Proceedings of the 2013 IEEE Conference on Computer Vision. Sydney, NSW: IEEE, 2013. 481-488 http://ieeexplore.ieee.org/document/6751169/
27
Whiteley W. A matroid on hypergraphs, with applications in scene analysis and geometry. Discrete & Computational Geometry, 1989, 4(1): 75-95.
28
Whiteley W. Parallel redrawing of configurations in 3-space. 1986.
29
Servatius B, Whiteley W. Constraining plane configurations in computer-aided design:combinatorics of directions and lengths. SIAM Journal on Discrete Mathematics, 1999, 12(1): 136-153. DOI:10.1137/S0895480196307342
30
Eren T, Whiteley W, Belhumeur P N. Using angle of arrival (bearing) information in network localization. In: Proceedings of the 45th IEEE Conference on Decision and Control. San Diego, CA: IEEE, 2006. 4676-4681 http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4177903
31
Eren T, Whiteley W, Morse A S, Belhumeur P N, Anderson B D O. Sensor and network topologies of formations with direction, bearing, and angle information between agents. In: Proceedings of the 42nd IEEE Conference on Decision and Control. Maui, HI: IEEE, 2003. 3064-3069 http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1273093
32
Jackson B, Jordán T. Graph theoretic techniques in the analysis of uniquely localizable sensor networks. 2009.
33
Jacobs D J, Hendrickson B. An algorithm for two-dimensional rigidity percolation:the pebble game. Journal of Computational Physics, 1997, 137(2): 346-365. DOI:10.1006/jcph.1997.5809
34
Kennedy R, Daniilidis K, Naroditsky O, Taylor C J. Identifying maximal rigid components in bearing-based localization. In: Proceedings of the 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems. Vilamoura-Algarve, Portugal: IEEE, 2012. 194-201 http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6386132