《GAMES203：三维重建和理解》2 配准（registration）

RhinoFreak · 发表于 2022-2-2 10:03

目录及序言：《GAMES203：三维重建和理解》0 目录及序言 - 知乎 (zhihu.com)
<hr/>1 配准问题的分类

配准（registration）的目标是对准那些扫描得到的图形数据片段，推断出相对的变换矩阵，把各个点云定义在局部坐标系下的数据合并到全局坐标系下，得到完整的三维模型点云数据。
可从三个方面粗略地分类配准问题：

完全重叠（full overlap）或部分重叠（partial overlap）：
- 完全重叠：待配准的数据一个是另一个的拷贝，或者被涵盖于其中。例如，扫描仿件数据，与原件数据配准，评估仿制精度。
- 部分重叠：例如，从不同视角扫描物理对象进行三维重建，得到的那些数据片段应该互相只有部分重叠。这种情况更常见。
全局（global）或局部（local）：
- 局部：待配准数据在初始时已经基本上对准了，位姿（pose）相近。例如，扫描同一物理对象时，前后两帧的数据存在连续性，它们的形态理应不会相差太大。
- 全局：点云初始位姿任意。例如，相机在一个位置扫描得到几帧数据，然后转移到另一个地方去扫描，那么这两批数据之间就有可能位姿差别很大了。
成对（pairwise）或多对（multiple）：
- 成对：配准两片点云；
- 多对：同时配准多片点云；三维重建显然需要配准多片点云。

2 ICP 算法成对配准点云

假设点云的区域完全重叠，是刚性形状（rigid shapes），能通过刚体变换相联系，初始时位姿相近，基本对准了，此时可使用 ICP 算法（iterative closest point，迭代最近点）进行配准，它以最优化衡量点云之间对准程度的函数为目标，采用迭代的方法，使得点云契合程度不断趋于局部最优。
虽然配准时一般不知道点云的对应关系，但是在输入时要求数据基本对准了，所以可以假设当前两片点云之间距离最近的点是对应的。

刚体变换（rigid transformation）是指仅对目标进行平移、旋转，而不改变目标形状的变换方式；

ICP 算法示意

给定点云 $p_i$ 和 $q_i$ ，构建误差函数 $E=\underset{i}{\sum}\left(R \cdot p_i+t-q_i\right)^2$ ，ICP 算法通过迭代寻找一个旋转变换和平移变换 $T$ ，极小化误差函数的值。当两片点云的位姿足够接近，即残差（residual）平方和小于某个阈值时，可认为算法收敛，停止迭代。

作为局部算法，依赖好的初始位姿；
每次迭代分为两步：寻找最近点、寻找刚体变换；
计算点的最近邻：
- 现在有很多优秀的软件库，调用相关算法能方便地得到两片点云之间的最近点对，以前则手动构建一些层次化的数据结构（例如 k 维树）以计算点的最近邻；
- 有些算法不需要寻找最近邻，例如有的算法直接构建了平方距离场（squared distance field）；
一般不会让点云中的所有点都参与配准计算，而是从中抽样一些点，于是，抽样的方式也会影响配准的效果。比起均匀抽样，更多地抽样那些制约潜在不稳定变换的样本点，配准的效果会更好，参见《Geometrically Stable Sampling for the ICP Algorithm》。

2.1 从最优化的视角看待 ICP 算法

配准问题可视为一个最小二乘（least squares）问题，计算刚体变换以极小化 $F(\alpha)=\underset{i}{\sum}d^2\left(\alpha\left(x_i^0\right),\,\Phi\right)$ 的值，其中 $d^2\left(\alpha\left(x_i^0\right),\,\Phi\right)$ 是 $\alpha\left(x_i^0\right)$ 与 $\Phi$ 之间的平方距离。

由于刚体变换中的旋转分量的行列之间满足正交性质，所以该问题是一个受限非线性最优化问题（constrained nonlinear optimization problem）；

ICP 算法先优化刚体变换，再优化最近点对，交替地进行极小化，一直在减小目标函数的值，是线性收敛的，一般可以收敛到局部最优；

需要注意的是，一直减小目标函数的值并不保证可以收敛到局部最优。例如，对于 $f(x)=x^2$ ，迭代形式 $x_i=3+\frac{1}{i}$ 可以一直减小目标函数的值，但是不能收敛到 $x=0$ 。

一般采用高斯-牛顿算法求解非线性最小二乘问题；

配准时使用高斯-牛顿算法迭代，等同于极小化点与点所在面之间的距离，而不再是点与点之间的距离；
- 极小化点与面之间距离的算法变体参见《Object modeling by registration of multiple range images》
- 点与面之间距离的计算可参考《Geometry of the Squared Distance Function to Curves and Surfaces》；
当误差函数的初始值较小，即两片点云在初始状态已经粗略匹配了时，ICP 算法使用高斯-牛顿算法进行迭代的收敛速度会非常快，接近于二阶收敛，而当误差函数的初始值较大时，一般也会比基础的 ICP 算法收敛速度略微好一点，但基础的 ICP 算法会更稳健（robust）；
当误差函数的初始值较大时，一般需要对迭代中黑塞矩阵等的计算做一些特别的处理，例如 BFGS 算法。

2.2 补充：牛顿法（Newton method）与高斯-牛顿算法（Gauss–Newton algorithm）

牛顿法（Newton method）是指在搜索区间上逐次构造新的、与所寻求函数相应的二次函数，并用一系列二次函数的极小点逐步逼近原寻求函数极小点的一种方法。

一维情况的牛顿法示例：
设实值函数 $f(x)$ 在区间 $\left[a,\,b\right]$ 只有一个严格局部极小值，在所讨论点处 $f(x_k)$ 、和 $f''(x_k)$ 都存在，则在的领域有二阶泰勒展开： $f(x_k+t)=f(x_k)+f'(x_k)t+\frac{1}{2}f''(x_k)t^2+o\left(t^2\right) \\$
于是，构造二次函数：
$\varphi(t)=f(x_k)+f'(x_k)t+\frac{1}{2}f''(x_k)t^2 \\$
通过使 $\frac{\mathrm{d}\varphi(t)}{\mathrm{d}t}=0$ ，解得新函数 $\varphi(t)$ 的极小点 $t=-\frac{f'(x_k)}{f''(x_k)}$ ，
那么，便可以用新函数的极小点作为原寻求函数极小点的近似值： $x_{k+1}=x_k+t=x_k-\left[f''(x_k)\right]^{-1}f'(x_k)\\$
如果这个近似值不满足预先提出的精度要求，则在 $x=x_{k+1}$ 处再构造一个二次函数，重复该步骤。
将牛顿法推广到高维，则此时是梯度（gradient）、是雅可比矩阵（Jacobian matrix）， $\left[f''(x_k)\right]^{-1}$ 是黑塞矩阵（Hessian matrix）的逆矩阵或穆尔-彭罗斯伪逆（Moore-Penrose pseudoinverse）。

高斯-牛顿算法（Gauss–Newton algorithm）在牛顿法的基础上修改而来，只能用来求解最小二乘问题，因为它基于最小二乘的假设，对二阶导数（黑塞矩阵）做了近似。
给定 $m$ 个维函数 $r=\left(r_1,\,\cdots,\,r_m\right)$ 和个变量 $\beta=\left(\beta_1,\,\cdots,\,\beta_n\right)$ ，其中 $r_i:\mathbb{R}^n\to\mathbb{R},,i=1,\cdots,m$ ， $m \ge n$ ，求解最小二乘问题，极小化目标函数 $S(\beta)=\overset{m}{\underset{i=1}{\sum}}r_i\left(\beta\right)^2$ 的值，牛顿法构造的迭代形式如下：

$\beta^{(k+1)}=\beta^{(k)}-H^{-1}\cdot \nabla S\left(\beta^{(k)}\right) \\$
其中， $H$ 是 $S(\beta)$ 的黑塞矩阵，于是有：

$\begin{align} &\because S=\overset{m}{\underset{i=1}{\sum}}r_i^2 \\ &\therefore \nabla S_j =2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}\\ &\therefore H_{jk} =\frac{\partial \nabla S_j}{\partial \beta_k}\\ &\quad\quad\quad =2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\cdot\frac{\partial r_i}{\partial \beta_k}+\color{Red}{{r_i\cdot \frac{\partial^2 r_i}{\partial \beta_j\,\partial \beta_k}}}\right) \quad \color{Green}{\text{高斯-牛顿法忽略了红色项}} \\ &\quad\quad\quad \approx 2\sum_{i=1}^m\frac{\partial r_i}{\partial \beta_j}\cdot\frac{\partial r_i}{\partial \beta_k} \quad\quad \color{Green}{\text{用两个雅可比矩阵的乘积去近似黑塞矩阵}}\\ &\quad\quad\quad= 2\sum_{i=1}^m J_{ij}J_{ik}\\ &\therefore \nabla S=2J_r^T r,\,H\approx2J_r^TJ_r \end{align} \\$
于是得到了高斯-牛顿算法构造的迭代形式：

$\beta^{(k+1)}=\beta^{(k)}-\left(2J_r^TJ_r\right)^{-1}\cdot 2J_r^T r\left(\beta^{(k)}\right)\\ \text{整理得：}\beta^{(k+1)}=\beta^{(k)}-\left(J_r^TJ_r\right)^{-1}J_r^T \cdot r\left(\beta^{(k)}\right) \\$

可发现 $\left(J_r^TJ_r\right)^{-1}J_r^T$ 是矩阵 $J_r$ 的左伪逆矩阵。

2.3 图形数据部分重叠时的配准

待配准图形数据只有部分区域重叠

当待配准的两片点云只有部分区域重叠时：

可以将 ICP 算法构造的误差函数中的 $l_2$ 范数（平方距离， $\rho_2(t)=t^2$ ）换成稳健范数（robust norm，例如 $\rho_{GM}(t)=\frac{t^2}{\sigma^2+t^2}$ 、 $\rho_{1}(t)=|t|$ ），然后还是用高斯-牛顿算法求解；
用得比较多的，更易于解释的方法是用 IRLS（iteratively reweighted least squares，迭代重加权最小二乘）求解。原本的 ICP 算法可看作所有最近邻点对的权重都是，而 IRLS 在每次迭代后，根据残差调整每个点对的权重，当前残差越大，即当前点与它的最近点的距离越远，则令该点对的权重 $w_i$ 越小。
也有做法是通过双向修剪（bi-directional pruning），剔除不重叠部分的那些多余的点：计算点的最近邻，再计算的最近邻，如果和距离太远则剔除，否则认为和对应。

双向修剪（bi-directional pruning）示意

ICP 算法的一些其它变体参见《Efficient Variants of the ICP Algorithm》。
3 全局匹配（global matching）

全局匹配算法不要求数据有好的起始位姿，于是就不能假设点云之间最近的点是对应的了，而是需要首先提取一些在刚体变换下保持不变的几何特征量，然后认为几何特征相同的点是对应点，在这些特征点之间建立初步的对应关系，最后根据一些刚性约束，筛选出好的一致的对应关系，拟合刚体变换。

全局匹配示意：先提取特征，再匹配特征，最后提取相对位姿

现在基于深度学习的方式的大致思路是先用神经网路提取特征，然后用相联系的神经网路进行匹配。

基于深度学习的方法示意：使用神经网路提取特征，再用另一个神经网络进行匹配

3.1 特征描述子（feature desciptor）

全局匹配算法首先根据特征描述子（feature desciptor）提取特征。
一些人工设计的特征描述子举例：

自旋图像（spin images），参见《Spin-Images: A Representation for 3-D Surface Matching》；
- 是对该点几何信息的详细描述，编码了当前点与所在曲面片（patch）内其它点的距离，以及两者间连线与当前点表面法线的相对角度，记录为附近点到当前点法线的距离和到当前点所在切平面的距离，
- 区分度大，但计算和存储的开销也较大；
积分不变量（integral invariants），参见《Integral Invariants for Robust Geometry Processing》；
- 记录了一个以当前点为球心的一定大小的球中，有多大比例的面积（体积）是在模型内部；
- 对尺度变换、噪声不敏感；

作为特征的积分不变量示意

3D SIFT，参见《A 3-Dimensional Sift Descriptor and Its Application to Action Recognition》；
Patch features，可参见《Salient Geometric Features for Partial Shape Matching and Similarity》；

如果追求稳健（robust），那么最好使用人工设计的特征描述子，而不是深度学习得到特征，因为深度学习存在泛化误差。
3.2 对应关系的一致性（consistency）

因为刚体变换不改变形状，所以对应关系联系两个特征点应该是一致的。

例如，如果对应关系 $q_1 \leftrightarrow q_2$ 和 $q_1' \leftrightarrow q_2'$ 能确定一个刚体变换，那么它们两者之间应该满足以下约束条件：

几何一致（geometric consistency），变换前后对应的距离和角度应该不变：
$\|\boldsymbol{p}(q_1)-\boldsymbol{p}(q_1')\|=\|\boldsymbol{p}(q_2)-\boldsymbol{p}(q_2')\| \\ \angle\left( \boldsymbol{n}(q_1),\,\boldsymbol{n}(q_1')\right) = \angle\left( \boldsymbol{n}(q_2),\,\boldsymbol{n}(q_2')\right)\\ \angle\left(\boldsymbol{n}(q_1),\,\boldsymbol{p}(q_1)\boldsymbol{p}(q_1')\right) = \angle\left(\boldsymbol{n}(q_2),\,\boldsymbol{p}(q_2)\boldsymbol{p}(q_2')\right)\\ \angle\left(\boldsymbol{n}(q_1'),\,\boldsymbol{p}(q_1)\boldsymbol{p}(q_1')\right) = \angle\left(\boldsymbol{n}(q_2'),\,\boldsymbol{p}(q_2)\boldsymbol{p}(q_2')\right) \\$
其中 $\boldsymbol{p}(q)$ 是特征点的位置， $\boldsymbol{n}(q)$ 是特征点所在表面的法线；
根据特征描述子提取的不变特征一致（consistency in descriptors）：
$\boldsymbol{f}\left(q_1\right)=\boldsymbol{f}\left(q_2\right) \\ \boldsymbol{f}\left(q_1'\right)=\boldsymbol{f}\left(q_2'\right) \\$
其中 $\boldsymbol{f}(q)$ 是根据特征描述子 $\boldsymbol{f}$ 在特征点处提取到的不变特征；

3.3 筛选对应关系的方法

根据特征描述子提取特征后，一般可以在两片待配准的两片点云之间建立大量的对应关系，这些对应关系有好有坏，目标是找到这些对应关系的一个子集，使子集包含的对应关系能被一个刚体变换所拟合，而寻找子集算法的核心思想是根据刚性约束条件筛选出一致的对应关系。
和二维计算机视觉中相应方法的核心思想相同。
3.3.1 随机抽样一致性（RANSAC）算法

RANSAC（random sample consensus）算法的基本原理是反复抽样数据、拟合变换，筛选出能拟合尽可能多的样本数据的变换，忽略那些异常值的影响。它的基本步骤如下：

抽样三个特征点对的位置，检查是否满足距离约束；
如果满足距离约束条件，则拟合一个刚体变换；
检查有多少其它的特征点对与当前拟合出的刚体变换一致；如果超过阈值，则停止迭代，否则回到第一步继续抽样；

如果既知道特征点的位置，又知道点所在表面的法线方向，则可略微修改 RANSAC 算法，一次抽样两个对应关系的数据，即可拟合一个刚体变换，此时除了检查距离约束之外，还需要检查角度约束。
3.3.2 霍夫变换（Hough transform）

霍夫变换（Hough transform）是基于投票原理的参数方法，它的基本原理是：因为每一个刚体变换的参数可以对应于参数空间中的一个点，于是适当地量化参数空间，参数空间内每个点的初始值为 $0$ ，对于每一个当前拟合的刚体变换，都使参数空间对应点的数值加，即向这个点代表的刚体变换投一票。在投票结束时，参数空间中数值最大的点便对应于所寻求到的刚体变换。

当霍夫变换拟合刚体变换应用于一个几何模型本身时，可以检测模型是否对称，参见《Partial and Approximate Symmetry Detection for 3D Geometry》；

把霍夫变换拟合刚体变换应用于一个几何模型本身时，可以发现参数空间中的极大值对应于把模型做对称变换

3.3.3 基于谱的匹配方法（spectral approach）

基于谱的方法的基本思想是建立一个图的邻接矩阵，它的节点代表潜在的对应关系，而边上的权重代表潜在的对应关系之间是否一致。正确的特征对应关系之间理应两两一致，能形成一个强连接的聚类簇，而不正确的对应关系应该只会恰巧与其它对应关系一致，不太可能建立强连接的聚类簇。提取最大团（clique），拟合刚体变换完成配准。

简单地说，一个给定图（graph）的团（clique）是它的一个完全子图。如果一个团不是其它任一团的真子集，则称该团为给定图的极大团（maximal clique），并称顶点最多的极大团为图的最大团（maximum clique）。

基于谱的匹配方法示意

于是，可以通过邻接矩阵的主特征向量得到主聚类簇，加上一些其它的约束条件，筛选得到正确的特征对应关系。
参见《A Spectral Technique for Correspondence Problems Using Pairwise Constraints》。
3.3.4 混合方法（hybrid method）

例如，Yang 等人将基于谱的匹配和迭代重加权最小二乘结合在一起，提出了一种混合的、稳健的几何匹配算法，参见《Extreme Relative Pose Estimation for RGB-D Scans via Scene Completion》。

迭代重加权最小二乘相当于在求解如下最优化问题：

其中， ;
- 对应关系中不正确的比例不能超过一半；
基于谱的匹配相当于在求解如下最优化问题：
- 可以容忍更多比例的不正确对应关系，但正确值和异常值分离得不彻底；
将二者结合后，相当于在求解如下最优化问题：
$\begin{align} &\underset{\{x_c\},\,R,\,t}{\text{maximize}} \sum_{c,\,c'\in C} w_\gamma(c,\,c')x_cx_{c'}\left(\delta-r_\left(R, \,t\right)(c)-r_\left(R, \,t\right)(c')\right)\\ &\text{subject to }\sum_c x_c^2=1 \end{align} \\$
其中， $w_\gamma(c,\,c')$ 代表建立特征 $c$ 与 $c^{\prime}$ 之间对应关系的一致性程度， $r_\left(R, \,t\right)(c)$ 和 $r_\left(R, \,t\right)(c')$ 是回归误差；

3.3.5 基于学习的方法（learning-based methods）

例如，Yue Wang 和 Solomon 提出了一种基于深度学习的配准算法，基本流程是首先用动态图卷积神经网络（DGCNN）将未对齐的点云数据嵌入（embed）到一个共同的空间之中，然后用一个基于注意力的模块结合指针网络（pointer network）预测两片点云的近似匹配，最后用一个奇异值分解模块提取刚体变换，得到最终的结果，参见《Deep Closest Point: Learning Representations for Point Cloud Registration》。

Deep Closest Point 算法采用的神经网络结构示意

4 同时配准多个图形数据

联合成对配准（joint pairwise registration）
ICP 算法成对配准的推广，基本思想是修改误差函数，同时极小化多片点云对应点之间的距离。
- 需要指明哪些扫描数据的哪些区域是互相重叠的；
- 当有大量扫描数据需要配准时，算法的收敛速度会比较慢；

同时配准和重建（simultaneous registration and reconstruction）：
假设输入的数据基本对准了，算法的基本思想是把空间划分成预定义分辨率的网格，在每个网格内部根据扫描数据拟合出潜在的曲面（latent surface），然后在各个网格中配准扫描数据与潜在曲面，再又优化潜在曲面，如此交替地优化，最终可配准所有的扫描数据，并且也得到了重建的曲面，参见《High Quality Pose Estimation by Aligning Multiple Scans to a Latent Map》；
- 扫描数据与潜在曲面做配准，隐含着在不同扫描数据之间做配准；
- 通过空间网格的划分自动决定了哪些数据之间需要做匹配，于是不需要指明扫描数据的重叠区域；

同时配准和重建，优化过程中的两个步骤

5 其它的配准方法

非刚体配准（non-rigid registration）：
上文介绍的算法都假设输入数据通过刚体变换相联系，而非刚性配准允许输入数据之间存在形变；
基于概率模型的方法：
- 把扫描数据视为概率分布，通过对准概率分布来实现配准；
- 在医疗图像中用得比较多；
其它的基于学习的方法；

参见《Registration of 3D point clouds and meshes: A survey from rigid to nonrigid》。
阅读材料

《A Method for Registration of 3-D Shapes》Paul J. Besl and Neil D. McKay. IEEE Transactions on Pattern Recognition and Machine Intelligence. 1992.
《Efficient Variants of the ICP Algorithm》Szymon Rusinkiewicz and Marc Levoy. 2001.
《Geometry and Convergence Analysis of Algorithms for Registration of 3D Shapes》Helmut Pottmann, Qixing Huang, Yongliang Yang, and Shimin Hu. International Journal of Computer Vision. 2006.
《High Quality Pose Estimation by Aligning Multiple Scans to a Latent Map》Qixing Huang and Dragomir Anguelov. IEEE International Conference on Robotics and Automation 2010.
《DynamicFusion: Reconstruction and Tracking of Non-rigid Scenes in Real-Time》Richard Newcombe, Dieter Fox, and Steve Seitz. CVPR 2015.

其它材料：

《数字几何处理研究进展》胡事民, 杨永亮, 来煜坤. 计算机学报, 2009.
《计算机视觉基础》，由鲁鹏教授开设，录播视频。

本节录播视频

		自动登录	找回密码
密码			立即注册

《GAMES203：三维重建和理解》2 配准（registration）

本帖子中包含更多资源

浏览过的版块