0 引言

通过配准建立三维颅骨模型点对点的精确对应是进行颅骨形态统计分析、颅骨面貌复原的前提和基础。但受采集设备不同以及被采集者头部姿态不统一的影响，颅骨模型在位置、姿态和分辨率上存在较大差异。再加上颅骨形状和拓扑结构复杂，个体差异大，又经常存在牙齿、部分头骨缺失等情况，使得三维颅骨模型的自动配准和稠密点对应一直是一个难题。

目前有很多可用于一般三维模型的自动配准方法^[1-3]，包括刚性和非刚性配准方法。在非刚性配准中，应用最广泛的是迭代最近点(Iterative Closest Point, ICP)算法^[4]，但由于颅骨间的差异是非刚性的，ICP算法难以达到准确的配准结果。在非刚性配准中，基于特征点的配准^[5]是常用的方法，但如何自动获得三维颅骨的特征点同样是一个难题^[6]。

近年来，有学者针对三维颅骨的自动配准进行了研究，提出了一些方法^[7-13]。其中大部分方法^[7-9]采用ICP算法进行颅骨的粗配准，在此基础上计算两组对应特征点，之后采用薄板样条函数(Thin Plate Spline, TPS)算法实现精确配准。这些方法的区别在于特征点的选择以及配准次数。如文献[7]通过执行三次ICP和TPS实现两个三维颅骨模型的配准，其中前两次通过ICP粗配准找到脊线的对应关系，然后采用TPS变形实现颅骨曲率变化大的区域的精确配准，最后一次则通过ICP算法在颅骨光滑区域找到对应点，并再次采用TPS变形实现这些区域的精确配准。文献[8-9]通过ICP粗配准后随机采样对应点作为TPS的特征点，然后进行TPS变形，并迭代多次。文献[10]则是在ICP粗配准后求解解剖学上的对应点作为TPS的特征点。因TPS是全局变形算法，导致每一次的TPS变形将影响前一次的配准结果。为此，文献[11]在ICP粗配准后采用具有紧支撑的径向基函数(Compact Supported Radial Basis Function, CSRBF)代替TPS进行精确配准，从而只对配准误差大的局部区域进行调整。文献[12]则采用一致点漂移(Coherent Point Drift, CPD)算法进行精确配准。但这些基于ICP进行粗配准的方法在两个颅骨初始位置和姿态相差较大时或颅骨存在较多缺失时得不到正确的结果。最近，文献[13]提出了一种3D颅骨体模型的配准方法，首先采用主成分分析(Principal Component Analysis, PCA)对颅骨模型进行坐标校正，然后采用ICP进行粗配准，最后通过非刚性的基于体的变形把采用四面体网格模型表示的颅骨模板配准到体素化的颅骨模型上，该方法可对颅骨初始姿态和位置进行调整，但当颅骨缺失较多，尤其不满足对称性时配准效果不佳。

为实现不同姿态和分辨率以及存在缺损的三维颅骨的自动配准，本文提出了一种基于边缘对应的三维颅骨非刚性自动配准方法。与常用的基于特征点的配准方法不同，该方法无需在待配准颅骨和参考颅骨模型上自动标定特征点，而是通过颅骨边缘类型识别建立两个颅骨模型在边缘上的对应关系，继而获得整个颅骨稠密点云的对应关系。

1 边缘提取及分类

本文的三维颅骨模型采用三角网格模型进行表示。在三角网格模型中，一条边通常被两个三角形共同使用，如果一条边只被一个三角形使用，那么这条边称为边界边，多条边界边首尾相连构成一个孔洞边缘，而位于边界上的顶点称为边缘点。

为了给颅骨边缘自动识别提供统计数据，本文从颅面数据库中随机选择了150个三维颅骨模型，然后对其中每一个模型，根据上述方法提取边缘，并交互地标注边缘类型。其中边缘类型包括6类：眼眶、鼻框、颞骨边缘、上颌骨边缘、下颌骨边缘以及其他，如图 1所示。

之后统计除其他边缘外的每一类边缘的边缘长度均值和标准差(见表 1)，以及各边缘之间最短距离的均值和标准差(见表 2)。其中一条边缘的长度为其包含的边界边的长度总和，而两条边缘Bnd_i和Bnd_j之间最短距离dis=min{distance(p_m, q_n)|p_m∈Bnd_i, q_n∈Bnd_j}，distance(p_m, q_n)为边缘点p_m和q_n的欧氏距离。相对于各边缘中心点之间的距离，各边缘间最短距离受颅骨尺寸和形状个体差异的影响小，更适合用于区分颅骨边缘类型。

在边缘长度以及边缘间最近距离统计数据的基础上，可对三维颅骨模型的边缘类型进行自动识别。从表 1可知，眼眶、鼻框和颞骨边缘长度相近，但上颌骨边缘和下颌骨边缘长度与其余边缘长度差异很大，因此可直接根据边缘长度区分出上颌骨边缘和下颌骨边缘。此外，通过边缘长度还可剔除绝大多数其他边缘。之后，再结合边缘间最短距离统计信息对剩余边缘的类型进行识别。具体步骤如下：

1) 如果三维颅骨模型包含上颌骨边缘，则先任意确定其左右方位。

2) 对剩下的每一条待识别边缘，首先利用表 1的统计信息，根据其边缘长度length确定其可能的边缘类型，即如果|length-avglength_type|≤1.96 stdlength_type，则该边缘可为类型type，其中type为眼眶、鼻框、或颞骨边缘。注意，每一条边缘可能有多种类型。并且，本文进一步把type进行细化，使其分为左右。即如果一条边缘可能的类型是眼眶(颞骨边缘)，则该边缘的类型可以是左眼眶(左颞骨边缘)或是右眼眶(右颞骨边缘)。

3) 对每一种类型(左眼眶、右眼眶、鼻框、左颞骨边缘、右颞骨边缘)，通过步骤2)后可获得对应的一个边缘集合，该集合中每条边缘都可以是这种类型。

4) 对每一种类型，从对应的边缘集合中任意取出一条边缘，可获得若干个边缘组合。

5) 采用两步法从所有边缘组合中筛选正确的边缘组合。首先剔除元素有重复的边缘组合；然后在剩下的每个组合中添加上/下颌骨边缘(如果上/下颌骨边缘存在)，之后计算所有边缘间的最短距离，并判断每个距离值dis是否满足表 2中对应的距离条件：dis-avgdis≤1.96 stddis，满足则记1分，否则记0分; 最终得分最高的组合即为最可能的边缘分类结果。

值得注意的是，由于颅骨的对称性，最终得分最高的通常有两个组合，这两个组合左右对称。而根据上述方法无法确定正确的左右方位，需在下一步的姿态校正中进行确定。此外，由于本方法根据边缘间最短距离确定边缘类型，因此，当三维颅骨中仅包含一条除上/下颌骨边缘之外的边缘时，无法确定其类型。

2 姿态校正

由于三维颅骨模型可通过不同设备或软件获得，因此不同三维颅骨模型的初始位置和姿态可能存在较大差异。为了消除这样的差异，在待配准颅骨的各边缘类型识别的基础上，根据待配准颅骨边缘中心点与参考颅骨模型中对应的边缘中心点，对待配准颅骨进行旋转平移变换，使其与参考颅骨模型位置和姿态相近。但该变换的求解需要至少三对对应点，因此本算法要求待配准颅骨至少包含三条非其他类型的边缘。

此外，为确定正确的左右方位，对两个得分最高的分类结果，依次利用上述方法根据对应的边缘中心点对计算刚性变换矩阵，对待配准颅骨进行坐标变换，并求待配准颅骨上的每一个顶点到参考颅骨上的最近点，以两者之间的欧氏距离为该点的误差，累计所有顶点的误差。比较两次的误差总和，其中误差总和小的对应着正确的左右方位。

图 2显示了姿态校正前后待配准颅骨(浅灰色显示)和参考颅骨(深灰色显示)之间的位置关系。其中：图 2(a)为两个颅骨的初始姿态；图 2(b)为左右方位错误的姿态调整结果；图 2(c)为左右方位正确时的最终结果。可以看出即使两个颅骨模型的初始姿态相差很大，通过调整后，可得到相近的姿态和位置, 从而实现了两个颅骨模型的预配准。

3 精确配准

本文通过两次CPD变换实现两个颅骨的精确配准，进而获得两者的对应关系。CPD算法由Myronenko等^[14]提出，它将两个点集的匹配看作概率密度估计问题，参考点集作为混合高斯模型的中心，目标点集为观测数据点，通过最大化观测数据点的似然概率估计点集之间的几何变换。这种方法在含有噪声和溢出点以及两个点集的顶点数量相差较大时有较强的鲁棒性，因此可用于不同分辨率颅骨和存在缺损的情况。

假设两个点集，其中参考点集Y=(y₁, y₂, …, y_M)^T(表示M×D的矩阵)，目标点集X=(x₁, x₂, …, x_N)^T(表示N×D的矩阵)，M与N可以不相等，D是点的维数，T(Y, θ)是几何变换后的参考点集，T为几何变换算子，θ为几何变换参数。将T(Y, θ)看作混合高斯模型的中心，x为该混合高斯概率模型产生的数据点，其似然概率为：

$ p(x) = \sum\limits_{m = 1}^{M + 1} {p(m)p(x|m)} $

其中p(x|M+1)=1/N，用于考虑缺失点，$p(x|m) = \frac{1}{{{{(2\pi {\sigma ^2})}^{D/2}}}}\exp ( - \frac{{{{\left\| {x - T({y_m},\theta )} \right\|}^2}}}{{2{\sigma ^2}}}),m = 1,2,...,M$, p(M+1)=w(0≤w≤1)，表示缺失部分的比例，p(m)=(1-w)/M，即所有高斯混合模型分量具有相同的先验概率。因此，观察数据点集的联合似然概率为：

$ p(X) = \prod\limits_{n = 1}^N {p({x_n})} $

(1)

对式(1)取负对数，得到负对数似然函数，它是几何变换参数θ和方差σ²的函数：

$ E(\theta ,{\sigma ^2}) = - \sum\limits_{n = 1}^N {\ln } \sum\limits_{m = 1}^{M + 1} {p(m)p({x_n}|m)} $

(2)

通过最大化式(1)，即最小化式(2)可以估计θ和σ²。文献[14]利用最大期望(Expectation-Maximization, EM)算法求解该最优化问题。首先在E-步骤中利用现有的参数估计值θ_old和σ_old²，根据贝叶斯原理计算后验概率p_old(m|x_n)(n=1, 2, …, N)；然后在M-步骤中通过最小化目标函数Q(θ, σ²)更新参数值，重复迭代E-步骤和M-步骤直到收敛或达到用户设定迭代次数。其中M-步骤中的目标函数Q(θ, σ²)是式(2)的上限，具体表达式为：

$ \begin{array}{l} Q(\theta ,{\sigma ^2}) = {\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^M {\frac{{{p_{{\rm{old}}}}(m|{x_n})}}{{2{\sigma ^2}}}\left\| {{x_n} - T({y_m},\theta )} \right\|} } ^2} + \\ \;\;\;D/2 \cdot \ln {\sigma ^2}\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^M {{p_{{\rm{old}}}}} } (m|{x_n}) \end{array} $

其中：

$ {p_{{\rm{old}}}}(m|{x_n}) = \frac{{\exp (\frac{{{{\left\| {{x_n} - T({y_m},{\theta _{{\rm{old}}}})} \right\|}^2}}}{{ - 2{\sigma _{{\rm{old}}}}^2}})}}{{\sum\nolimits_{k = 1}^M {\exp (\frac{{{{\left\| {{x_n} - T({y_k},{\theta _{{\rm{old}}}})} \right\|}^2}}}{{ - 2{\sigma _{{\rm{old}}}}^2}}) + c} }} $

其中$c = {(2\pi {\sigma _{{\rm{old}}}}^2)^{\frac{D}{2}}}\frac{w}{{1 - w}}\frac{M}{N}$

几何变换T可以是刚性变换、仿射变换、或非刚性变换。本文中采用非刚性变换，具体公式和算法流程见文献[14]。

在第一次CPD配准时，两个点集由各自的边缘点构成，但并不是所有的边缘点，而是两个颅骨共有的边缘上的顶点。此外，为避免距离各边缘较远的颅顶在CPD配准中产生错误的变形，本文在参考颅骨顶部随机采样了十余个顶点，然后对其中每个顶点p_s，求其在待配准颅骨上的最近点q_t，但由于待配准颅骨可能存在缺损，因此对每个q_t，反过来求其在参考颅骨上的最近点p_s^*，之后把所有顶点q_t和p_s^*分别加入到目标点集和参考点集中。基于这两个点集采用CPD算法进行非刚性变换，然后把该变换作用到参考颅骨上。经过变换后的参考颅骨在共有边缘处与待配准颅骨能够很好地匹配上，但在距离边缘较远的区域，则未能完全匹配上。为此，本文在整个颅骨区域构建新的目标点集和参考点集，并再次执行CPD配准。为获得新的点集，同时避免两个颅骨因分辨率不同带来的影响，对待配准颅骨上的每一个顶点求其在参考颅骨三角网格上的最近点作为对应点，并以两点之间的距离作为配准误差。为了让新点集中的点尽可能落在配准误差大的区域，对所有对应点对按其配准误差从大到小进行排序，然后从中采样出目标点集和参考点集。此外，为避免点集中的点过于集中，采样时设定每个点集中的顶点之间的距离不小于阈值D。D的大小决定了点集中顶点的个数，D越大则点集中的顶点越少，CPD的配准时间越短，配准准确度也会相应地降低。

图 3显示了图 2中经姿态校正后的两个颅骨的第一次CPD的配准结果，其中图 3(a)显示的是用于CPD配准的两个点集，其中黑色显示的为参考点集，灰色显示的是目标点集；图 3(b)显示的是经CPD配准后的两个点集；图 3(c)显示了配准后的两个颅骨。图 4显示了对应的第二次CPD的配准结果，其中图 4(a)和(b)显示了配准前后的两个点集，图 4(c)是两个颅骨的最终的配准结果。可看出两个颅骨的所有区域基本都匹配上了。

4 实验结果及分析

本文采用C++、Matlab和OpenGL实现了本文方法，所有实验在一台CPU为Intel Xeon E5-2620 v3, 2.4 GHz，内存为32 GB的Dell precision T7810工作站上完成。参考颅骨及三个测试颅骨模型的数据情况见表 3。其中颅骨1和参考颅骨的分辨率相当；颅骨2的分辨率约为参考颅骨的2倍；颅骨3的分辨率约为参考颅骨的4倍。此外，颅骨1形态完整并且与参考颅骨初始位置和姿态相近；颅骨2形态完整但其初始位置和姿态与参考颅骨相差较大；颅骨3存在很大缺失，同时其初始位置和姿态与参考颅骨相差较大。

在所有实验中，第一次CPD的参数为：高斯核宽度β=4.0，规格化权重λ=3.0，噪声权重w=0.1，最大迭代次数iter=100，容忍度tol=1E-10；第二次CPD的参数大部分和第一次相同，不同的是：β=2.0，w=0.05；阈值D为5 mm。

为了与本文方法进行比较，我们实现了目前三维颅骨自动配准中最常用的ICP和TPS相结合的方法(简称ICP+TPS方法^[8])，以及在ICP+TPS前，先采用PCA进行姿态调整的方法(简称PCA+ICP+TPS方法)。其中ICP的迭代次数为100次，随机采样顶点数为10 000；TPS的特征点对数为1 000。采用这三种方法对三个测试颅骨一一进行配准，配准结果分别如图 5~7所示。图中还分别显示了各方法的中间结果，例如对于PCA +ICP+TPS方法，分别显示了经过PCA之后的结果、经过ICP之后的结果以及最终的结果。

从图 5中可看出，当两个颅骨的初始位置和姿态相近时，三种方法都能取得很好的结果，并且ICP+TPS方法与PCA+ICP+TPS方法的结果相差不大，而本文方法在各边缘附近区域的匹配效果要优于另两种方法。而当两个颅骨的初始位置和姿态相差较远时，如图 6所示，则ICP+TPS方法会得到错误的配准结果，但PCA+ICP+TPS方法和本文方法都能得到很好的结果。通过对比这两种方法的最终结果，再次证实本文方法对边缘附近区域的配准更准确。而当两个颅骨的初始位置和姿态相差较远，同时待配准颅骨存在较大缺失时，如图 7所示，ICP+TPS方法和PCA+ICP+TPS方法都会失效，但本文方法依然能得到不错的结果。

本文还统计了这三种方法对三个测试颅骨的配准时间，包括总时间以及分步骤时间(括号内数字)，如表 4所示。其中对于ICP+TPS方法，分步骤时间分别对应ICP配准时间和TPS配准时间；对于PCA+ICP+TPS方法，分步骤时间分别对应PCA坐标校正时间、ICP配准时间和TPS配准时间；对于本文方法，分步骤时间则分别对应边缘提取和分类时间、坐标校正时间、第一次CPD配准时间和第二次CPD配准时间。从表 4可知，当待配准颅骨完整时，三种方法配准时间相差不大；当待配准颅骨不完整，同时与参考颅骨姿态相差较大时，本文方法明显提高了配准效率。

5 结语

三维颅骨模型配准是三维颅骨形态分析以及三维颅骨面貌复原的前提。本文提出了一种三维颅骨网格模型的自动配准方法，在颅骨边缘类型识别的基础上实现颅骨姿态校正，并进一步结合CPD算法逐步实现了颅骨边缘的配准和整个颅骨的配准。本文方法不仅可以用于初始位置和姿态以及分辨率相差较大的颅骨的配准，还可用于有较大缺失的颅骨的配准，从而弥补了现有方法的不足。

但本文方法在姿态校正时要求待配准颅骨至少包含三条边缘。当颅骨边缘数量不满足条件时，本文方法会失效。因此，如何对任意形状的待配准颅骨进行姿态校正，或直接找到它在参考颅骨中的对应区域是下一步工作需要解决的重点问题。此外，目前我们依赖边缘长度统计值和边缘间最短距离统计值对边缘类型进行识别，该方法不能用于仅有一条边缘(除上/下颌骨边缘外)的识别。在后续工作中，还需研究其他的三维模型形状分析方法，从而采用更先进的特征表示和分类方法。