点云配准算法研究已久，其应用领域涉及曲面匹配^[1-2]、目标识别^[3-4]、姿态估计^[5]以及文物复原^[6]等多个方面。颅骨配准是点云配准算法的一个重要应用领域，是颅面复原^[7]的一个重要步骤。颅骨配准就是从已有的颅骨数据库中找出一个或多个与待复原颅骨U最为相似的颅骨S的过程。找到的颅骨S的面貌即可作为待复原颅骨U的参考面貌，从而为颅面复原提供可能的依据。

对于颅骨的点云数据模型，可以采用点云配准的方法来实现颅骨配准。点云配准方法中应用最多的是基于特征的配准方法，包括基于全局特征的配准方法^[8]和基于局部特征的配准方法^[9-10]两大类。全局特征描述的是整个点云模型，而局部特征只描述特征点的邻域特征，与全局特征相比，局部特征更适用于部分覆盖的复杂点云的配准。由于颅骨数据模型复杂，且不同个体的颅骨差异很小，因此对其配准精度的要求较高，一般采用基于局部特征的由粗到细的方法来实现颅骨的精确配准。

粗配准就是将两个位于不同坐标系中的颅骨进行粗略对齐的过程，通常采用点云的特征实现，如法向和曲率^[11]特征以及由积分不变量^[12]计算的或凹或凸的特征区域等。细配准就是在粗配准的基础上，将两个颅骨进行进一步细对齐的过程。目前应用最为广泛的点云细配准算法是由Besl等人^[13]提出的最近点迭代(ICP)算法。该算法步骤简单，易于实现，但要求两个待配准的点云之间存在包含关系。鉴于此，国内外学者提出了许多改进的ICP算法，如王欣等人^[14]提出了基于点云边界特征点的改进ICP算法，提高了逆向工程中点云数据配准的效率和精度；Li等人^[15]提出了一种基于动态调整因子的ICP算法，在不影响配准精度和收敛方向的情况下，可以大大提高算法的配准速度；Mavridis等人^[16]提出了一种基于混合优化系统的稀疏ICP算法，提高了点云配准的精度和速度；Du等人^[17]提出了尺度ICP算法，解决了含尺度因素的点云配准问题。

以上这些算法在点云的配准精度、速度以及尺度因素等方面有了一定程度的改进，但是对颅骨数据的配准结果却并不十分理想。针对颅骨点云配准中的配准精度和收敛速度的问题，提出一种由粗到细的优化配准方法。首先采用基于局部特征的点云配准算法实现颅骨粗配准，然后通过在ICP算法中加入高斯概率模型和动态迭代系数的方式来改进ICP算法，并将其应用到颅骨点云模型的细配准中，不仅可以解决由噪声引起的配准效果不佳的问题，而且可以大大提高算法的收敛速度，从而实现颅骨的快速精确配准。

局部特征是指一种快速、鲁棒的点云特征描述子，这里包括局部深度、法线的偏角和点云密度等3个特征，它们都具有旋转和平移不变的特性。

设点云$\boldsymbol{P} $由$\left\{ {{\boldsymbol{p}_1}, {\boldsymbol{p}_2}, \cdots, {\boldsymbol{p}_N}} \right\} $等 $N $个点组成，其中任意一点${\boldsymbol{p}_i} $的邻域定义为以${\boldsymbol{p}_i} $为球心$r $为半径的球体，记作$\boldsymbol{p}_n^i = \left\{ {\boldsymbol{p}_i^j|j = 1, 2, \cdots, k} \right\} $。令${\boldsymbol{n}_i} $和$\boldsymbol{n}_i^j $分别表示点${\boldsymbol{p}_i} $和$\boldsymbol{p}_i^j $的法向量，则设置${\boldsymbol{n}_i} $为局部坐标轴。沿局部坐标轴的正方向，定义2维平面$\boldsymbol{L} $为投影平面，如图 1(a)所示。将${\boldsymbol{p}_i} $邻域内的所有点投影到$\boldsymbol{L} $上，则得到一个新的点集$\boldsymbol{p'}_n^i = \left\{ {\boldsymbol{p'}_i^j|j = 1, 2, \cdots, k} \right\}$。定义$\boldsymbol{p}_i^j $和$\boldsymbol{p'}_i^j $之间的距离为局部深度，那么局部深度的取值范围应该为[0，2 $r $]。局部深度的定义式为

$ {d_j} = r-{\boldsymbol{n}_i} \cdot \left( {\boldsymbol{p}_i^j-{\boldsymbol{p}_i}} \right) $

(1)

计算法线的偏角时，首先采用基于局部最小二乘估计的方法计算点的法线。对于点云$\boldsymbol{P} $中任意一点${\boldsymbol{p}_i} $，设其邻域$\boldsymbol{p}_n^i $的质心为${\boldsymbol{\bar p}_i} $，那么${\boldsymbol{p}_i} $的协方差矩阵为

$ Cov\left( {{\mathit{\boldsymbol{p}}_i}} \right) = {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{p}}_i^1 - {{\mathit{\boldsymbol{\bar p}}}_i}}\\ \vdots \\ {\mathit{\boldsymbol{p}}_i^k - {{\mathit{\boldsymbol{\bar p}}}_i}} \end{array}} \right]^{\rm{T}}} \cdot \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{p}}_i^1 - {{\mathit{\boldsymbol{\bar p}}}_i}}\\ \vdots \\ {\mathit{\boldsymbol{p}}_i^k - {{\mathit{\boldsymbol{\bar p}}}_i}} \end{array}} \right] $

(2)

将$Cov\left( {{\boldsymbol{p}_i}} \right) $的最小特征值作为${\boldsymbol{p}_i} $的特征向量，如图 1(b)所示，那么法线的偏角可表示为

$ {\theta _j} = \arccos \left( {{\mathit{\boldsymbol{n}}_i} \cdot \mathit{\boldsymbol{n}}_i^j} \right) $

(3)

式中，${\theta _j} \in [0, \pi] $。

将3D点投影到2D平面上是一种有效的描述3D点云局部几何形状的方法，可描述点的分布信息，如图 1(c)所示。首先在2D投影平面$\boldsymbol{L} $和${\boldsymbol{n}_i} $的交点处绘制一个半径为r的圆；然后将圆用相同距离的宽度分成几个环；最后通过计算落在每个环内的点的数目来生成直方图。

分布在每个环内的点的密度用点${{\boldsymbol{p'}}_i}\; $到$\boldsymbol{p'}_i^j $的距离来计算，即水平投影距离$\rho $，其计算式为

$ {\rho _j} = \sqrt {{{\left\| {{{\mathit{\boldsymbol{p'}}}_i} - \mathit{\boldsymbol{p}}_i^{'j}} \right\|}^2} - {\mathit{\boldsymbol{n}}_i} \cdot {{\left. {\left( {{{\mathit{\boldsymbol{p'}}}_i} - \mathit{\boldsymbol{p}}_i^{'j}} \right)} \right)}^2}} $

(4)

式中，${\left\| {{{\boldsymbol{p'}}_i}-\boldsymbol{p'}_i^j} \right\|^2} $表示${\boldsymbol{p}_i} $和$\boldsymbol{p}_i^j $之间的距离，${\rho _j} $的取值范围为$[0, r] $。

在计算了点的3个局部几何特征后，即得到3个子直方图，将这3个子直方图合成为一个直方图，就得到了局部特征的描述子。

粗配准采用上面介绍的局部特征描述子来实现，主要包含局部特征提取、相关性计算和外点删减等3个步骤。

通常初始的颅骨点云数据模型的数据量都很大，使得其计算的时间复杂度和空间复杂度都较高。因此首先采用实时压缩算法对颅骨的源点云${\boldsymbol{P}_S} $和目标点云${\boldsymbol{P}_T} $进行均匀采样，得到两个相应的点集：$\boldsymbol{M} = \left\{ {{\boldsymbol{m}_i}|i = 1, 2, \cdots, {N_M}} \right\} $和$\boldsymbol{D} = \left\{ {{\boldsymbol{d}_i}|i = 1, 2, \cdots, {N_D}} \right\} $，${N_M} $和$ {N_D}$分别表示点集$\boldsymbol{M} $和$\boldsymbol{D} $的大小。然后计算点集$\boldsymbol{M} $和$\boldsymbol{D} $的局部特征描述子，分别记为${\boldsymbol{F}_M} = \left\{ {\boldsymbol{f}_M^i|i = 1, 2, \cdots, {N_M}} \right\} $和${\boldsymbol{F}_D} = \left\{ {\boldsymbol{f}_D^i|i = 1, 2, \cdots, {N_D}} \right\} $。

两个点集$\boldsymbol{M} $和$\boldsymbol{D} $的刚体变换可以通过计算点到点的相关性来求解。这里通过求解特征集${\boldsymbol{F}_M} $和${\boldsymbol{F}_D} $的相关性来计算刚体变换。

对于特征$\boldsymbol{f}_M^i \in {\boldsymbol{F}_M} $，在${\boldsymbol{F}_D} $中寻找其$ k$最近特征为$\left\{ {\boldsymbol{f}_D^{i1}, \boldsymbol{f}_D^{i2}, \cdots, \boldsymbol{f}_D^{ik}} \right\} $，那么其相关点$\left\{ {{\boldsymbol{d}_{i1}}, {\boldsymbol{d}_{i2}}, \cdots, {\boldsymbol{d}_{ik}}} \right\} $就是点$\boldsymbol{p}_M^i $的$k $相关候选点，记作${\boldsymbol{e}_i} = \left\{ {\left( {{\boldsymbol{m}_i}, {\boldsymbol{d}_{ij}}} \right)|j = 1, 2, \cdots, k} \right\} $。那么最终得到点集$\boldsymbol{M} $和$\boldsymbol{D} $的相关候选点集为$\boldsymbol{E} = \left\{ {{e_i}|i = 1, 2, \cdots, {N_M}} \right\} $，$\boldsymbol{E} $中含$k \cdot {N_M} $个相关点。

目前已经建立了点集$\boldsymbol{M} $和$\boldsymbol{D} $的相关候选点集$\boldsymbol{E} $。由于待配准的点云大多是部分覆盖，所以并不是$\boldsymbol{M} $中的所有的点都能在$\boldsymbol{D} $中找到相关点，因此需要从$\boldsymbol{E} $中删除外点。理想情况下，两个点集中正确配准的点的距离为0，但是实际并非如此，因此设置阈值${d_f} $来移除外点。$ {d_f}$的计算式为

$ {d_f} = {\mu _f} - \alpha \cdot {\sigma _f} $

(5)

式中，${\mu _f} $表示均值，${\sigma _f} $表示集合$\boldsymbol{E} $中点对的局部特征描述子间的欧几里德距离的标准偏差；$\alpha $用于控制输出的相关点集的大小，是个常数，通过多次实验的建议取值范围为$\alpha $∈[0.5，1.5]。

通过删除外点就得到了相关性点集$\boldsymbol{C} = \left\{ {\left( {\boldsymbol{c}_M^i, \boldsymbol{c}_D^i} \right)|\boldsymbol{c}_M^i \in {\boldsymbol{P}_M}, \boldsymbol{c}_T^i \in {\boldsymbol{P}_D}, i = 1, 2, \cdots, l} \right\} $，$\boldsymbol{C} $包含$l $对相关点，由此实现了点集$\boldsymbol{M} $和$\boldsymbol{D} $粗配准。

细配准阶段采用改进的ICP算法实现，首先在ICP算法中加入高斯概率模型^[18]以提高算法的抗噪性，然后加入动态迭代系数，在不影响配准精度和收敛方向的情况下，提高算法的收敛速度。

对于含噪声的点集$\boldsymbol{M} = \{ {\boldsymbol{m}_i}\} _{i = 1}^{{N_M}} $和$\boldsymbol{D} = \{ {\boldsymbol{d}_j}\} _{j = 1}^{{N_D}} $，可将其配准问题转化为概率密度估计的问题。根据全概率定理，点集$\boldsymbol{M} $的概率公式可表示为

$ p\left( \mathit{\boldsymbol{m}} \right) = \sum\limits_{j = 1}^{{N_y}} {p\left( {{\mathit{\boldsymbol{d}}_j}} \right)p\left( {\mathit{\boldsymbol{m}}\left| {{\mathit{\boldsymbol{d}}_j}} \right.} \right)} $

(6)

若点集$\boldsymbol{M} $和$\boldsymbol{D} $都满足高斯概率模型，那么它们之间的关系可描述为

$ p\left( {\mathit{\boldsymbol{m}}\left| \mathit{\boldsymbol{d}} \right.} \right) = \frac{1}{{{{\left( {2{\rm{ \mathsf{ π} }}{\sigma ^2}} \right)}^{n/2}}}}{\exp ^{ - \frac{{\left\| {\mathit{\boldsymbol{T}}\left( \mathit{\boldsymbol{m}} \right) - \mathit{\boldsymbol{d}}} \right\|_2^2}}{{2{\sigma ^2}}}}} $

(7)

式中，$n $表示点集的维数，${\sigma ^2} $表示高斯概率模型的方差，$\boldsymbol{T}\left( \boldsymbol{m} \right) $表示$\boldsymbol{m} $的刚体变换，$\left\| {\boldsymbol{T}\left( \boldsymbol{m} \right)-\boldsymbol{d}} \right\|_2^2 $表示点集$\boldsymbol{T}\left( \boldsymbol{m} \right) $和$\boldsymbol{d} $之间的距离。

通常3维扫描仪获取的点云数据都是杂乱无章的，需要对其进行模型简化，使其在数据量上适当降低并满足均匀分布。因此对于点集$\boldsymbol{D} $，有

$ p\left( \mathit{\boldsymbol{d}} \right) = \frac{1}{{{N_D}}} $

(8)

目标函数定义为

$ \begin{array}{*{20}{c}} {F\left( {\mathit{\boldsymbol{m}};\mathit{\boldsymbol{R}},\mathit{\boldsymbol{t}},{\sigma ^2}} \right) = - \sum\limits_{i = 1}^{{N_M}} {\log p\left( {{\mathit{\boldsymbol{m}}_i}\left| {\mathit{\boldsymbol{R}},\mathit{\boldsymbol{t}},{\sigma ^2}} \right.} \right)} = }\\ { - \sum\limits_{i = 1}^{{N_M}} {\log } \sum\limits_{j = 1}^{{N_D}} {p\left( {{\mathit{\boldsymbol{m}}_j}} \right)p\left( {{\mathit{\boldsymbol{m}}_i}\left| {{\mathit{\boldsymbol{d}}_j}\mathit{\boldsymbol{R}},\mathit{\boldsymbol{t}},{\sigma ^2}} \right.} \right)} } \end{array} $

(9)

首先建立刚体变换的新目标函数为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}\left( {{{\mathit{\boldsymbol{\tilde R}}}_k},{{\mathit{\boldsymbol{\tilde t}}}_k},{{\tilde \sigma }_k}} \right) = \sum\limits_{i = 1}^{{N_M}} {\sum\limits_{j = 1}^{{N_D}} {p\left( {{\mathit{\boldsymbol{d}}_j}\left| {{\mathit{\boldsymbol{m}}_i},{\mathit{\boldsymbol{R}}_{k - 1}},{\mathit{\boldsymbol{t}}_{k - 1}},\sigma _{k - 1}^2} \right.} \right)} } \times }\\ {\log \left( {p\left( {{\mathit{\boldsymbol{d}}_j}} \right)p\left( {{\mathit{\boldsymbol{m}}_i}\left| {{\mathit{\boldsymbol{d}}_j},{{\mathit{\boldsymbol{\tilde R}}}_k},{{\mathit{\boldsymbol{\tilde t}}}_k},\tilde \sigma _k^2} \right.} \right)} \right)} \end{array} $

(10)

令$p_i^{k-1} = p\left( {{\boldsymbol{d}_{{c_k}\left( i \right)}}|{\boldsymbol{m}_i}, {{\boldsymbol{\tilde R}}_k}, {{\boldsymbol{\tilde t}}_k}, \tilde \sigma _k^2} \right) $，则式(10) 可简化为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{M}}\left( {{{\mathit{\boldsymbol{\tilde R}}}_k},{{\mathit{\boldsymbol{\tilde t}}}_k},{{\tilde \sigma }_k}} \right) = \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1} \times \frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2}}{{2{{\tilde \sigma }^2}}} - } }\\ {\sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times \log \left( {{N_D}{{\left( {2{\rm{ \mathsf{ π} }}\tilde \sigma _k^2} \right)}^{\frac{n}{2}}}} \right)} \end{array} $

(11)

式中，只有$-\sum\limits_{i = 1}^{{N_M}} {\boldsymbol{p}_i^{k-1}} \times \frac{{\left\| {{{\boldsymbol{\tilde R}}_k}{\boldsymbol{m}_i} + {{\boldsymbol{\tilde t}}_k}-{\boldsymbol{d}_{{c_k}\left( i \right)}}} \right\|_2^2}}{{2{{\tilde \sigma }^2}}} $与旋转矩阵和平移矢量有关，因此最大化式(11) 的问题就等价于最大化式(11) 的第1部分的问题，因此求解目标进一步转换为

$ \begin{array}{*{20}{c}} {\left( {{\mathit{\boldsymbol{R}}_k},{\mathit{\boldsymbol{t}}_k}} \right) = \mathop {\arg \min }\limits_{\tilde R_k^{\rm{T}}{{\tilde R}_k} = {I_n},\det \left( {\mathit{\boldsymbol{\tilde R}}} \right) = 1,{{\tilde t}_k}} \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1}} \times }\\ {\frac{{\left\| {{{\mathit{\boldsymbol{\tilde R}}}_k}{\mathit{\boldsymbol{m}}_i} + {{\mathit{\boldsymbol{\tilde t}}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2}}{{2{{\tilde \sigma }^2}}}} \end{array} $

(12)

定义均方根误差为

$ RMS = {\left( {\sum\limits_{i = 1}^{{N_M}} {{\mathit{\boldsymbol{p}}_i} \times \left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2} } \right)^{\frac{1}{2}}} $

(13)

定义方差为

$ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2 = \sum\limits_{i = 1}^{{N_M}} {\mathit{\boldsymbol{p}}_i^{k - 1} \times \frac{{\left\| {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k} - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2}}{n}} $

(14)

那么高斯概率方差的更新公式为

$ \sigma _k^2 = \left\{ \begin{array}{l} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k - 1}^2/\lambda \;\;\;\;\;\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k - 1}^2/\lambda > \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2\\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2\;\;\;\;\;\;\;\;\;\;\;\;\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _{k - 1}^2/\lambda < \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}} \over \sigma } _k^2 \end{array} \right. $

(15)

式中，$\lambda $称作退火系数。当$\lambda = 1 $时，方差不会更新，即ICP算法。当$\lambda $取值较大时，可以有效地提高迭代速度，但容易陷入局部最小值。通过多次实验，一般建议$\lambda $的取值在1 2之间。

因此后验高斯概率的计算公式为

$ \begin{array}{*{20}{c}} {p\left( {{\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}\left| {{\mathit{\boldsymbol{m}}_i},{\mathit{\boldsymbol{R}}_k},{\mathit{\boldsymbol{t}}_k},\sigma _k^2} \right.} \right) = }\\ {\frac{{1/{{\left( {2{\rm{ \mathsf{ π} }}\sigma _k^2} \right)}^{\frac{n}{2}}}{{\exp }^{ - \frac{{\left\| {\left( {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k}} \right) - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2}}{{2\sigma _k^2}}}}}}{{\sum\limits_{i = 1}^{{N_M}} {1/{{\left( {2{\rm{ \mathsf{ π} }}\sigma _k^2} \right)}^{\frac{n}{2}}}{{\exp }^{ - \frac{{\left\| {\left( {{\mathit{\boldsymbol{R}}_k}{\mathit{\boldsymbol{m}}_i} + {\mathit{\boldsymbol{t}}_k}} \right) - {\mathit{\boldsymbol{d}}_{{c_k}\left( i \right)}}} \right\|_2^2}}{{2\sigma _k^2}}}}} }}} \end{array} $

(16)

动态迭代系数是一个可以自动调整刚体变换参数的整数值，可以在不影响ICP算法的配准精度和收敛方向的情况下，大大提高算法的收敛速度。

在ICP算法中加入动态迭代系数$h $的步骤如下：

1) 计算刚体变换矢量$\boldsymbol{q} = {\left[{R|t} \right]^{\rm{T}}} $以及${\boldsymbol{q}_{k-1}} $的相邻两次迭代的变化量$\Delta {\boldsymbol{q}_{k-1}} $；

2) 用刚体变换矢量$\Delta {\boldsymbol{q}_{k-1}} $更新基本ICP算法中的${\boldsymbol{M}_k} $共$h $次，即执行${\boldsymbol{M}_k} = \Delta {\boldsymbol{q}_{k-1}}\left( {{M_{k-1}}} \right) $共$h $次。

通常，动态迭代系数$h $取大于等于0的整数。当$h $不同时，算法的收敛速度也不同，如图 2所示。随着动态迭代系数$h $的增大，收敛速度也会越来越快，但是当$h $增大到一定的程度，收敛曲线会出现震荡，可能不再收敛。通过多次实验验证，建议$h $的取值在1 4之间。

由3.1节和3.2节的方法描述，改进的ICP算法可描述如下：

1) 给定刚体变换初值$\boldsymbol{q} = {\left[{{\boldsymbol{R}_0}, {\boldsymbol{t}_0}} \right]^{\rm{T}}} $，${\boldsymbol{R}_0} $为初始旋转矩阵，${\boldsymbol{t}_0} $为初始平移矢量；

2) 令${\boldsymbol{M}_0} = {\boldsymbol{R}_0}\boldsymbol{M} + {\boldsymbol{t}_0} $，退火系数$\lambda \in \left( {1, 2} \right] $，初始概率为$p_i^0 = \frac{1}{{{N_x}}}, i = 1, 2, \cdots, {N_x} $，迭代次数$k = 0 $，动态迭代系数$h = 0 $；

3) 令$k = k + 1 $，采用ICP算法^[13]建立点集的相关性${c_k}\left( i \right) \in \left\{ {1, 2, \cdots, {N_D}} \right\} $；

4) 利用奇异值分解(SVD)法计算旋转矩阵${\boldsymbol{R}_k} $和平移矢量${\boldsymbol{t}_k} $；

5) 计算平移矢量${\boldsymbol{t}_k} = \frac{1}{{{N_M}}}\left( {-{\boldsymbol{R}_k}\sum\limits_{i = 1}^{{N_M}} {{m_i} + \sum\limits_{i = 1}^{{N_M}} {{\boldsymbol{d}_{{c_k}\left( i \right)}}} } } \right) $，以及${\boldsymbol{q}_k} $相邻两次迭代的变化量$\Delta {\boldsymbol{q}_k} $；

6) 计算${\boldsymbol{M}_k} = {\boldsymbol{R}_k}\boldsymbol{M} + {\boldsymbol{t}_k} $；

7) 判断均方根误差$RMS $，若$RM{S_k}-RM{S_{k-1}} > 0 $，则执行$h = h + 1 $操作，否则执行$h = 0 $操作；

8) 判断动态迭代系数$h $，若$h > 0 $，则执行${M_k} = \Delta {q_k}\left( {{M_k}} \right) $共$h $次来更新点集${M_k} $；

9) 计算两组点集的方差，并更新高斯概率的方差；

10) 计算后验高斯概率，直到满足$\left| {RM{S_k}-RM{S_{k-1}}} \right| < {\varepsilon _2} $或$k > Ste{p_{\max }} $算法终止，否则转到步骤3)，这里${\varepsilon _2} $和$Ste{p_{\max }} $是预先设置的阈值。

实验1采用的点云数据源于Stanford 3D Scanning Repository^[19]，如图 3所示。图 3(a)为初始的兔子点云模型，图 3(b)为初始的龙点云模型。首先通过局部特征提取、相关性计算和外点删减等步骤实现粗配准，粗配准结果如图 4(a)所示；然后再采用提出的改进ICP算法实现点云细配准，细配准结果如图 4(b)所示。通过反复实验测试，在细配准阶段动态迭代系数$h = 3 $时的配准结果最佳。

从图 4的配准结果可见，提出的基于局部特征的点云配准算法具有良好的配准效果。

为了进一步验证该算法的性能，在细配准阶段，再分别采用ICP算法和PICP算法^[18]进行细配准。PICP算法是一种抗噪性较强的点云配准算法，具有较高的配准精度。ICP算法、PICP算法和提出的改进ICP算法的配准误差、迭代次数和耗时等参数如表 1所示。

从表 1的配准结果可见，ICP算法、PICP算法和改进的ICP算法均能实现点云的细配准，但是PICP算法和改进的ICP算法的配准精度明显要高、耗时明显要短。特别是改进的ICP算法，其配准精度跟PCIP算法不相上下，比ICP算法提高了约30 %；配准速度比ICP算法和PICP算法分别提高了约60 %和30 %。因此说，提出的基于局部区域的点云配准算法是一种精度更高、速度更快的点云配准算法。

实验2采用西北大学可视化技术研究所采集的308套完整的CT扫描的颅骨点云数据模型。图 5为两个待配准的颅骨，图 5(a)为一个未知颅骨U，图 5(b)为一个参考颅骨S。对于这两个待配准颅骨，首先采用基于局部特征的点云配准算法实现颅骨的粗配准，粗配准结果如图 6(a)所示；然后采用提出的改进ICP算法实现颅骨的细配准，细配准结果如图 6(b)所示。在细配准阶段，通过反复实验验证，动态迭代系数$h = 3 $时的配准效果最佳。

从图 6的配准结果来看，基于特征区域的粗配准算法可以将两个颅骨初步对齐，改进的ICP算法则实现了两个颅骨的精确细配准。为了进一步验证改进ICP算法的性能，细配准过程再分别采用ICP算法和PICP算法^[18]实现。ICP算法、PICP算法和改进ICP算法的配准误差、迭代次数以及耗时等参数如表 2所示。

从表 2的配准结果来看，改进ICP算法的配准精度和收敛速度都较高。跟ICP算法相比，改进ICP算法的配准精度和收敛速度分别提高了约30 %和60 %；跟PICP算法相比，改进ICP算法的配准速度提高了约50 %，配准精度不相上下。因此说改进ICP算法一种精度更高、速度更快的颅骨点云模型细配准算法，提出的基于局部特征的点云配准算法是一种有效的颅骨配准算法。

颅骨配准是计算机辅助虚拟颅面复原技术的重要研究内容之一，其配准结果的准确与否将直接影响到颅面复原的准确性。针对颅骨配准中的配准精度和收敛速度的问题，提出了一种颅骨点云模型的局部特征配准方法。首先采用基于局部特征的配准方法实现颅骨粗配准，然后再采用改进的ICP算法实现颅骨细配准，其细配准阶段的配准精度和收敛速度比ICP算法分别提高了约30 %和60 %，实现了颅骨的最终精确配准。该算法不仅适用于公共点云数据模型的配准，而且对颅骨点云模型的配准效果尤为精确。但是，算法设计中未考虑到尺度因素对配准结果的影响。在今后的研究中，要综合多种因素，进一步优化颅骨配准算法，并利用其配准结果实现颅骨面貌复原以及颅骨软组织厚度的评价方法，实现颅骨面貌的虚拟复原研究。