遗存颅骨身份认定在医学、灾难调查、古人类考古等^[1-3]方面有重要的应用，一直以来是国内外研究热点问题。颅面复原是颅骨身份认定的重要途径之一。基于知识的颅面复原是目前颅面复原的一种主要方法，该方法的基本思路是：对待复原颅骨U，从已有颅面数据库中找出与待复原颅骨最为相似的颅骨S，颅骨S的面貌即可作为待复原颅骨U的参考面貌；再通过颅骨U和S之间的特征差别，对颅骨S已有的面皮进行统计形变，最终形成待复原颅骨U的面貌。从数据库中寻找与待复原颅骨U相似颅骨的过程即3维颅骨配准。3维颅骨配准属于3维点云配准的一个应用领域。随着海量3维模型的产生，3维点云配准已经成为研究热点。3维点云配准是将两个不同坐标系中的3维空间数据集转换到统一坐标系统中的计算过程。配准的任务就是求解在不同视角下3维数据坐标点之间的转换关系^[4]。按照3维点云之间的几何变换关系将配准算法分为基于刚体变换的配准方法和基于非刚体变换的配准方法。本文只针对刚体变换的配准方法。针对3维刚体模型的配准问题，一般分为粗配准和精配准两个过程。

粗配准过程就是使得两个位于不同坐标系中的3维点云的相同部位能够大致的匹配和重合。在粗配准过程中，模型的特征描述是匹配中的关键问题。法向和曲率^[5]是常用的局部特征描述方法，具有旋转平移不变性，但是离散曲率对噪声敏感，会带来一定的误差。而积分不变量^[6]是一种描述曲面内在性质的几何不变量，能够避免对噪声敏感的缺点，能准确地描述曲面的局部特征，因而得到了较为广泛的应用，常用于初始粗对齐阶段的特征点或特征区域的提取，很多研究者已经成功地将积分不变量应用到了3维物体配准算法中。

在细配准阶段，运用最广泛的是由Besl等人^[7]提出的最近点迭代 (ICP) 算法。ICP算法步骤简单，易于实现，但是要求两个待匹配的点集间存在着包含关系。若不满足该条件，将会影响ICP的收敛结果，甚至产生错误的匹配。此外，ICP算法对两组点云的相对初始位置要求较高，即两组点云间的初始位置不能相差太大，否则效率低下，并有可能会产生错误的配准结果。鉴于此，国内外的学者提出了很多改进的ICP算法。文献[8]提出了一种分数ICP算法 (FICP)，有效地剔除了异常点，并提高了ICP算法的鲁棒性。文献[9]提出了一种改进的k-d树遍历方法来加快最近点的寻找过程，有效提高了算法的收敛速度。文献[10]基于点云的边界特征点，提出了一种改进的ICP算法，有效地提高了算法的效率和精度。至今为止，ICP算法仍旧是配准算法中研究的重点^[11-12]。

颅骨是非常复杂的3维点云数据，且不同人的颅骨差异非常小，因此其配准的精度要求较高。目前颅骨配准方法大多是采用特征点标定的方法^[13-14]。由于颅骨模型的结构复杂，凹面空洞较多，3维颅骨模型的特征点自动化标定及匹配一直是研究难点，国内外在这方面的研究相对较少, 点的自动标定还不够成熟，大部分都是2维图形特征点标定^[15-16]，而且手动标定和半自动标定特征点有操作者的主观因素，影响配准结果的有效性。

本文提出一种基于特征区域和改进ICP算法3维颅骨匹配算法：1) 根据积分不变量计算颅骨所有点的凹凸性，并将点分为三大类，即凹顶点、凸顶点和平面顶点；2) 根据连通性对凸凹顶点进行K-means聚类，将聚为一类的顶点记为凸的或凹的特征区域；3) 运用主成分分析法计算特征区域相似性，并采取穷举法计算特征区域的最优匹配，对可能的匹配计算相似变换，将颅骨进行粗对齐；4) 采用改进的ICP算法进行颅骨的细匹配。图 1为本文具体的流程图。

本文使用的数据为CT扫描的3维点云颅骨数据模型，是西北大学与陕西中医学院附属医院合作采集的，每一个志愿者都是颅骨和面皮整套数据，记录志愿者的年龄、性别、祖籍等完整数据。使用Geomagic软件直接重构出来的颅骨模型内部有很多噪声，只能进行手动去除，如图 2 (a) 所示噪声。针对该问题，项目组基于C++编写了颅骨重构程序，该系统自动提取颅骨的外轮廓进行颅骨重构。具体步骤为：1) 对采集到的颅骨CT数据进行去噪、二值化和边界提取处理，使其成为易于实验分析的单层颅骨轮廓数据；2) 使用8邻域搜索提取每一层颅骨轮廓；3) 整合每一层的颅骨轮廓，形成整体颅骨点云数据；4) 针对眼睛、鼻子等噪声较大的区域，结合k-d Tree使用8邻域搜索，最终得到干净的，仅有外表皮的颅骨数据。如图 2 (b) 所示。

去噪后的颅骨数据，一般为30万50万个点, 非常影响配准效率。本文使用一种结合点云曲率和法向的数据精简算法^[17]对颅骨点云进行简化。该算法通过提出一种能够确定最佳聚类数的改进K-means聚类算法构建点云数据点的空间拓扑关系；其次，综合考虑类内点云数据点的曲率及法向量方向，并将两者加权获得特征因子，从而通过特征因子来判断3维数据点是否为特征点；最后，通过设定特征因子的阈值来简化不同类内的点云数据，从而得到含特征的点云数据。如图 3所示，图 3(b)是图 3(a)简化后的结果，简化率为50%，当点数减少一半的时候，形状没有发生明显的变化。该算法很好地保留了数据的细节特征和几何形状。

由于每次数据采集的设备不同，个体的头颅摆放位置也不完全一致，造成采集得到的数据可能在不同的体坐标系，本文采用法兰克福坐标系进行坐标统一。法兰克福平面${P_{\rm{F}}}$由3点${L_p}$、${R_p}$、${M_p}$确定

$ {P_{\rm{F}}} = {M_p}{L_p} \times {M_p}{R_p} $

(1)

式中，${L_p}$、${R_p}$、${M_p}$3点分别表示左、右双耳孔中点和左眼眶下缘中点，另外选取眉心点为${V_p}$点，建立图 4(a)所示的法兰克福坐标系，颅骨模型如图 4(b)所示。

建立步骤如下：

1) 选取${L_p}$、${R_p}$、${M_p}$和${V_p}$4点。由${L_p}$、${R_p}$、${M_p}$ 3点可确定Frankfurt平面记作${P_{{\rm{laneF}}}}$；

2) 确定新坐标系原点以$\overline {{L_p}{R_p}} $为法向量且过${V_p}$的平面与${L_p}{R_p}$的交点为新的坐标原点，记为${O_n}$；

3) 确定新坐标轴${L_p}{R_p}$为$X$轴；过新坐标原点且以${P_{{\rm{laneF}}}}$法向作为其方向向量的射线作$Z$轴；过新坐标原点且同时垂直于$X$轴、$Z$轴的直线为$Y$轴；根据右手规则即可确定$Y$轴的正向。

统一坐标系后，再进行尺度归一化。设定所有颅骨模型的${L_p}-{R_p}$之间的距离为单位1，则对颅骨尺度变换

$ \left[\begin{array}{l} x'\\ y'\\ z' \end{array} \right] = \frac{1}{{{L_p} - {R_p}}}\left[\begin{array}{l} x\\ y\\ z \end{array} \right] $

(2)

式中，$\left( {x, y, z} \right)$是缩放前点的坐标，$(x\prime, y\prime, z\prime ')$是缩放后点的坐标。

由于颅骨表面有凹凸度不同的区域，会包含较多的特征点，本文根据颅骨表面的这一特点来实现颅骨的粗配准。

本文使用体积积分不变量来对每个点进行凹凸性判断。如图 5所示，$\boldsymbol{S}$为3维实体模型$D$的表面曲面，以$\boldsymbol{S}$上任意一点$p$为中心，以$r$为半径的球体为该点$p$的球邻域，则该点$p$的体积积分不变量公式为

$ {V_r}\left( p \right) = {\smallint _{{B_r}(p)}}{I_D}\left( x \right){\rm{d}}x $

(3)

式中，${B_r}(p)$表示以当前点$p$为中心，以尺度$r$为半径的球邻域。

积分不变量的几何意义为：球${B_r}(p)$在曲面外侧部分的体积，如图 5所示。${V_r}(p)$的值与曲面在点$p$附近的凹凸程度有关，与$r$的大小成正比，且${V_r}(p)$不受${B_r}(p)$内部的噪声点的影响。因此，${V_r}(p)$反映了点$p$邻域曲面的凹凸程度，并可进行不同半径尺度的计算。

曲面上顶点的凹凸性可以根据积分不变量的计算来进行判定。这里定义${\phi _r}\left( p \right)$来描述顶点$p$在半径$r$的区域内的凹凸性，定义式为

$ {\varphi _r}\left( p \right) = \frac{3}{{4{\rm{\pi }}{r^3}}}{V_r}\left( p \right)-\frac{1}{2} $

(4)

当${\phi _r}\left( p \right) > 0$时，${V_r}\left( p \right) > \frac{2}{3}{\rm{\pi }}{r^3}$，$p$为凸顶点；当${\phi _r}\left( p \right) < 0$时，${V_r}\left( p \right) < \frac{2}{3}{\rm{\pi }}{r^3}$，$p$为凹顶点；当${\phi _r}\left( p \right) = 0$时，${V_r}\left( p \right) = \frac{2}{3}{\rm{\pi }}{r^3}, {\phi _r}\left( p \right) \in [-\frac{1}{2}, \frac{1}{2}]$，$p$为平面顶点。在实际计算中，需要设置阈值$\varepsilon $，当${\phi _r}\left( p \right) > \varepsilon $时，$p$为凸顶点；当${\phi _r}\left( p \right) <-\varepsilon $时，$p$为凹顶点；当$-\varepsilon \le {\phi _r}\left( p \right) \le \varepsilon $时，$p$为平面顶点。

特征区域指的是一个连通的凹的或凸的或平面的局部区域。根据2.1节计算体积积分不变量，利用K-means分类方法把点分为凸特征点和凹特征点，提取特征区域的算法如下：1) 计算颅骨表面每一个顶点的体积积分不变量${V_r}(p)$和${\phi _r}\left( p \right)$。

2) 根据${\phi _r}\left( p \right)$对颅骨表面的顶点进行分类，将满足${\phi _r}\left( p \right)$的凸顶点记入到集合${\boldsymbol{\phi} _a}$中，将满足${\phi _r}\left( p \right) <-\varepsilon $的凹顶点记入到集合${\boldsymbol{\phi} _b}$中。

3) 按照顶点的连通性，分别对集合${\boldsymbol{\phi} _a}$和${\boldsymbol{\phi} _b}$中的顶点进行K-Means聚类，将颅骨表面划分为多个或凹或凸的特征区域。若某个特征区域顶点数少于给定阈值，可去掉该区域。

设待复原颅骨和相似颅骨的局部区域集合分别为$\boldsymbol{R}_u$和$\boldsymbol{R}_s$，$\boldsymbol{R}_u$和$\boldsymbol{R}_s$分别属于不同的颅骨表面${\phi _u}$和${\phi _s}$的区域，若区域集合$\boldsymbol{R}_u$和$\boldsymbol{R}_s$中的区域同时满足：

1) 类型相同，即都为同凸或者同凹的区域，不考虑平面区域的情况；

2) 具有相近的主成分；

3) 面积大小接近。

则认为它们是相似区域。

本文采用主成分分析法 (PCA) 来比较特征区域的相似性。粗配准步骤如下：

1) 计算特征区域顶点集的主成分信息。设特征区域$\boldsymbol{R}$，包含$n$个顶点，记为${p_i} \in \boldsymbol{R}(i = 1, 2, \ldots, n)$，顶点的坐标记为$({x_i}, {y_i}, {z_i})$，则$\boldsymbol{R}$的质心为

$ \bar p = \frac{1}{n}\sum\limits_{i = 1}^n {{p_i}} $

(5)

$\boldsymbol{R}$的3×3的协方差矩阵$\boldsymbol{M}$为

$ \boldsymbol{M} = \sum\limits_{i = 1}^n {\left[{\left( {{p_i}-\bar p} \right)\cdot{{\left( {{p_i}-\bar p} \right)}^{\rm{T}}}} \right]} $

(6)

计算$\boldsymbol{M}$的特征值为${\lambda ^R}_0 \ge {\lambda ^R}_1 \ge {\lambda ^R}_2$，可定义为区域$\boldsymbol{R}$的主成分，对应的特征向量为$\boldsymbol{n}_0^R,\boldsymbol{n}_1^R,\boldsymbol{n}_2^R$，定义为区域$\boldsymbol{R}$的主方向。

2) 比较特征区域的相似性。定义尺寸特征$\phi \left( \boldsymbol{R} \right)$和各向异性特征$A(\boldsymbol{R})$为

$ \phi \left( \boldsymbol{R} \right) = {\left( {\lambda _0^R + \lambda _1^R + \lambda _2^R} \right)^{1/2}} $

(7)

$ A\left( \boldsymbol{R} \right) = {\left| {\frac{{{\lambda ^R}_1}}{{{\lambda ^R}_0}}} \right|^{1/2}} $

(8)

设有两个局部区域${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$分别属于集合$\boldsymbol{R}_u$和$\boldsymbol{R}_s$，定义$\phi \left( {{\boldsymbol{R}_1}, {\boldsymbol{R}_2}} \right)$和$A\left( {{\boldsymbol{R}_1}, {\boldsymbol{R}_2}} \right)$为

$ \phi \left( {{\mathit{\boldsymbol{R}}_1}, {\mathit{\boldsymbol{R}}_2}} \right) = \left| {\frac{{S\left( {{\mathit{\boldsymbol{R}}_1}} \right)-S\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}{{S\left( {{\mathit{\boldsymbol{R}}_1}} \right) + S\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}} \right| $

(9)

$ A\left( {{\mathit{\boldsymbol{R}}_1}, {\mathit{\boldsymbol{R}}_2}} \right) = \left| {\frac{{A\left( {{\mathit{\boldsymbol{R}}_1}} \right)-A\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}{{A\left( {{\mathit{\boldsymbol{R}}_1}} \right) + A\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}} \right| $

(10)

若$\phi \left( {{\boldsymbol{R}_1}, {\boldsymbol{R}_2}} \right) < {\varepsilon _\phi }$且$A\left( {{\boldsymbol{R}_1}, {\boldsymbol{R}_2}} \right) < {\varepsilon _A}$，则认为${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$具有相近的主成分。

设${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$的面积分别为$N({\boldsymbol{R}_1})$和$N({\boldsymbol{R}_2})$，则定义$N({\boldsymbol{R}_1}, {\boldsymbol{R}_2})$为

$ N\left( {{\mathit{\boldsymbol{R}}_1}, {\mathit{\boldsymbol{R}}_2}} \right) = \left| {\frac{{N\left( {{\mathit{\boldsymbol{R}}_1}} \right)-N\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}{{N\left( {{\mathit{\boldsymbol{R}}_1}} \right) + N\left( {{\mathit{\boldsymbol{R}}_2}} \right)}}} \right| $

(11)

若$N\left( {{\boldsymbol{R}_1}, {\boldsymbol{R}_2}} \right) < {\varepsilon _n}$，则认为${\boldsymbol{R}_1}$和${\boldsymbol{R}_2}$的面积大小相近，${\varepsilon _\phi }, {\varepsilon _A}, {\varepsilon _n}$为给定阈值。

3) 根据相似特征计算的结果，一般情况下，得到的相似区域对中会存在一些伪匹配对，为减少寻找最优匹配对的搜索空间，提高匹配的准确性，采用基于距离主方向约束的聚类算法排除伪匹配区域对。由于颅骨表面凹凸区域数量不多，故此处采用穷举法寻找最优匹配，再对每一个可能的匹配计算其3维变换，将两个颅骨表面粗略对齐。

在颅骨表面的细匹配中，ICP算法是比较有效的，能完成精细配准。但是其匹配过程需要不断迭代，直到满足终止条件为止，因此匹配速度较慢。针对该问题，根据文献[7]的经典ICP算法，这里提出一种改进的ICP算法，即在迭代过程中加入迭代因子$h$，在不影响配准精度的情况下，可以有效地提高迭代的收敛速度。

为解决文献[7]的ICP算法在高精度颅骨匹配过程中的时间消耗问题，首先优化寻找相关点的原始相关点，提出一种新的寻找相关点的搜索方法。该方法是在基本ICP算法的基础上，加入一个动态的迭代因子$h$，以此提高迭代速度。

设${\boldsymbol{P}_L} = \left\{ {{p_{li}}} \right\}_{i = 1}^{Lp},{\boldsymbol{P}_R} = \left\{ {{p_{ri}}} \right\}_{i = 1}^{Rp}$是待配准的两组颅骨点云数据，点${p_{li}}$在点集$\boldsymbol{P}_L$中，ICP算法首先在点集$\boldsymbol{P}_R$寻找中寻找它的最近点$y_i$作为${p_{li}}$的相关点。让$Y = \left\{ {{y_i}} \right\}_{i = 1}^{Rp}$作为点集$\boldsymbol{P}_L$的相关点，即最小化目标函数

$ f\left( q \right) = \frac{1}{{Lp}}\sum\limits_{i = 0}^{Lp} {{{\left\| {{y_i}-R\left( {{\boldsymbol{q}_R}} \right){p_{li}}-{\boldsymbol{q}_T}} \right\|}^2}} $

(12)

式中，$\boldsymbol{q}_R$是旋转矢量，$\boldsymbol{q}_T$表示平移矢量。式 (12) 函数解决目标函数最小化并且这个过程表示为$d$，且

$ \left( {q, d} \right) = Q\left( {{\boldsymbol{P}_L}, {\boldsymbol{P}_R}} \right) $

(13)

式中，$Q$表示最小化操作，字母$d$表示配错的平方。

字母$\boldsymbol{q}$是点集$\boldsymbol{P}_L$和$\boldsymbol{P}_R$之间的刚体转换矢量即$\boldsymbol{q} = {\left[{{q_R}|{q_T}} \right]^{\rm{T}}}$即

$ q = {\left[{{q_0}\;\;{q_x}\;\;{q_y}\;\;{q_z}\;\;{t_x}\;\;{t_y}\;\;{t_z}} \right]^{\rm{T}}} $

(14)

式中，参数向量中四元数参数满足约束条件为

$ {q^2}_0 + {q^2}_x + {q^2}_y + {q^2}_z = 0 $

迭代初始值为${\boldsymbol{P}_0} = {\boldsymbol{P}_L}$，颅骨刚体转换向量为

$ {\boldsymbol{q}_0} = {\left[ {1\;\;0\;\;0\;\;0\;\;0\;\;0\;\;0} \right]^{\rm{T}}} $

在基本ICP算法中加入动态迭代因子$h$的步骤如下：

1) 计算颅骨转换矢量${\boldsymbol{q}_k} = {\left[{{\boldsymbol{q}_R}|{\boldsymbol{q}_T}} \right]^{\rm{T}}}$，以及$\boldsymbol{q}_k$的相邻两次迭代的变化量$\Delta {q_k}$；

2) 用刚体变换变化量$\Delta {q_k}$更新经典ICP算法，经典ICP使用配准参数得到从$p_0$到新位置中的${\boldsymbol{P}_{k + 1}}$转换，即${\boldsymbol{P}_{k + 1}} = {q_k}({\boldsymbol{P}_0})$。改进的ICP则另执行${\boldsymbol{P}_{k + 1}} = \Delta {q_k}({\boldsymbol{P}_{k + 1}})$共$h$次，$h$通常是一个整数。

加入动态迭代因子$h$后的改进ICP算法的步骤如下：

设初值：令${\boldsymbol{P}_0} = {\boldsymbol{R}^*}_0\boldsymbol{P} + {\boldsymbol{T}_0}$刚体变换矢量$\boldsymbol{q} = {\left[{{\boldsymbol{R}_0}|{\boldsymbol{T}_0}} \right]^{\rm{T}}}$，动态调整因子初值$h = 0$。这里${{\boldsymbol{R}_0}}$, ${{\boldsymbol{T}_0}}$是粗配准后的初始刚体变换值。

1) 计算相关点${Y_k} = C({\boldsymbol{p}_{lk}}, {\boldsymbol{p}_{rk}})$，式中，$C$根据点集${\boldsymbol{p}_{lk}}$中点的坐标，在曲面上搜索相应最近点点集${\boldsymbol{p}_{rk}}$。

2) 计算配准参数$(\Delta {q_k}, {d_k}) = Q({P_{lk}}, {Y_k})$，为了得到$\Delta {q_k}$，将原始ICP算法中的$\left( {{q_k}, {d_k}} \right) = Q({p_0}, {Y_k})$转换成$(\Delta {q_k}, {d_k}) = Q({P_{lk}}, {Y_k})$。这表明新算法从转换从${p_{lk}}$到$Y_k$，而原始的ICP算法转换从$p_0$到$Y_k$。

3) 使用配准参数$\Delta {q_k}$到${p_{lk}}$得到一个新位置，应用配准参数$\Delta {q_k}$到${p_{lk}}$，${P_{l(k + 1)}} = \Delta {q_k}({p_{lk}})$。本算法计算$\Delta {q_k}$直接从${p_{lk}}$到$Y_k$，然而原始ICP从$P_0$到$Y_k$计算$q_k$。

4) 计算连续两次颅骨变换后误差的平方差$\Delta {d_k}$，即$\Delta {d_k} = {d_k}-{d_{k + 1}}$，如果$\Delta {d_k} > 0$, 执行操作$h = h + 1$，否则执行$h = 0$。

5) 判断动态调整因子$h$, 如果$h > 0$，则用刚体变换矢量$\Delta {q_k}$更新基本ICP算法中的${q_i}(k + 1)$共$h$次，通常$h$是一个整数。

6) 判断，若满足迭代终止条件${d_k}-{d_{k + 1}} < \tau $则停止，否则$k = k + 1$转到步骤1)。

实验使用西北大学可视化所308套完整的CT扫描数据。对于一个未知颅骨U，在颅骨库中配出一个或者几个相似颅骨。经多次实验，当动态因子$h$=3时，迭代次数和迭代时间都是最好的，$h$=1相当于经典ICP算法，$h$大于3时，配准效果不佳。这是由于随着迭代因子的增大，收敛速率会越来越快。但是当迭代因子增加到一个确定的值，当再增大迭代因子时收敛速率不会再加快并且出现震荡现象，即匹配错误率变大。图 6—图 9显示了待配准颅骨与颅骨库中的一个较相似的颅骨S详细的配准过程。表 1为待配准颅骨U和颅骨S配准的详细信息。

为了验证本文所提算法的普适性，采用了兵马俑3号坑百戏俑的5号俑右腿和右脚匹配为补充，将两个断裂面配准。由于断裂面比较复杂，凹凸曲面较多，不均匀且没有规律，其中有些地方已经磨损或者缺失。经过多次实验，$h$=4最为配准效果最好且配准时间也比较快，具体过程如图 10—图 13所示。表 2为$h$=4时6次配准的详细信息。

对于公共数据集，本文选择兔子做实验，多次实验得到$h$=3时效果最佳。图 14—图 17为配准的详细过程。表 3是算法运行分析表。

针对从颅骨寻找相似颅骨的问题，提出一种基于特征区域和改进ICP算法的颅骨配准算法。针对ICP算法的改进效果，本文使用了3套数据，即颅骨数据、兵马俑数据、公共数据集。从实验1(图 8、图 9) 的配准效果可以看出，改进ICP算法比经典ICP算法配准效果稍微好一些。而且改进ICP算法对兵马俑的数据同样适用实验2(图 12、图 13)，配准效果比原始ICP配准效果好。最后使用的公共数据集兔子，从实验3(图 16、图 17) 的对比可以看出，改进ICP算法配准效果更好。

从3个运行分析表可以看出，虽然改进ICP算法在迭代次数上与经典ICP算法相比提高不大，但改进ICP算法在配准时间比经典ICP算法快很多，这主要是因为加入了迭代因子$h$，从而进行下一次迭代时是在前一次迭代的基础上进行，节约了计算转换矩阵的时间。这样也便于从颅骨库中快速找到相似颅骨。对于兵马俑数据同样适用，可以快速地从兵马俑碎片中找到与其对应的断裂面。

本文提出了一种基于特征区域和改进ICP算法的颅骨表面配准方法。首先，将颅骨模型去噪简化并归一化，通过计算体积积分不变量，将颅骨划分为多个或凹或凸的特征区域，使用K-Means聚类分相似区域。再通过主成分分析法来计算相似的特征区域，鉴于特征区域的个数不多，本文采用穷举法寻找特征区域间的最优匹配，通过对每一个可能的匹配计算3维变换，将两个颅骨面粗略对齐，最后采用加入迭代因子的ICP算法进行颅骨表面的细配准。实验采用西北大学可视化所颅骨库中的颅骨进行配准，并使用兵马俑碎块数据和公共数据集进一步验证本文算法。本文提出的配准算法通过计算点的体积积分不变量并聚类形成特征区域，避免了人工手动标定特征点，在细配准阶段取通过加入迭代因子取得了较高的迭代收敛速度。实验结果表明，对不同的数据集需要设定不同的迭代因子，对迭代因子的自动反馈设定还需要进一步研究。因此该算法对颅骨的配准效果较好，数量较多的兵马俑碎块断裂面的自动匹配的精度上还存在一些问题，因此有待进一步的深入研究。