最近3维(${\rm 3D}$)模型越来越引起大家的关注。它可以广泛地应用在娱乐业的${\rm 3D}$游戏或者${\rm 3D}$电影^[1]、${\rm 3D}$测量^[2]、${\rm 3D}$上颌牙齿曲线扫描^[3]、机器人视觉引导^[4]等。点云数据的获取通常有以下几种方式：接触式^[5]、非接触式^[6-7]或者断层扫描^[8]等。然而，无论采用哪种方式，均需要利用从多个角度采集的点云数据才能构建物体完整的3维模型。不同的采集角度意味着点云数据存在于不同的坐标系，因而需要对不同坐标系下的点云数据进行拼合处理，其实质就是将来自不同坐标系的点云数据放在同一个坐标系下。

对于大多数商用级的${\rm 3D}$激光扫描仪系统，拼合的过程大都采用非自动化的方式。比如可以通过将标记贴在被扫描的物体上，或通过预先校准两次扫描之间的浮动来实现。另外，也可以采用自动拼合的方式。然而，自动拼合算法由于稳定性差、准确度低、低速等缺点，其应用受到一定限制。典型的点云自动拼合方式可以分为两类：直接基于点云数据本身的^[9]和基于几何特征的^[10-11]的拼合方式。

直接的拼合方式，就是将不同视角的点云直接配准在一起。一种典型的配准方式就是最近点迭代(${\rm ICP}$)算法^[12]。${\rm ICP}$算法不需指定任何特征，但是如果两片点云数据没有较好的初始位置，根据距离最近原则确定的匹配点会包含很多错误的对应关系，导致最终往往只能得到局部最优解，而不是全局最优解。而且${\rm ICP}$算法要求待拼接的一片点云完全包含另外一片，对部分重叠或者不重叠的问题无法解决。${\rm Chen}$等人^[13]提出使用两个法向量的距离来代替两个点的距离作为最近点搜索的评价函数。但此方式计算效率低，因为它需要求解非线性的最小二乘问题，并且它要求两片点云足够靠近。为了解决非线性最小二乘问题，${\rm Masuda}$等人^[14]引入了最小中值平方估计(${\rm LMedS}$)策略，但是它的每次迭代都需要重采样。

另一种拼合方式是基于特征点的，特征点可以是洞、边沿、拐角、法向量等。${\rm Jiang}$等人^[15]提出根据$k$维角不变特征选择对应点，这种特征由点法向量及其邻域点法向量之间的夹角形成的$k$维角度组成，它对于缩放和旋转变换是不变的，并且适用于具有较小曲率的表面。辛伟等人^[16]通过采用不规则曲面点与其一邻域点重心间的距离作为刚性特征对点云进行拼合，此方法对于点云的初始位置有一定的要求。${\rm Chen}$等人^[17]提出了一种基于主曲率分布的局部表面描述因子来实现点云匹配的方法。${\rm Zhu}$等人^[18]提出了基于曲率和法向量的${\rm hash}$拼合算法。为了使得拼合更加的准确，${\rm Akca}$等人^[19]提出了归一化其他的一些信息比如：密度，颜色，质地等。然而，这些信息对采集设备有很高的要求。基于特征点的方式没有任何的初始位置限制可以很容易地实现拼合并且达到较高的拼合精度，并且适用不同密度的点云数据，然而却无法适用在没有明显特征的物体上。

最近，在这两类方法的研究上取得了重要进展。${\rm Chen}$等人^[20]提出了${\rm HT-ICP}$算法。该算法改变传统${\rm ICP}$中的原始非线性优化问题，通过参数化具有扭曲指数的刚性运动将欧几里得变换导出为矩阵方程(其可以利用高斯消去策略方便地求解)。为了解决传统${\rm ICP}$算法部分重叠的问题，提出了一种基于估计的重叠率和可接受的间隔来去除在非重叠区的错误点对的策略。与传统的${\rm ICP}$及其变种相比，该算法同时提高了拼合的速度和精度。${\rm He}$等人^[21]提出了一种基于特征点的算法，该算法适用于部分重叠的点云数据的拼合。在此算法中，对应点关系的建立基于两个步骤：1)通过计算具有最大局部曲率变化的特征点之间的${\rm Hausdorff}$曲率距离来获得初始匹配点；2)基于局部表面的圆周形状特征从初始匹配点获得精确匹配点。基于准确的匹配点，该算法使用四元素迭代方式获得欧几里得变换。${\rm Papazov}$等人^[22]提出了一种拼合算法，基于反向距离将拼合问题转化为两个点云非线性代价函数的最小化问题，这可以显著减小噪声和边界的影响。在不需要任何初始对准的情况下，使用随机方法使成本函数全局最小化，该方法使用广义${\rm BSP}$树并且允许在复杂形状的搜索空间(例如旋转空间)上最小化非线性标量场。

在基于特征的拼合方法中，通常需要用到微分信息。但是这类方法对噪声非常敏感，容易产生错误的匹配点对，同时对重叠率有较高要求。为此，提出一种逐层筛选的匹配点对计算策略，并利用霍夫变换确定最优的坐标变换参数，最终实现点云拼合。首先利用点云的曲率相似度确定初始匹配点对的集合，然后利用刚体不变量特征相似度对初始匹配点对进行筛选。由于曲率和刚体不变特征都是基于局部求导的，对噪声很敏感，所以匹配点对数组中会存在很多错误点对。为了减小噪声对匹配点对的影响，利用持续特征直方图(${\rm PFH}$)对错误的匹配点对进行进一步的筛选。为了进一步提高匹配的正确率，消除错误匹配点对对坐标变换参数估计的影响，利用霍夫变换实现聚类分析，从而估算出最优的坐标变换参数并实现点云拼接。实验结果表明，本文算法不仅能适用于具有任意初始相对位置的部分重叠点云的拼接，而且可以取得很高的拼合精度和噪声鲁棒性。

基于刚体不变量特征值和持续特征直方图进行点云拼合，经过霍夫变换得到变换参数。在这个过程中，曲率特征和刚体不变量特征是基于球面映射的^[23]，通过此方法来提取初始匹配点对，然后使用直方图特征值^[24]对匹配点对进行筛选，最后通过霍夫变换进一步筛选出正确的匹配点对，从而实现点云拼合。

对于规则排列的数据点可以进行二次参数曲面^[25]拟合来估算法向量和曲率。由于很难对散乱无序的点云数据进行参数化，因此采用局部切平面参数化^[26]来解决该问题，即用投影到局部切平面的点代替原始数据点来进行参数化。然而，该方法计算出来的法向量方向往往是不一致的，所以需要对法向量方向进行调整。为了解决这一问题，${\rm Hoppe}$等人^[27]提出一种法向量方向“传播”策略：选取点云中的一点作为调整的初始点，按照一定的优先原则对该点的某个近邻点进行调整；然后将已完成法向量调整的点的法向量方向传播给邻近点的法向量，直到调整完所有数据点的法向量方向为止。本文在选择近邻点进行法向量方向调整时采用广度优先结合$cost$权值的优先原则。其中，$cost = 1-\left| {\mathit{\boldsymbol{N}}\left( {{p_i}} \right) \cdot \mathit{\boldsymbol{N}}\left( {{q_j}} \right)} \right|$，权值越小优先级越高。在此过程中，一个点$q_{j}∈Knn(p_{i})(Knn(p_{i})$表示点$p_{i}$的$k$-近邻像素组成的集合)是否需要进行法向量的调整依赖于两个法向量$\mathit{\boldsymbol{N}}(p_{i})$和$\mathit{\boldsymbol{N}}(q_{j})$之间的内积，如果$\mathit{\boldsymbol{N}}(p_{i})·\mathit{\boldsymbol{N}}(q_{j})<0$，则需要进行法向量的调整。

对点云$\mathit{\boldsymbol{P}}$中任意点$p_{i}$，其$r$-闭球域的球面映射为

$ {\mathit{\boldsymbol{G}}_r}\left( {{p_i}} \right) = \left\{ {\mathit{\boldsymbol{N}}\left( q \right)\left| {q \in \mathit{\boldsymbol{R}}\left( {{p_i}} \right)} \right.} \right\} $

(1)

式中，$\mathit{\boldsymbol{R}}\left( {{p_i}} \right) = \left\{ {q|{{\left\| {q-{p_i}} \right\|}_2} < r, q \in \mathit{\boldsymbol{P}}} \right\}$是点$ p_i$的$r$-闭球域，$\mathit{\boldsymbol{N}}(q)$是$q$点处的法向量。以$\mathit{\boldsymbol{G}}_{r}(p_{i})$中所有点的偏差向量为列向量构成如下$3×K$的偏差矩阵

$ \mathit{\boldsymbol{D}}\left( {{p_i}} \right) = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{N}}\left( {{q_1}} \right) - g\left( {{p_i}} \right)}\\ {\mathit{\boldsymbol{N}}\left( {{q_2}} \right) - g\left( {{p_i}} \right)}\\ \vdots \\ {\mathit{\boldsymbol{N}}\left( {{q_K}} \right) - g\left( {{p_i}} \right)} \end{array}} \right] $

(2)

式中，$g\left( {{p_i}} \right) = \frac{1}{K}\sum\limits_{k = 1}^K {\mathit{\boldsymbol{N}}\left( {{q_k}} \right)} $是$\mathit{\boldsymbol{G}}_{r}(p_{i})$的质心，$K$为$r$-闭球域中点的个数。从偏差矩阵的转置矩阵$\mathit{\boldsymbol{D}}{\left( {{p_i}} \right)^{\rm{T}}}$可以获得其奇异值$σ_{1}≤σ_{2}≤σ_{3}$，及其对应的右奇异矢量$\mathit{\boldsymbol{v}}_{1}、\mathit{\boldsymbol{v}}_{2}、\mathit{\boldsymbol{v}}_{3}$。考虑到这3个奇异值对刚体变换来说是不变量，将不变特征标识定义为

$ f\left( {{p_i}} \right) = \left\{ \begin{array}{l} {\left[ {\frac{{{\sigma _1}}}{{{\sigma _3}}},\frac{{{\sigma _2}}}{{{\sigma _3}}}} \right]^{\rm{T}}}\;\;\;\;\;{\sigma _3} \ne 0\\ {\left[ {1,1} \right]^{\rm{T}}}\;\;\;\;\;\;\;\;\;\;\;\;{\sigma _3} = 0 \end{array} \right. $

(3)

这里，利用奇异值的比值来定义不变特征标识，不仅仅使得它与坐标系无关，而且也可以消除点数量的影响。它可以很好地体现出曲面的特性，而且减少噪声的影响，增强鲁棒性。

假设目标点云$\mathit{\boldsymbol{P}}=\{p_{i}, i=1, 2, …, n_{1}\}$和浮动点云$\mathit{\boldsymbol{Q}}=\{q_{j}, i=1, 2, …, n_{2}\}$是两个来自不同视角的点云数据。实现$\mathit{\boldsymbol{P}} $和$ \mathit{\boldsymbol{Q}}$对齐配准的关键是在两片部分重叠的点云中找到正确的匹配点对。假设$p_{i}∈\mathit{\boldsymbol{P}}$和$q_{j}∈\mathit{\boldsymbol{Q}}$是一对匹配点对，那么这两个点的邻域曲面也应该相似，即满足

$ \begin{array}{l} {\mathop{\rm sgn}} \left( {{C_g}\left( {{p_i}} \right)} \right) = {\mathop{\rm sgn}} \left( {{C_g}\left( {{q_j}} \right)} \right)\\ {\mathop{\rm sgn}} \left( {{C_a}\left( {{p_i}} \right)} \right) = {\mathop{\rm sgn}} \left( {{C_a}\left( {{q_j}} \right)} \right) \end{array} $

(4)

式中, $C_{a}$和$C_{g}$分别是点$p_{i}$或者$q_{j}$处的平均曲率和高斯曲率。由于曲面平坦区域的特征不够明显，容易造成误匹配，因此在寻找对应点对时设置一个充分小的阈值$ε_{1}>0$用以排除掉平坦区域的点，即

$ \begin{array}{l} \left| {{C_a}\left( {{p_i}} \right) \cdot {C_g}\left( {{p_i}} \right)} \right| > {\varepsilon _1}\\ \left| {{C_a}\left( {{q_j}} \right) \cdot {C_g}\left( {{q_j}} \right)} \right| > {\varepsilon _1} \end{array} $

(5)

为了进一步剔除错误的匹配点对，本文引入点曲率相似度过滤对应点对，即

$ \begin{array}{l} \left| {\frac{{{K_1}\left( {{p_i}} \right) - {K_1}\left( {{q_j}} \right)}}{{{K_1}\left( {{p_i}} \right) + {K_1}\left( {{q_j}} \right)}}} \right| \le {\varepsilon _2}\\ \left| {\frac{{{K_2}\left( {{p_i}} \right) - {K_2}\left( {{q_j}} \right)}}{{{K_2}\left( {{p_i}} \right) + {K_2}\left( {{q_j}} \right)}}} \right| \le {\varepsilon _2} \end{array} $

(6)

式中，$ε_{2}>0$是另一个预设值，$K_{1}$和$K_{2}$分别是点$p_{i}$或者$q_{j}$处的主曲率。以上公式通常被称为曲率相似性^[28]。

从理论上来说，如果两个点是一对正确的匹配点对，那么它们不仅仅在曲率上有相似性，而且其各自邻域内对应点的不变量特征标识也应该相似。利用点云特征不变量标识间的归一化互相关系数(${\rm NCC}$)来度量邻域标识相似度，从而对匹配点对进一步过滤。设$Ncc(p_{i}, q_{j})$为衡量点对($p_{i}, q_{j}$)的邻域标识相似度的归一化互相关系数，具体公式为

$ \begin{array}{*{20}{c}} {Ncc\left( {{p_i},{q_j}} \right) = }\\ {\frac{{\sum\limits_{l = 1}^k {\left( {f\left( {{p_l}} \right) - {c_{f{p_i}}}} \right) \cdot \left( {f\left( {{q_l}} \right) - {c_{f{q_j}}}} \right)} }}{{\sqrt {\sum\limits_{l = 1}^k {{{\left\| {\left( {f\left( {{p_l}} \right) - {c_{f{p_i}}}} \right)} \right\|}^2}} } \sqrt {\sum\limits_{l = 1}^k {{{\left\| {\left( {f\left( {{q_l}} \right) - {c_{f{q_j}}}} \right)} \right\|}^2}} } }}} \end{array} $

(7)

式中，$f(p_{l})$和$f(q_{l})$是$点p_{l}∈Knn(p_{i})$和$q_{l}∈Knn(q_{j})$的刚体不变特征，$c_{fp_{i}}$和$c_{fq_{j}}$分别是$Knn(p_{i})$和$Knn(q_{j})$刚体不变特征的质心。显然，$Ncc(p_{i}, q_{j})$值越大则点对$(p_{i}, q_{j})$的邻域标识相似度越高。本文设定一个阈值$ε_{3}>0$，当$Ncc(p_{i}, q_{j})>ε_{3}$时，认为点对$(p_{i}, q_{j})$为一组正确匹配点对。首先通过点的曲率相似度初步排除高斯曲率和平均曲率相差较大的点，然后引入邻域标识相似度进一步剔除错误匹配点对，提高匹配正确率。

使用持续特征直方图(${\rm PFH}$)对匹配点对进行进一步的筛选。设点对$ (p_{i}, p_{j})(i≠j, i < j)$为点云数据中位于点$p$的$k$-邻域中的点对，其对应的法向量为$({\mathit{\boldsymbol{n}}_\mathit{\boldsymbol{i}}}, {\mathit{\boldsymbol{n}}_\mathit{\boldsymbol{j}}})$。选择其中一个为源点$p_{s}$，另一个为目标点$p_{t}$。选择源点$p_{s}$的具体公式为

$ \left( {{p_s},{p_t}} \right) = \left\{ \begin{array}{l} \left( {{p_i},{p_j}} \right)\;\;\;\;\alpha \le \beta \\ \left( {{p_j},{p_i}} \right)\;\;\;\;\alpha > \beta \end{array} \right. $

(8)

式中，$α=〈n_{i}, p_{j}-p_{i}〉，β=〈n_{j}, p_{i}-p_{j}〉$。

为了计算持续特征量，定义如图 1所示的达布架构^[24]。

图 1中，$p_{s}$为坐标原点，3个坐标轴定义为

$ \mathit{\boldsymbol{u}} = {\mathit{\boldsymbol{n}}_s},\mathit{\boldsymbol{v}} = \left( {{\mathit{\boldsymbol{p}}_t} - {\mathit{\boldsymbol{p}}_s}} \right) \times \mathit{\boldsymbol{u}},\mathit{\boldsymbol{w}} = \mathit{\boldsymbol{u}} \times \mathit{\boldsymbol{v}} $

(9)

利用上述达布架构，计算点对$(p_{i}, p_{j})$的持续特征量为

$ {p_f} = \sum\limits_{i = 1}^4 {step\left( {{s_i},{f_i}} \right) \cdot {2^{i - 1}}} $

(10)

式中

$ step\left( {{s_i},{f_i}} \right) = \left\{ \begin{array}{l} 0\;\;\;\;{f_i} < {s_i}\\ 1\;\;\;\;{f_i} > {s_i} \end{array} \right. $

(11)

$ \begin{array}{*{20}{c}} {{f_1} = \left\langle {\mathit{\boldsymbol{v}},{\mathit{\boldsymbol{n}}_t}} \right\rangle }\\ {{f_2} = \left\| {{p_t} - {p_s}} \right\|}\\ {{f_3} = \left\langle {\mathit{\boldsymbol{u}},{p_t} - {p_s}} \right\rangle /{f_2}}\\ {{f_4} = a\tan \left( {\left\langle {\mathit{\boldsymbol{w}},{\mathit{\boldsymbol{n}}_t}} \right\rangle ,\left\langle {\mathit{\boldsymbol{u}},{\mathit{\boldsymbol{n}}_t}} \right\rangle } \right)} \end{array} $

(12)

因为$f_{1}，f_{3}$是归一化向量间的点积，所以它们实际上就是3维向量间夹角的余弦值，因而其范围是$±1$。类似地，$f_{4}$是${\mathit{\boldsymbol{n}}_t}$在($\mathit{\boldsymbol{u}}, \mathit{\boldsymbol{w}}$)平面上的投影与$\mathit{\boldsymbol{u}}$轴之间的夹角，所以它的取值范围是$±{\rm \mathsf π}$。

显然，通过式(10)可以将点对$(p_{i}, p_{j})$的4个特征二值化，而式(11)的作用则是将经过二值化处理的4个的特征编码成0到15之间的整数。通过这样的处理，可以很方便地构造一个16-${\rm bin}$的直方图${\rm PFH}$，用以对点云数据点的邻域中所有点对的特征出现的频数进行统计。

利用${\rm PFH}$进一步剔除一些错误的匹配点。设$(p, q)$为经过1.2.2节和1.2.3节中筛选过程得到的一对匹配点，首先分别计算这两点处的$PFH(p)$和$PFH(q)$，然后计算这两个直方图间的欧氏距离$||PFH(p)-PFH(q)||_{2}$，设定一个阈值$ε$，若$||PFH(p)-PFH(q)||_{2}<ε$，则保留该匹配点，否则剔除该匹配点。

经过基于曲率相似度、邻域标识相似度、持续直方图的筛选之后，可以得到一个初步的候选匹配点对集合

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{Match}} = \left\{ {\left( {{{p'}_i},{{p'}_i}} \right)\left| {{{p'}_i} \in \mathit{\boldsymbol{P}},{{p'}_i} \in \mathit{\boldsymbol{Q}}} \right.,} \right.}\\ {\left. {i = 1,2, \cdots ,{N_m}} \right\}} \end{array} $

(13)

式中，$N_{m}$是点对数组$\mathit{\boldsymbol{Match}}$中点对的数量，$(p′_{i}, p′_{i})$是浮动点云$\mathit{\boldsymbol{Q}}$和目标点云$\mathit{\boldsymbol{P}}$中的对应点对。利用候选匹配点对集合可以找到坐标变换的参数从而实现点云拼接。

一个匹配点对对应的是一组坐标系的变换参数，$N_{m}$个点对对应的是$N_{m}$个坐标系的变换参数。虽然经过二次筛选，在匹配点对中仍然存在很多错误的匹配点对。这些错误的匹配点对来自于$\mathit{\boldsymbol{P}}$和$\mathit{\boldsymbol{Q}}$(如图 2)的重叠区域或者非重叠区域。由于存在错误的匹配点对，由不同的点对计算得到的变换参数会离散地分布在一定范围内而不是集中到一点。因此，需要从目标点云和浮动点云中的$N_{m}$离散分布的变换参数中提取正确的变换参数。由于在匹配点集中绝大部分点对是正确匹配的，因此由这些点对得到的变换参数应该是一致或近似一致的。基于此，并受基于霍夫变换的线检测思想的启发，本文将这些坐标变换参数在一个合适的空间表达出来，从而使得由正确的匹配点对得到的变换参数在此空间中会聚集在一起。这样，目标点云和浮动点云间的坐标变换参数即为具有最大计数元胞所对应的参数。这个策略可以有效地解决错误的点对对最终的变换参数的影响。

通过偏差矩阵$\mathit{\boldsymbol{D}}(p_{i})$来定义任意点$p_{i}∈\mathit{\boldsymbol{P}}$处的局部坐标系$\mathit{\boldsymbol{C}}(p_{i})=\{p_{i}, \mathit{\boldsymbol{x}}_{pi}, \mathit{\boldsymbol{y}}_{pi}, \mathit{\boldsymbol{z}}_{pi}\}$：用$\mathit{\boldsymbol{D}}(p_{i})$的3个右奇异矢量$\mathit{\boldsymbol{v}}_{1}, \mathit{\boldsymbol{v}}_{2}, \mathit{\boldsymbol{v}}_{3}$构造$\mathit{\boldsymbol{C}}(p_{i})$的3个轴，$\mathit{\boldsymbol{z}}$轴方向的单位向量$\mathit{\boldsymbol{z}}_{pi}$需要满足条件

$ \begin{array}{*{20}{c}} {\left| {{\mathit{\boldsymbol{z}}_{{p_i}}} \cdot \mathit{\boldsymbol{N}}\left( {{p_i}} \right)} \right| = \max \left\{ {\left| {{\mathit{\boldsymbol{v}}_1} \cdot \mathit{\boldsymbol{N}}\left( {{p_i}} \right)} \right|,} \right.}\\ {\left. {\left| {{\mathit{\boldsymbol{v}}_2} \cdot \mathit{\boldsymbol{N}}\left( {{p_i}} \right)} \right|,\left| {{\mathit{\boldsymbol{v}}_3} \cdot \mathit{\boldsymbol{N}}\left( {{p_i}} \right)} \right|} \right\}} \end{array} $

(14)

式中，$x$轴方向的单位向量$\mathit{\boldsymbol{x}}_{p_{i}}$对应剩下的奇异值较大的奇异矢量。由于$\mathit{\boldsymbol{z}}_{pi}$和$\mathit{\boldsymbol{x}}_{p_{i}}$对应的奇异向量既可能同向也可能反向，所以通过公式进行约束, 即

$ \left\{ \begin{array}{l} \left( {g\left( {{p_i}} \right) - {p_i}} \right) \cdot {\mathit{\boldsymbol{x}}_{{p_i}}} \ge 0\\ \left( {g\left( {{p_i}} \right) - {p_i}} \right) \cdot {\mathit{\boldsymbol{z}}_{{p_i}}} \ge 0 \end{array} \right. $

(15)

式中，$g(p_{i})$是球面映射$G_{r}(p_{i})$的质心。确定$\mathit{\boldsymbol{z}}_{pi}$和$\mathit{\boldsymbol{x}}_{p_{i}}$的方向后，可以根据$\mathit{\boldsymbol{y}}_{p_{i}}=\mathit{\boldsymbol{z}}_{p_{i}}×\mathit{\boldsymbol{x}}_{p_{i}}$计算得到$y$轴方向的单位矢量。使用该方法可以构造点$p_{i}$处的局部坐标系$\mathit{\boldsymbol{C}}(p_{i})$，同理也可以构造点$q_{j}∈\mathit{\boldsymbol{Q}}$处的局部坐标系。

假设$(p, q)$是一对匹配点，可以通过一个3×3的旋转矩阵$\mathit{\boldsymbol{R}}$和一个3维平移向量$\mathit{\boldsymbol{t}}$将点$q$变换到点$p$，即

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{R}} = \left( {{\mathit{\boldsymbol{x}}_p},{\mathit{\boldsymbol{y}}_p},{\mathit{\boldsymbol{z}}_p}} \right){\left( {{\mathit{\boldsymbol{x}}_q},{\mathit{\boldsymbol{y}}_q},{\mathit{\boldsymbol{z}}_q}} \right)^{ - 1}}\\ \mathit{\boldsymbol{t}} = \mathit{\boldsymbol{p}} - \mathit{\boldsymbol{Rq}} \end{array} \right. $

(16)

式中，$(x_{p}, y_{p}, z_{p})$和$(x_{q}, y_{q}, z_{q})$分别对应的是点$p$和$q$的局部坐标系，$\mathit{\boldsymbol{R}} \in {\bf{R}^{3 \times 3}}, \mathit{\boldsymbol{t}} \in {\bf{R}^3}$。为了实现匹配点对间变换关系的参数化，需要将旋转矩阵表达成$(ω, θ)$的形式(式中，$θ∈[0, {\rm π}]$为旋转角度，$\mathit{\boldsymbol{\omega}}=[ω_{x}, ω_{y}, ω_{z}]^{\rm T}$为旋转轴)。这样，任何一个旋转矩阵可以表示为由$(ω, θ)$张开的空间中以原点为球心、以$π$为半径的球内或表面上的一点。为了由旋转矩阵$\mathit{\boldsymbol{R}}$求出$(ω, θ)$用于霍夫变换，用四元素来表示旋转变换，即

$ {\mathit{\boldsymbol{Q}}_r} = \left( {{q_0},v} \right) = {q_0} + {q_1}\mathit{\boldsymbol{i}} + {q_2}\mathit{\boldsymbol{j}} + {q_3}\mathit{\boldsymbol{k}} $

(17)

式中，$q_0^2 + q_1^2 + q_2^2 + q_3^2 = 1$。由此，旋转矩阵可以表达为

$ \mathit{\boldsymbol{R}} = \left[ {\begin{array}{*{20}{c}} \begin{array}{l} q_0^2 + q_1^2 - q_2^2 - q_3^2\\ 2\left( {{q_1}{q_2} + {q_0}{q_3}} \right)\\ 2\left( {{q_1}{q_3} - {q_0}{q_2}} \right) \end{array}&\begin{array}{l} 2\left( {{q_1}{q_2} - {q_0}{q_3}} \right)\\ q_0^2 - q_1^2 + q_2^2 - q_3^2\\ 2\left( {{q_2}{q_3} + {q_0}{q_1}} \right) \end{array}&\begin{array}{l} 2\left( {{q_1}{q_3} + {q_0}{q_2}} \right)\\ 2\left( {{q_2}{q_3} - {q_0}{q_1}} \right)\\ q_0^2 - q_1^2 - q_2^2 + q_3^2 \end{array} \end{array}} \right] $

(18)

结合式(16)(18)可以得到参数$q_{0}, q_{1}, q_{2}, q_{3}$。另一方面，也可以使用旋转参数来表示四元数^[29]，即

$ {\mathit{\boldsymbol{Q}}_r} = \cos \left( {\frac{\theta }{2}} \right) + \omega \sin \left( {\frac{\theta }{2}} \right) $

(19)

旋转参数$ω$以及$θ$的计算公式为

$ \theta = 2{\cos ^{ - 1}}{q_0},\;\;\;\;\omega = \left\{ \begin{array}{l} \frac{v}{{\sin \left( {\theta /2} \right)}}\;\;\;\;\;\;\theta \ne 0\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(20)

利用上述公式可以得到集合$\mathit{\boldsymbol{Match}}$中所有点对的旋转参数$(ω, θ)$以及平移向量参数$\mathit{\boldsymbol{t}}$。下一步的任务就是研究这些参数在由$(ω, θ)$和$\mathit{\boldsymbol{t}}$张开的空间中的分布情况。为了精确地统计出参数空间中每一位置所聚集的点数，这两个空间都应该被分割为小的单元，通过统计每个单元内点的个数，将参数空间中包含点数最多的单元的坐标作为点云坐标变换的参数。

如前所述，通过对参数空间中不同单元内点的计数来确定点在参数空间中点的密度分布。因此，密度分布的精度取决于单元的大小以及空间划分的均匀性。实际上，需要划分两类参数空间：一类是球体表面的划分(即：相应于旋转轴$ω$)，另一类就是欧几里得空间的划分(即相应于$θ$的1维空间的划分以及相应于$\mathit{\boldsymbol{t}}$的3维空间的划分)。

本文使用一种基于正多面体的方法对球体表面进行划分。这种方式的优点是可以对单元无限细分而不改变它们的形状。而且，在划分过程中该方法可以确保均匀性。与其他多面体相比，正二十面体具有以下优点：在所有已知的正多面体中，它拥有最多面的，它的面都是正三角形，4个细分以后的子三角形具有更相似的形状。基于以上考虑，本文采用与球体内接的正二十面体对球体表面进行划分，并通过细分正二十面体的表面实现对球面的进一步细分(如图 3所示)。图 3(a)给出了初始三角形以及对它们进行细分的示意图。由该图可以看出，在细分过程中，该方法利用三角形各边的中点将其划分成4个子三角形。图 3(b)给出了为子三角形分配索引的方式。

假设将被划分的球体中心与坐标原点$O$一致。可以通过下面3个不等式判断参数空间(此处指的是旋转矩阵)中的点$x$在球面上的投影是否落在$△ABC$内，即

$ \mathit{\boldsymbol{a}} \cdot \mathit{\boldsymbol{x}} > 0,\mathit{\boldsymbol{b}} \cdot \mathit{\boldsymbol{x}} > 0,\mathit{\boldsymbol{c}} \cdot \mathit{\boldsymbol{x}} > 0 $

(21)

式中，$\mathit{\boldsymbol{a}} = \mathit{\boldsymbol{OA}} \times \mathit{\boldsymbol{OB}}, \mathit{\boldsymbol{b}} = \mathit{\boldsymbol{OB}} \times \mathit{\boldsymbol{OC}}, \mathit{\boldsymbol{c}} = \mathit{\boldsymbol{OC}} \times \mathit{\boldsymbol{OA}}$分别是$△ABC3$个顶点与原点$O$连线的向量。只有当式(21)中的3个判别式同时满足的时候$x$的投影才会落在$△ABC$内。如果将向量$\mathit{\boldsymbol{a, b, c}}$结合在一起，形成判别矩阵$\mathit{\boldsymbol{J}}=(\mathit{\boldsymbol{a, b, c}})$，$△ABC$的4个子三角形的判别矩阵$\mathit{\boldsymbol{J}}_{i}, (i=0, 1, 2, 3)$可以通过以下递归方式得到

$ {\mathit{\boldsymbol{J}}_i} = \mathit{\boldsymbol{J}}{\mathit{\boldsymbol{M}}_i}\;\;\;\left( {i = 0,1,2,3} \right) $

(22)

式中$\mathit{\boldsymbol{M}}_{i}, i=0, 1, 2, 3$是4个常量矩阵

$ {\mathit{\boldsymbol{M}}_0} = \left[ {\begin{array}{*{20}{c}} 1&{ - 1}&0\\ 0&1&0\\ 0&{ - 1}&1 \end{array}} \right]\;\;\;\;\;{\mathit{\boldsymbol{M}}_1} = \left[ {\begin{array}{*{20}{c}} 1&0&{ - 1}\\ 0&1&{ - 1}\\ 0&0&1 \end{array}} \right] $

$ {\mathit{\boldsymbol{M}}_2} = \left[ {\begin{array}{*{20}{c}} 1&0&0\\ { - 1}&1&0\\ { - 1}&0&1 \end{array}} \right]\;\;\;\;{\mathit{\boldsymbol{M}}_3} = \left[ {\begin{array}{*{20}{c}} { - 1}&1&1\\ 1&{ - 1}&1\\ 1&1&{ - 1} \end{array}} \right] $

如上所述，本文对球面的划分始于内接正二十面体，并以如图 3所示的递归的方式进行细分，直到满足收敛条件为止。如果细分次数是$m_{1}$次，那么将得到$20×4^{m_{1}}$个三角形。在细分的过程中落入同一个三角形的点会被记录到数组里。

可以将旋转角$θ$和平移向量$\mathit{\boldsymbol{t}}$认为是$k$维欧几里得空间中的点($θ$是$k=1，\mathit{\boldsymbol{t}}$是$k=3$)。要得到这两类点的密度分布需要对$k$维欧几里得空间进行划分。本文利用垂直于每个轴的超平面对单元的连续平分来实现对欧几里得空间的划分。细分前后的单元都是多面体，其$k$个表面均与坐标轴垂直。初始单元由参数点在各坐标轴的分布范围来确定。对初始单元的划分首先沿着第1个坐标轴进行，得到两个子单元，然后沿着第2个坐标轴将两个子单元划分成4个子单元，以此类推。若沿所有的坐标轴均划分$m_{2}$次，则最终得到单元数是$2^{km_{2}}$个。

通过对$θ, ω$和$\mathit{\boldsymbol{t}}$对应的参数空间的划分，不同视角的点云数据间的变换参数估计为具有最多计数的单元的坐标。在整个估算的过程中，先估算旋转角度θ，如果$θ≈0$或者$θ≈π$，对浮动点云数据点施加旋转操作以规避在这两点的不确定性。估算旋转角$θ$后再估计旋转轴$ω$以及平移矩阵$\mathit{\boldsymbol{t}}$。

为了验证本文算法，利用本文算法对来自于斯坦福大学官网的点云数据进行了拼接实验。该点云数据共有35 947个点。将整份点云数据分割为部分重叠的两份点云数据，图 4给出了两份部分重叠的点云数据，左图是目标点云数据$\mathit{\boldsymbol{P}}$(25 141点)，右图是经过任意的旋转和平移操作后得到的浮动点云数据$\mathit{\boldsymbol{Q}}$(18 475点)，所以重叠率为22%。

如前所述，在进行微分信息计算之前，法向量需要调整为一致向外。图 5给出了对图 4中目标点云$\mathit{\boldsymbol{P}}$和浮动点云$\mathit{\boldsymbol{Q}}$的法向量进行调整前后的对比图。图 5(a)(b)分别给出了点云数据$\mathit{\boldsymbol{P}}$的法向量在调整前、后的状态；图 5(c)(d)分别给出了点云数据$\mathit{\boldsymbol{Q}}$的法向量在调整前、后的状态。从图 5(a)(c)可以看到调整前两片点云的法向量的方向是不一致的，相反从图 5 (b)(d)可以看出，利用本文算法调整后点云的法向量基本一致朝外了。这个结果说明本文提出的广度优先结合${\rm cost}$权值的方法可以有效地实现法向量方向的调整。

由曲率相似度确定的初始匹配点对的筛选是基于刚体不变特征的。因此，该算法的性能取决于不同点云的特征值的不变性。图 6给出的是两片点云数据的刚体不变特征和曲率特征的示意图，图 6(a)给出的是两片不同点云的刚体不变量特征值的分布情况，其中横轴和纵轴分别为$σ_{1}/σ_{3}$，$σ_{2}/σ_{3}$，蓝色圆圈代表的是目标点云的刚体不变量特征，红色的点代表的是浮动点云的刚体不变量特征。图 6(b)给出了两片不同点云的曲率特征值的分布情况。其中，横轴和纵轴分别为高斯曲率和最大主曲率的绝对值，蓝色的圆圈代表目标点云的曲率特征，红色的点代表浮动点云的曲率特征。图 6(a)(b)两图中的左图均为右图的局部放大，从局部放大图可以明显看出，有部分点的不变特征量标识完全重合，基本上没有两点的曲率特征是完全一致的，所以本文使用刚体不变量标识为特征进行筛选将具有更强的鲁棒性，会获得更好的效果。

由于基于刚体不变特征得到的匹配点对中还存在大量错误的匹配点对，为了进一步提高匹配点对的正确率，利用了PFH算法对匹配点对进行进一步的剔除。为了验证该方法的合理性，图 7给出了$\mathit{\boldsymbol{P}}$和$\mathit{\boldsymbol{Q}}$经过PFH过滤后得到的3组正确的匹配点对及其相应的PFH示意图。图 7(a)给出的是匹配点对的对应关系示意图，图 7(b)(c)(d)给出的分别是匹配点对的PFH特征值，图中褐红色表示点云$\mathit{\boldsymbol{P}}$，浅蓝色表示点云$\mathit{\boldsymbol{Q}}$，横轴表示的是bin序数，纵轴表示的是每个bin中点数的比例。由该图可以看出正确匹配点对的PFH特征值是近似相同的。

为了验证本文算法对噪声的鲁棒性，对目标点云中点的3个坐标值分别添加3组相互独立，且均值为零，标准方差为点平均间距的1%5%高斯噪声。图 8给出了在曲率相似度匹配基础上经过刚体不变量筛选(蓝线)以及在此基础上进一步经过PFH筛选(红线)后的正确率随噪声大小的变化曲线。这里，匹配正确率定义为$A = \frac{{{N_1}}}{N}$，其中，$N$为匹配点对数量，$N_{1}$为正确的匹配点对数量。由该图可以看出使用PFH作为补充性剔除算法可以进一步增加算法鲁棒性，显著提升匹配点对的正确率，平均可以提升2倍左右。

经过上面的PFH过滤后得到了一组新的匹配点对，下面需要使用霍夫变换，求取两片点云数据间的旋转和平移参数。为了利用霍夫变换进行参数检测，需要对参数空间进行划分。本文将旋转角$θ$划分成256个子单元，对旋转轴$ω$的球面共进行了6次球面递归划分，得到81 920个子三角形。对平移向量$\mathit{\boldsymbol{t}}$的立方体沿着每个坐标轴均匀划分4次，得到4 096个子立方体。基于以上划分，图 9 (a)给出了旋转参数($ω, θ$)在以原点为中心点，半径为${\rm π}$的球体内和球面上的分布情况，而9(b)给出的则是平移参数$\mathit{\boldsymbol{t}}$在3维欧几里得空间中的分布情况。由该图可以看出, 两个空间内都存在明显的高密度区，这与本文前面的推断也是吻合的。参数空间中这些点的聚集不仅仅表明了本文所提出的基于曲率、刚体不变量特征值、PFH匹配点对提取方法的优越性，同时也验证了本文所提出的基于霍夫变换的参数检测策略的合理性。

为了研究霍夫变换对提升算法抗噪性能的效果，将其与经典的Ransac^[30]算法进行了比较。其中，对球面空间的匹配参数设定随机采样次数$k$=1 000，点到空间直线的距离阈值设为0.08。利用Ransac算法在不同的噪声条件下各运行了10次，其正确率如表 1所示。由表 1可以看出在不同噪声的情况下，使用霍夫变换得到的匹配点对正确率更高，还可以看出Ransac算法各次运行的结果各不相同的，重复性较差，这是因为Ransac算法具有一定的随机性。而从图 10可以看出, 霍夫变换除了具有更高的正确率外，抗噪能力也是强于Ransac算法，当噪声达到3%以后，Ransac的正确率急速下降，而霍夫变换在噪声5%的时候依然还保持着50%的正确率。表 2的第1行是Ransac算法10次运行时间的平均值，从表 2可以看出霍夫变换的运行时间也远远少于Ransac方法。综上所述，Hough在正确率、稳定性、抗噪性、运行速度上都要优于Ransac方法。

图 11给出了不同方法对含2%噪声的点云数据进行拼接的结果。显然使用霍夫变换的方法能得到更佳的拼合结果。从图 10和图 11也可以看出, 参数聚类在点云拼接中的重要性和有效性。

提出了一种针对散乱点云数据的拼合算法。该算法首先利用曲率相似度建立初始匹配点对。然后，分别利用刚体不变特征以及持续特征直方图(PFH)对初始匹配点对进行连续两次筛选形成新的匹配点对。通过基于PFH的二次筛选，匹配点对的正确率提升2倍左右。为了避免误匹配对点云拼接的影响，本文提出对坐标变换进行参数化处理并利用霍夫变换对坐标变换参数进行最优估计，从而最终实现点云的拼合。实验结果表明, 该方法可以应用在对任何部分重叠或者全部重叠的具有任意相对位置的点云数据的拼接中，并且保持着较高的精度和较强的鲁棒性。由于本文算法需要分步实现，比较耗时，下一步将研究点云数据精简策略及其相应的拼合方法，以进一步提高算法的效率。