论文翻译—ICP算法论文（节选1~4章）

标题：A Method for Registration of 3-D Shapes
作者： J. Besl, Neil D. McKay
目录：

• Ⅰ Introduction
• Ⅱ Literature Review
• Ⅲ Mathematical Preliminaries
  • A. Point to Parametric Entity Distance
  • B. Point to Implicit Entity Distance
  • C. Corresponding Point Set Registration
• Ⅳ The Iterative Closest Point Algorithm
  • A. ICP Algorithm Statement
  • B. Convergence Theorem
  • C. An Accelerated ICP Algorithm
  • D. Alternative Minimization Approaches
• Ⅴ The Set of Initial Registrations
  • A. Initial States for Global Matching
  • B. A Counter Example
  • C. Local Shape Matching
• Ⅵ. Experimental Results
  • A. Point Set Matching
  • B. Curve Matching
  • C. Surface Matching
    • 1 A Bezier Surface Patch
    • 2 The NRCC African Mask
    • 3 Terrain Data
• Ⅶ. Conclusions
• Ⅷ. Future Directions

一种用于3D形状配准的方法

摘要——本文描述了一种通用的、独立于表达的方法，用于准确且高效计算的 3D 形状（包括自由曲线和曲面）配准。该方法基于迭代最近点（iterative closest point, ICP）算法处理整个6自由度，ICP算法只需要一个找到几何实体上离给定点最近的点的过程。ICP算法总是单调地收敛到均方距离度量的最小值，且经验显示在前面的迭代中收敛速率很快。因此对于具有一定“形状复杂性”级别的特定类别的物体，给定足够的初始旋转和平移的集合，是可以通过测试每个初始配准来全局地在整个6自由度上最小化均方距离度量。例如，一个给定形状的“模型”和一个表达模型形状主要部分的感知的“数据”形状，能够以允许给定的模型复杂度，通过测试一个初始平移和一组相对较小的旋转，在几分钟内配准起来。该方法的一个重要应用是将来自不固定刚性物体的感知数据与理想几何模型的先验进行配准，来做形状检查。我们描述的方法也能够被用来决定基本问题，如，不同几何表达的一致（形状等价）以及估计未知对应的点集间的运动。实验结果展示了配准算法在点集，曲线和曲面上的能力。

关键词——自由曲线匹配，自由曲面匹配，运动估计，位姿估计，四元数，3D配准

1. 引言

对自由曲线，表面和点集的全局和局部形状匹配度量在【3】中有提到，它试图形式化和统一计算机视觉中关键问题的描述：给定一个在传感器坐标系下的3D数据，这描述了一个数据形状与模型形状的对应，并且给定一个以不同几何形状表达的模型坐标系中的模型形状，估计最优的旋转和平移，从而对齐，或者配准的模型形状与数据形状，这是通过最小化形状间的距离实现的，且因此能够通过均方误差度量来确定形状间的等价。许多应用的感兴趣的关键问题是：将来自距离图像的分段区域与在计算机辅助设计（computer-aided-design(CAD)）模型上的B样条表面进行匹配是否可行？本文提出了一个自由表面匹配问题的解决方案，正如在【3】和【5】中定义的作为一个简单、通用且统一方法的特殊情况，它可以泛化到 $n$ 维并为 1）未知对应的点集匹配问题；2）自由曲线匹配问题；提供解决方案。算法不需要额外的特征，不需要曲线或表面的导数，且除了去除统计的异常值，不需要对3D数据进行预处理。

    我们提出的方法的主要应用是将来自未固定的刚性物体的数字化数据与理想的几何模型先验进行匹配，来做形状检查。在浅景深范围内使用高精度非接触式测量设备来检测形状时，不同感测点的不确定度变化不大。因此，为简单起见并与基于数千个数字化点的检测应用相关，不考虑点之间不同的不确定性的情况。类似地，统计离群值的去除被认为是一个预处理步骤，可能是最好的实现，也不会被处理。在检查应用的背景下，高精度非接触性测量设备不生成坏的数据的假设是合理的，因为一些传感器有能力拒绝高度不确定的测量值。

    我们提出的形状配准算法可以在以下的几何数据的表达中被使用：

1）点集
2）线段集（折线）；
3）隐式曲线：$\vec{g}(x, y, z) = 0$
4）参数化曲线：(x(u), y(u), z(u))
5）三角集合（刻面）
6）隐式表面：$g(x, y, z) = 0$
7）参数化表面：(x(u, v), y(u, v), z(u, v))

这覆盖了要使用3D形状配准方法的大部分应用。其他表达是通过提供评估给定形状上与给定数字化点最近的点的过程来处理的。

    本文结构如下：首先回顾一些文献中的相关论文。其次会提到的几何表示涵盖了计算形状上到给定点的最近点的数学预备知识。然后，描述了迭代最近点（ICP）算法，以及关于它单调收敛性的定理的证明。最后，基于点集，曲线，表面的实验结果展现了ICP配准算法的能力。

2. 文章预览

    在 3D 自由形状的配准（姿势估计、对齐、运动估计）领域,发表的工作相对较少。大部分解决全局形状匹配或配准的现有文献，已经解决了有限的形状类别，即，1）多面体模型，2)分段（超​​）二次模型【2】【27】，3）未知对应（关系）的点集。读者可以参考【6】和【14】了解1985年以前这些领域的工作。对于以下未讨论的其他近期相关工作的样本，见【8，10，12，13，19，20，24，26，34，35，37，39，44，46，48，53，58，59】。

    从历史上看，使用 3-D 数据的自由形状匹配最早是由 Faugeras 和他在 INRIA 的团队完成的【18】，他们在 1980 年代初期展示了与雷诺汽车零件（转向节）的有效匹配。这项工作在计算机视觉社区中推广了使用四元数进行相应 3-D 点集的最小二乘配准。奇异值分解（SVD）算法的替代使用【23，1，49】在这个时间范围内并不广为人知。这项工作的主要限制是它依赖于在自由形状中可能存在的相当大平面区域。

    Schwartz 和 Sharir【50】在1985年年底开发了一个不需要特征提取的自由形状空间曲线匹配问题的解决方案。他们使用非四元数方法来计算最小二乘旋转矩阵。该方法在合理的质量曲线数据上工作很好，但在有很大噪声的曲线（数据）上有困难，因为该方法使用曲线的弧长采样来获得对应的点集。

    Haralick等【28】使用稳健方法结合最小二乘 SVD 配准方法解决了 3-D 点集姿态估计问题，这为最小二乘四元数或 SVD 点集匹配提供了稳健的统计替代方案。该算法能够处理统计异常值并且在理论上，只要正交矩阵的行列式严格为正数，就可以替换我们基于四元数的算法。最近的会议论文集【47】涵盖了关于该主题的新的贡献。

    Horn【31】衍生出了最小二乘四元数匹配的 Faugeras’s 方法【18】的可替代的公式，它使用 $4\times 4$ 矩阵最大特征值代替最小特征值。Horn【30】和 Brou【11】也开发出了扩展高斯图像（extended Gaussian image, EGI）的方法，该方法允许基于表面法线直方图的凸和受限非凸形形状集的匹配。

    Taubin【55】在隐式代数非平面3-D曲线和曲面估计（应用于无需特征提取的位置估计）领域做了一些有趣的工作。他描述了一个使用近似距离度量高达10度的使用隐式代数近似数据点的方法（he describes a method of approximating data points with implicit algebraic forms up to the tenth degree using an approximate distance metric）。全局形状（没被遮挡的形状）可以根据广义特征值进行识别，并且可以恢复配准变换。该方法被证明可用于完整的平面曲线和空间曲线形状，但尚不清楚其有效性能否很好地推广到更复杂的表面，例如地形数据或人脸。Taubin 曾表示近似距离拟合的数值方法倾向于在十度以上分解（tend to break down above the tenth degree）。他后来【56】通过研究基于广义形状多项式的形状匹配，扩展了他在形状描述方面的工作。这表明了一些有趣的理论结果，但是在复杂表面上的实际使用仍有待证明。

    Szeliski【54】描述了一种从稀疏范围数据中估计运动的方法，没有点之间的对应关系，也没有特征提取。他的主要目标是创建一种方法来估计观察者在同一地形的两个距离图像帧之间的运动。给定一帧的点集，他应用平滑假设来创建点的平滑样条近似。然后，应用传统的最速下降算法来旋转并平移第二个数据集，使得它最小化点和表面之间的协方差加权的 $z$ 差异的总和。他的方法是基于规则的 $xy$ 网格结构，且没有计算真实的3-D点到表面的距离。最速下降方法是一个比我们下面提出的ICP算法到达局部极小值更慢的选项。Szeliski 使用最优贝叶斯数学，让他能够从模拟测距仪中降低更远距离的噪声值，在数据点中的不确定性从前景到背景变化很大。对于浅景深的高精度传感器，点之间的不确定性变化要小几个数量级，并且不太受关注。Szeliski 提供了合成地形数据和块（block）的实验结果。地形数据运动测试是一个沿着一个轴的简单的平移：一个1-D的相关问题。他的块测试涉及到6自由度，但是块是一个非常简单的形状。总之，这项工作提出了一些有趣的创新点，但是实验结果在应用中没有说服力。

    Horn 和 Harris【33】也解决了在给定相同地形的有序数字化距离图像帧的情况下，估计观察者的精确的刚体运动的问题。他们描述了一个射程率约束方程和海拔率约束方程。结果是一种非迭代最小二乘法，只要数据帧之间的运动相对较小，就可以提供六自由度运动估计。这个方法比 Szeliski 提出的方法要快得多，但是尚不清楚这种方法是否可以推广到形状的任意旋转和平移。

    Kamgar-Parsi等【36】也描述了一种无需提取特征的多重叠距离图像配准方法。该方法使用*工作地非常好。使用 2.5-D 距离数据的水平集（level sets），但基本上仅限于平面中的三个自由度，因为这项工作旨在拼凑地形地图数据。

    Li【38】解决了任意旋转和平移的自由形式表面匹配问题。他的方法为数据和模型形状形成了基本表面区域的属性关系图，然后使用不精确的方法执行图匹配，该方法允许属性和图邻接关系的可变性。这似乎是一个合理的方法，但它依赖于基于导数的数量的提取（relies on extraction of derivative-based quantities）。实验结果显示为咖啡杯和雷诺汽车零件；见Wong等【60】其它的相关工作使用用于3-D匹配的属性图。

    Gilbert 和 Foo【21】和 Gilbert等【22】的工作是相关的，它解决了两个目标形状之间距离的计算。此类方法可以作为下面描述的类似形状匹配技术的基础。然而，他们的方法的主要不便之处在于，必须将对象形状分解为凸的亚体（convex subbodies），这个问题通常对于任意 CAD 模型或数字化 3-D 数据来说并非微不足道。

3. 数学预备

    在这一章中，我们描述了在上面列出的各种几何表达上计算到给定点最近的点的方法。首先，覆盖基本几何实体，然后是参数实体，最后是隐式实体。读者可以参考 Mortenson【42】下面的一些项目以获取更多信息。

    两点 ${\vec{r}}_1=(x_1,y_1,z_1)$ 和 ${\vec{r}}_2=(x_2,y_2,z_2)$ 之间的欧式距离 $d({\vec{r}}_1,{\vec{r}}_2)=\Vert {\vec{r}}_1,{\vec{r}}_2\Vert =\sqrt{(x_2-x_1{)}^2+(y_2-y_1{)}^2+(z_2-z_1{)}^2}$。令 $A={\vec{a}}_i,(i=1,...,N_a)$ 表示具有 $N_a$ 个点的点集。定义点 $\vec{p}$ 和点集 $A$ 之间的距离

d ( p ⃗ , A ) = m i n i ∈ { 1 , . . . , N a } d ( p ⃗ , a ⃗ i )   (1) d(\vec{p}, A) = \underset{i \in \{1,...,N_a\}}{min} d(\vec{p}, \vec{a}_i)~\tag{1} d(p ​,A)=i∈{ 1,...,Na​}min​d(p ​,a i​) (1)

$A$ 的最近点 ${\vec{a}}_j$ 满足等式 $d(\vec{p},{\vec{a}}_j)=d(\vec{p},A)$。

    令 $l$ 为连接两个点 ${\vec{r}}_1$ 和 ${\vec{r}}_2$的线段。定义点 $\vec{r}$ 与线段 $l$ 之间的距离

$$d(\vec{p},l)=\underset{u+v=1}{\min }\Vert u{\vec{r}}_1+v{\vec{r}}_2-\vec{p}\Vert \text{\hspace{1em}(2)}$$

其中 $u\in [0,1],v\in [0,1]$。需要的闭式计算很直接。令 L = { l i } ,   ( i = 1 , . . . , N l ) L = \{l_i\},~(i = 1,...,N_l) L={ li​}, (i=1,...,Nl​) 为 $N_l$ 条线段 $l_i$ 的集合。定义点 $\vec{p}$ 和线段集 $L$ 的距离

d ( p ⃗ , L ) = m i n i ∈ { 1 , . . . , N l } d ( p ⃗ , l i )   (3) d(\vec{p}, L) = \underset{i \in \{1, ... , N_l\}}{min} d(\vec{p}, l_i)~\tag{3} d(p ​,L)=i∈{ 1,...,Nl​}min​d(p ​,li​) (3)

    令 $t$ 为由三个点 ${\vec{r}}_1=(x_1,y_1,z_1),r_2=(x_2,y_2,z_2),r_3=(x_3,y_3,z_3)$ 定义的三角形，定义点 p 到三角形 $t$ 间的距离

$$d(\vec{p},t)=\underset{u+v+w=1}{\min }\Vert u{\vec{r}}_1+v{\vec{r}}_2+w{\vec{r}}_3-\vec{p}\Vert \text{\hspace{1em}(4)}$$

其中，$u\in [0,1],v\in [0,1],w\in [0,1]$ 需要的闭式计算也很直接。令 T = { t i } ,   ( i = 1 , . . . , N t ) T = \{t_i\},~(i = 1, ... , N_t) T={ ti​}, (i=1,...,Nt​) 为 $N_t$ 个三角形 $t_i$ 的集合。定义点 $\vec{p}$ 和三角形集 $T$ 的距离

d ( p ⃗ , T ) = m i n i ∈ { 1 , . . . , N t } d ( p ⃗ , t i )   (5) d(\vec{p}, T) = \underset{i \in \{1, ... , N_t\}}{min} d(\vec{p}, t_i)~\tag{5} d(p ​,T)=i∈{ 1,...,Nt​}min​d(p ​,ti​) (5)

$y_i$ 到三角形集上的最近点满足等式 $d(\vec{p},{\vec{y}}_i)=d(\vec{p},T)$

3.1 点到参数化实体距离

在这一节中，参数曲线和参数表面到被视为单个的参数实体 $\vec{r}(\vec{u})$ 其中 $\vec{u}=u\in R^1$ 应替换为参数曲线，而 $\vec{u}=(u,v)\in R^2$ 应替换为参数曲面（R表示实线）。区间内曲线的求值域，但曲面的求值域可以是平面中的任意闭合连接区域。关于参数实体的更多信息，如 贝塞尔和B样条曲线和曲面，见【9，15-17，42，52】

    从一个给定点到一个参数实体的距离 $E$ 为

$$d(\vec{p},E)=\underset{\vec{r}(\vec{u})\in E}{\min }d(\vec{p},\vec{r}(\vec{u}))\text{\hspace{1em}(6)}$$

计算距离的计算不是闭式的，而是相对复杂的。下面描述了一个计算点到曲线和点到曲面距离的方法。一旦实现了单个实体的距离度量，参数化实体集再次变得简单明了。令 F = { E i } ,   ( i = 1 , . . . , N e ) F = \{E_i\},~(i = 1, ... ,N_e) F={ Ei​}, (i=1,...,Ne​) 表示 $N_e$ 个由 $E_i$ 表示的参数实体。定位点 $\vec{p}$ 到参数实体集 $F$ 的距离

d ( p ⃗ , F ) = m i n i ∈ { 1 , . . . , N e } d ( p ⃗ , E ⃗ i )   (7) d(\vec{p}, F) = \underset{i \in \{1, ... ,N_e\}}{min} d(\vec{p}, \vec{E}_i)~\tag{7} d(p ​,F)=i∈{ 1,...,Ne​}min​d(p ​,E i​) (7)

在参数实体集 $F$ 上的最近点 ${\vec{y}}_i$ 满足等式 $d(\vec{p},{\vec{y}}_j)=d(\vec{p},F)$

    计算从一个点到一个参数实体的距离的第一步是创建一个基于单纯形的近似（线段或三角形）。对于一个参数空间曲线 C = { r ⃗ ( u ) } C = \{\vec{r}(u)\} C={ r