随着计算机技术的发展, 颅面复原技术有了很大的进步.基于模板变形的方法是目前计算机辅助颅面复原主流方法之一.该方法需要在颅面数据库中选取与待复原颅骨最为相似的参考颅骨模型以及对应的面皮模型, 然后基于颅骨特征点得到参考颅骨向待复原颅骨的非刚性变换, 最后将该变换应用在参考面貌模型上, 实现待复原颅骨的面貌复原.

针对参考颅骨的选择, Liang等^[1]在实现颅面复原时提出了基于特征点的径向基函数(Radial basis function, RBF)插值方法, 但是并未对参考颅骨的选择做过多探讨. Pei等^[2]为了使复原之后的颅面具有更好的个性特征, 采用了人脸局部的模板, 但是重点放在了重构时局部模板的研究上.税午阳等^[3]使用配准的方法进行颅面复原, 但在最相似颅骨的选择时使用的是默认的最相似颅骨. Duan等^[4]使用形状参数回归建模进行三维人脸重建, 分别构建颅骨和颅面的统计学形状模型, 得出它们的对应关系, 待复原颅骨根据这种对应关系进行颅面复原, 并没有真正从数据库中选择出最相似的颅骨.朱新懿等^[5]虽然对参考颅骨的选择方法进行了探讨, 但主要是对颅骨上的特征点进行求距离运算, 在进行参考颅骨选择的时候度量方法过于简单, 且同一结果可能会选择产生多个不同的参考颅骨.

颅骨属于三维模型的一种, 因此三维模型检索的方法也同样可以用在颅骨上, 但是当前的三维模型中有很多是直接用三维特征^[6-9]来检索, 时间耗费很多.不仅如此, 颅骨虽然属于三维模型的一种, 但是当前的三维模型检索主要是针对不同类的模型进行检索, 例如很多文献^[10-12]使用的都是基于三维的数据测试集.而对于颅骨的检索, 主要是在颅骨这一类三维模型中进行检索, 各个颅骨之间总体形状相似但差异细微, 一般三维模型检索算法并未考虑该特征.

通过研究发现, 在基于三维检索的基础上, 进行三维转二维^[13-15]的方法可以明显降低检索时间; 同时根据颅骨本身具有的特点进行适用于颅骨特征的选择, 会显著提高颅骨检索的准确性.

针对通常的三维模型算法检索耗时多且颅骨之间差异的特殊性.本文提出一种使用三维颅骨的二维特征来检索最相似参考颅骨的方法.具体思路为: 1)将待复原颅骨和参考颅骨分别建立法兰克福坐标系, 分别获取颅骨$XOY$, $XOZ$, $YOZ$三个平面的正反面总共6个方向的正交投影, 同时获取三个坐标面的深度投影; 2)根据正交投影获得轮廓特征, 根据深度投影获得剖面特征; 3)使用轮廓信息得到夹角特征和距离特征(全局距离特征和局部距离特征), 使用剖面信息得到空域特征; 4)将各个特征进行融合加权, 比较待复原颅骨与数据库中参考颅骨的加权结果, 得出最相似的参考颅骨.具体算法流程如图 1所示.

图 1  本文方法的流程示意图

Fig. 1  Our workflow

下载: 全尺寸图片 幻灯片

1.   颅骨投影

1.1   颅骨数据坐标归一化

在颅面采集的时候, 不可避免地会有因采集时间和采集方法等造成的差异, 为了使得采集到的三维数据更加规范化, 更好地进行颅面相似度的比较, 引入颅面形态学中常用的法兰克福平面, 依据该平面建立如图 2的法兰克福坐标系.

图 2  颅骨的法兰克福坐标系

Fig. 2  Frankfurt coordinate of skull

下载: 全尺寸图片 幻灯片

假设有两个颅骨的大小不同, 但是它们的形状一样, 仍认为它们是相似性好的颅骨; 也就是说, 在颅骨相似性检索时, 应该优先考虑颅骨的形状而不是大小, 即形状优先.因为在基于模板变形的颅面复原中, 形状相同而大小不同的颅骨, 可以根据两个颅骨的大小比例程度, 对面皮进行相同比例的缩放来获得待复原颅骨的面貌复原.基于以上问题的考虑, 本文对颅骨进行归一化操作, 通过缩放, 规定每个颅骨中法兰克福坐标系的原点到左耳孔中心是一个单位长度.

1.2   三维颅骨的六视图正交投影

为了更有效地表达颅骨特征和进行高效检索, 对采集到的颅骨进行去噪等预处理之后, 对三维颅骨进行正交投影和深度投影.在正交投影时获得颅骨的轮廓特征, 包括颅骨轮廓的夹角属性和距离属性; 在深度投影时获取颅骨的剖面特征, 得到三维颅骨的空域属性; 最后对这些特征属性进行融合.特征提取与融合算法流程如图 3所示.

图 3  2维特征分析过程

Fig. 3  2D characteristic analysis process

下载: 全尺寸图片 幻灯片

基于人眼的视觉特征, 本文在提取轮廓特征时选取6个正交投影视图, 同时将模拟切割模型应用到投影面上来表示剖面属性, 使用距离-灰度映射的思想来获取三维颅骨剖面上的点, 也就是说, 一个点到坐标平面的距离越近, 该点就越有可能是剖面上面的点.因为剖面没有方向信息, 所以本文在提取剖面及剖面特征时使用三视图正交深度投影.

对三维颅骨数据分别沿着$X$, $Y$和$Z$轴的正负6个方向进行正交投影, 也就是把颅骨分别投影在$YOZ$, $XOZ$和$XOY$这三个面的正反面, 设点为三维颅骨上面的点, 则投影出六视图的过程如下:

1) 方向$X+: X_i\geq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(y_i, z_i)$;

2) 方向$X- : X_i\leq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(y_i, z_i)$;

3) 方向$Y+: Y_i\geq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(x_i, z_i)$;

4) 方向$Y- : Y_i\leq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(x_i, z_i)$;

5) 方向$Z+: Z_i\geq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(x_i, y_i)$;

6) 方向$Z- : Z_i\leq 0$投影: $P_i(X_i, Y_i, Z_i)\rightarrow p_i(x_i, y_i)$.

一般我们认为颅骨是左右对称的, 所以对应在法兰克福坐标系上时, 只需要考虑$X+$或者$X-$中的一个投影即可.本文选择投影在$X+$方向, 所以在具体计算时, 只考虑5个方向的投影, 即$X+$, $Y+$, $Y-$, $Z+$, $Z-$.

1.3   三维颅骨的深度投影

在对剖面特征进行表示时, 本文只考虑灰度比较小的投影点.将经过预处理之后的三维点云模型分别在$YOZ$, $XOZ$和$XOY$平面上进行投影, 在投影时, 用距离$d$表示投影点到投影面的距离, 每个投影点的灰度可以根据这个距离进行一一映射, 本文使用[0, 255]区间表示灰度级别.映射遵循的准则是每个点的灰度与它到投影面的距离成正比, 即距离投影面越近的那个点, 灰度值越小, 否则, 灰度值越大.

在遍历颅骨上的点时找到$X$, $Y$, $Z$坐标中正负的最大绝对值, 即$\left | X+\right |_m$和$\left | X-\right |_m$, $\left | Y+\right |_m$和$\left | Y-\right |_m$, $\left | Z+\right |_m$和$\left | Z-\right |_m$, 设

假设点$P$为颅骨上的点, 那么该颅骨的深度投影步骤如下:

1) 面$XOY$.若$P$到投影面的距离是$d_i=\left | z_i \right |$, 那么它的灰度$f_i=(d_i/M3)\times 255$, 即$P_i(X_i, Y_i, Z_i)$ $\rightarrow$ $P_i'(X_i, Y_i, f_i)$;

2) 面$XOZ$.若$P$到投影面的距离是$d_i=\left | y_i \right |$, 那么它的灰度$f_i=(d_i/M2)\times 255$, 即$P_i(X_i, Y_i, Z_i)$ $\rightarrow$ $P_i'(X_i, Z_i, f_i)$;

3) 面$YOZ$.若$P$到投影面的距离是$d_i=\left | x_i \right |$, 那么它的灰度$f_i=(d_i/M1)\times 255$, 即$P_i(X_i, Y_i, Z_i)$ $\rightarrow$ $P_i'(Y_i, Z_i, f_i)$.

2.   颅骨特征提取

2.1   投影距离特征

本文将投影的距离特征分为局部距离特征和全局距离特征.局部距离特征是指将该投影面分为扇形区域之后, 各个扇形区域的特征; 全局距离特征是指整个二维投影中每个扇形区域里凹凸中点之间的距离, 通过该距离来反映该二维投影整体的属性.

2.1.1   局部距离特征

对获得的六视图投影, 建立二维坐标系.以方向$Y$+, 即$XOZ$正面为例, 法兰克福坐标系的原点、$X$轴和$Z$轴分别为二维坐标系的原点和坐标轴.特征提取按照以下步骤进行操作.

步骤1. 将投影以坐标原点为中心, 划分为$N$个弧度相等的扇形;

步骤2. 在每个扇形中, 将弧度相同的点归为一类, 设一个扇形中弧度相同的点有$m$类, 则弧度相同的点可以表示成$C_j$ $(j=1, 2, \cdots, m)$;

步骤3. 遍历弧度相同的点, 计算到原点的距离, 其中距离最大的点即为投影轮廓上的点, 设为$MC_j$ $(j=1, 2, \cdots, m)$;

步骤4. 找到每个扇形区域中的最凸点和最凹点, 即投影轮廓点到原点距离的最近点和最远点.遍历一个扇形中的$MC_j$, 值最大的即是凸点, 值最小的即是凹点.这两个距离表示的是该扇形区域的距离特征;

步骤5. 局部距离特征的提取:扇形区域的两个特征点(凸点和凹点)到原点的距离分别记作$DisM_i$ $(i=1, 2, \cdots, N)$和$DisN_i$ $(i=1, 2, \cdots, N)$, 将所有扇形区域的凸点和凹点到原点的距离记录下来, 就组成了一个局部距离特征向量, 即$\boldsymbol{DisLoc_i}= (DisM_1, DisN_1, \cdots$, $DisM_N, $ $DisN_N)$, 其中, $i$表示三维颅骨的第$i$个正交投影.

2.1.2   全局距离特征

在上一节分割成$N$个扇形的基础上, 在每一个扇形区域中找到最凸点和最凹点所在的弧度, 设最凸点的弧度为$A(\theta _i)$, 最凹点的弧度为$V(\theta _i)$, 对它们求平均值, 即$Aver(\theta _i)= (A(\theta _i)$ $+$ $V(\theta _i))/2$.在每个扇形的轮廓点$MC_j$中找到弧度为$Aver(\theta _i)$的凹凸中点, 记作$P_i(X_i, Y_i)$ $(i$ $=$ $1$, $2$, $\cdots, N)$.

用$D_{ij}$表示凹凸边界的中点$P_i$到$P_j$的欧几里得距离, 分别计算每个凹凸中点到所有扇区中该类点的距离, 最后, 构建一个$N\times N$的距离矩阵, 显而易见, 这个矩阵是对称的, 而且它的对角线元素全部是0.

在矩阵$D$的基础上, 进行以下改进, 得到多尺度距离矩阵(Multiscale distance matrix, MDM)^[16]:

步骤1. 对矩阵$D$的第一列, 将列内元素进行循环的移动, 直至元素0在第一行, 然后将剩下的每一列按照第一列的方法进行循环移动, 最终使得整个矩阵的第一行都变成0.新产生的矩阵记作$D_m$;

步骤2. 在步骤1的基础上, 删除矩阵$D_m$的第一行和倒数的$\left\lfloor{(N-1)}/{2}\right\rfloor$行, 删除之后的矩阵就是多尺度距离矩阵.

图 4是一个当$N=4$的时候构建多尺度距离矩阵的例子.在$N=4$中, 每行都可获取该轮廓几何属性的确定范围.例如, 第一行全部是0, 表示的几何属性是每个凹凸中点到自身的距离; 第二行表示每个点与它最近邻扇区中凹凸中点的距离, 同时这一行能最精确的表示一个形状的几何属性; 第三行表示的是每个点到与它间隔一个扇区的凹凸中点的距离, 从第三行开始, 随着行号的增加, 剩下的行属性表示出的形状就变得模糊起来, 直至行号到达$\left\lfloor{N}/{2}\right\rfloor$时, 表示的形状属性是最不精确的.

图 4  多尺度距离矩阵的构建

Fig. 4  The generation of MDM

下载: 全尺寸图片 幻灯片

对于一个封闭的形状, $D_m$的第二行和最后一行是相同的, 它们描述的都是一个点与它最近邻点的距离.所以在$D_m$中, 大约一半的行是多余的.

图 5是一个从颅骨$Y+$方向投影轮廓中产生的$MDM$的例子, 其中, 按照上述规则在边界上采样64个点, 因此这个矩阵的大小是$32\times 64$.图 5 (b)是颅骨二维轮廓.在图 5 (e)的色谱中, 不同的颜色代表不同的数值, 而且随着颜色的加深, 数值变得越来越大.选择矩阵中的四行来展示颅骨的属性, 分别是第1, 8, 16, 32行, 分别对应图 5 (a)、图 5 (c)、图 5 (d)和图 5 (f).显然, 第一行很好地表示了颅骨的轮廓特征, 最后一行则几乎不能表示颅骨的轮廓特征.

图 5  颅骨二维轮廓的多尺度距离矩阵

Fig. 5  The MDM of skull 2D contour

下载: 全尺寸图片 幻灯片

2.2   夹角特征

在第2.1节获取距离特征的基础上, 将每一个扇形区域中比较大的边界所在的弧度设为: $B_i=(\alpha $, $2\alpha$, $\cdots, N\times\alpha)$, 将扇形区域的最凸点和最凹点同原点进行连线, 将两条连线与所在扇形区域弧度较大边界的夹角记为$DegA_i$和$DegV_i$，其中$i$表示第$i$个扇形区域.将所有最大夹角和最小夹角记录下来, 形成夹角特征向量, 即$\boldsymbol{ Deg}_i=(DegA_1, DegV_1$, $\cdots$, $DegA_N, $ $DegV_N)$, 其中, $i$表示三维颅骨的第$i$个正交投影.

2.3   剖面特征

在进行剖面特征提取时, 如果某个点到切割面的距离越近, 这个点越有可能落在切割面上, 因为本文使用距离-灰度进行映射, 所以在二维投影中, 灰度值越小的点, 越有可能落在剖面上.即该点越能反映剖面的信息.那么在分析剖面情况的时候, 应该注重灰度值小的那些点.

对剖面的特征提取算法步骤如下:

步骤1. 将剖面上的点根据灰度进行从小到大排序, 灰度值越小, 该点越有可能落在剖面上, 即灰度越小, 越能够反映出剖面的特征.选取灰度在前$\lambda$%的点.满足这个条件的灰度点$V_i(X_i, Y_i, f_i)$集合就组成了剖面点集.

步骤2. 找到剖面的最小外包圆, 将该剖面分割成$N$个扇形.

步骤3. 在扇形的灰度点集$V_i(X_i, Y_i, f_i)$中找距离圆心最近和最远的点, 分别记作$n_i$和$m_i$, 然后用扇形区域所包含点的统计信息来生成扇形格子的参数方程:

其中, $i$表示第$i$个扇形区域.

对所有的扇形区域做类似的变换, 得到所有的$N$个参数方程, 即$\{Z_0, Z_1, \cdots, Z_{N-1}\}$.

步骤4. 对各个参数方程做傅里叶变换, 使得每个参数方程由$N$个频率由低到高的三角函数之和组成.即

其中$Z(k)$是参数方程的另一种表现形式，将各个参数方程变换之后, 即:

步骤5. 本文通过一个$N$维旋转不变的特征描述子来描述一个剖面图, 因为傅里叶函数的$L_2$范数具有旋转不变的特征, 我们选取傅里叶函数的振幅值来对旋转不变性进行特征描述.一个剖面图的特征描述可以表示为

2.4   多特征加权融合

最相似颅骨检索指的是用户在输入待复原颅骨模型的时候, 在数据库中准确高效地检索出与之最相似的参考颅骨.本文使用距离特征和夹角特征来反映颅骨的轮廓信息, 使用距离-灰度映射的特征来表示颅骨的剖面信息.针对多种特征, 本文根据各个特征物理意义的不同, 对其进行恰当的加权融合, 使其能够更全面地表现出一个颅骨的特征.

将待复原颅骨和参考颅骨分别记做$X$和$Y$, 分别对$X$颅骨和$Y$颅骨计算其轮廓特征和剖面特征的差异.

2.4.1   对应特征的差异计算

图 6给出了各个特征进行加权融合的图解, 分别对$X$和$Y$两个颅骨的对应视图和剖面进行差异的计算.

图 6  $X$和$Y$两颅骨对应特征的加权差异

Fig. 6  The weighted differences of the corresponding features for $X$ skull and $Y$ skull

下载: 全尺寸图片 幻灯片

局部距离差异:

其中, $DisLoc_{Xi}$和$DisLoc_{Yi}$分别表示$X$和$Y$两个颅骨的第$i$个正交投影的局部距离差异.

全局距离差异:

其中, $MDM_X$和$MDM_Y$分别表示$X$和$Y$两个颅骨的投影轮廓所对应的多尺度距离矩阵.

角度差异:

其中, $Deg_{Xi}$和$Deg_{Yi}$分别表示$X$和$Y$两个颅骨的第$i$个正交投影的角度差异.

根据式(4)和式(5), 计算单幅视图轮廓的距离差异:

其中, $\varepsilon _1$和$\varepsilon_2$是权重因子, 满足$\varepsilon_1+\varepsilon _2=1$.

$X$和$Y$总的轮廓差异:

其中, $\omega _1$和$\omega _2$是权重因子, 满足$\omega _1+\omega_2=1$的条件.

在计算轮廓差异的过程中, 考虑到三维颅骨投影时不同方向上投影点的数量表示的信息量不同, 本算法利用投影点的总个数来避免这种不一致性.

其中, $SumNum$表示5个投影面中投影点的总个数, $Num_{Xi}$表示$X$模型的第$i$个投影图像中点的个数, $Num_{Yi}$表示$Y$模型的第$i$个投影图像中点的个数.

2.4.2   对应剖面特征差异

$X$和$Y$两颅骨的一组对应剖面的差异计算公式为

$X$和$Y$两颅骨的总的剖面差异可表示为

2.4.3   颅骨总体差异计算

根据以上分析, 融合轮廓差异和剖面差异来计算$X$和$Y$两个颅骨的总体差异, 公式为

其中, $\upsilon _1$和$\upsilon _2$是权重因子, 满足$\upsilon _1+\upsilon _2=1$.

2.5   ICP+TPS试验评估

在对参考颅骨和待复原颅骨的相似性度量时, 本文通过对以往配准方法^[17-19]的研究, 提出使用ICP^[20]结合TPS^[21]的方法进行评估.在确定配准点集后, 采用基于ICP算法的刚性配准和基于TPS的非刚性配准两个步骤结合来进行误差评估. ICP解决两个点集间的刚性配准, 通过迭代计算旋转和平移变换使得定义的度量函数最小. TPS能够保证配准后的数据光滑, 是一种基于配准点集对应关系的非刚性配准算法, 并且是一种全局配准算法.设$X=\{x_i\}$和$Y=\{y_i\}$分别是参考颅骨$Ref$ $=\{ref_i, i=1, \cdots, k\}$和待复原颅骨$Ques_i=$ $\{Ques_ i$, $i$ $=1, \cdots, m\}$的配准点集, 则通过TPS变换使得$f(x_i)=y_i$. TPS函数形式为

其中, $\pi (x)=\pi _0+\sum _{j=1}^k\pi _jx_i^{(j)}$是仿射变换; $R(x)=\sum _{i=1}^n\overline{\omega}ø (\left \| x-x_i \right \|)$是基函数的线性组合; $\left \| x-x_i \right \|$为欧氏距离, 引入$WX^{\rm T}=0 $条件作为约束, 求解线性方程

可以计算系数$\pi $和$\omega $.由式(15)推得参考颅骨配准后的结果为

本文为了验证颅骨配准算法的有效性, 定义配准度量误差为

其中, $dis(f(ref_i)-Ques)$是变形后的颅骨顶点到待复原颅骨的欧氏距离.

2.6   加权系数评估方法

本文在进行多特征融合的加权系数计算时, 使用基于颅骨样本的查询结果分类信息熵结合有监督学习来确定加权系数.系数确定步骤如下:

步骤1. 将颅骨数据库中的数据按照$Dist$值分为$K$类, 每类的颅骨记为$L_k$ $(k=1, 2, \cdots, K)$.在每类中进行样本数据的选择, 选择出的样本数据记为$\overline{L}$, $\overline{L}$中存在的$K$类颅骨满足同一类中颅骨之间的$Dist$值差异很小, 不同类之间的$Dist$值差异较大.

步骤2. 样本中第$i$类颅骨总数为$m_i$, 且第$i$类中第$j$个颅骨样本记为$x_{ij}$, 将颅骨样本表示为: $\overline{L}=$ $\{x_{ij}|x_{ij}$ $\in$ $L$, $i=1, 2, \cdots, K$, $j=1, 2, \cdots, m_i\}$, 则样本总数为$m_0=\sum _{i=1}^km_i $.

步骤3. 用$F_p$表示在颅骨中提取的特征, 其中$p$ $=$ $1$, $2, \cdots, M$, $M$表示总共的特征数.

步骤4. 对$\forall x_{ij}\in L$, 在$\overline{L}$中进行与特征的匹配, 满足匹配条件的颅骨总数表示为$S$, 其中匹配到的属于$L_k$类的颅骨结果数量为$n_{ij}^k$, 计算本次的检索信息熵.

其中, $y_{ij}^k=n_{ij}^k/S$, 信息熵对$F_p$特征的检索效率进行了度量.

步骤5. 计算$F_p$的平均信息熵.

平均信息熵的数值越小, 表明该特征对颅骨的描述能力越强, 区分度也随之增高.

步骤6. 计算与特征相对应的加权系数.

3.   实验结果分析

为了对本文算法的有效性进行验证, 颅骨数据采用西北大学可视化研究所的CT采集颅面数据库, 这些数据均经过前期预处理, 可以使得本文的效果更加显著.实验平台为Microsoft Windows 10 Professional, MATLAB R2013b, CPU为Intel Core i7-3770 3.40 GHz, RAM为4.00 GB.

经过前期的实验, 本文将$N$取值为72, 即投影和剖面被分割为72个扇区.获取剖面特征中的参数λ取9.20.单幅视图局部权重因子$\varepsilon _1=0.381$, 单幅视图全局权重因子$\varepsilon _2=0.619$, 单幅视图轮廓权重因子$\omega _1=0.758$, 单幅视图角度权重因子$\omega _2 =0.242$, 颅骨视图权重因子$\upsilon _1 = 0.348, $颅骨剖面权重因子$\upsilon _2$ $=$ $0.652.$

3.1   本文方法实验结果

图 7是输入Skull_112号颅骨, 使用本文算法从数据库中随机选取50套数据进行检索, 将检索出来的与112号颅骨最相似的4套颅骨, 按照相似度从大到小排序, 表 1给出了最相似的颅骨的$Dist$值.本文算法首先排除了最不相似的颅骨, 并且检索出最相似的颅骨数据, 检索的结果中同时包含Skull_112这个输入颅骨本身的数据, 图 8给出了检索出的相似颅骨的配准效果(随着颜色的加深表示配准误差$Dist$值越来越大, 浅灰色表示配准效果较好, 深灰色表示配准效果一般, 黑色表示配准效果较差).首先从直观上看, 获得了十分理想的实验效果.

图 7  模板颅骨Skull_112和检索出来的最相似的四个颅骨

Fig. 7  Template skull_112 and four most similar retrieval skull

下载: 全尺寸图片 幻灯片

表 1  本文方法与Skull_112颅骨最相似的四个颅骨的$Dist$

Table 1  The $Dist$ of four most similar skulls compared with Skull_112 using the method of this paper

参考颅骨编号	$Dist$
Skull_112	0
Skull_12	0.14
Skull_15	0.24
Skull_11	0.36

下载: 导出CSV 

| 显示表格

图 8  本文方法与Skull_112依次最相似的颅骨配准效果

Fig. 8  The successively most similar skull registration effect compared with Skull_112 using our method

下载: 全尺寸图片 幻灯片

3.2   与其他方法比较

文献[5]专门针对颅骨进行了相似性的计算; 文献[8]使用子空间特征向量进行基于三维特征的模型检索; 文献[11]采用国际标准数据测试集PSB, 使用Curvedness feature (CF)方法在多类模型间进行同类模型的检索; 文献[14]使用Principal manifold (PM)方法对模型的三维特征进行转换, 使用二维特征进行检索.

上述4种算法分别具有其代表性, 因此将本文算法分别与这4种算法进行对比.

3.2.1   检索效果比较

与第4.1节一样, 使用各个算法分别检索与Skull_112号最相似的4套颅骨数据, 如图 9所示.

图 9  各种方法检索出的与Skull_112依次最相似的4个颅骨

Fig. 9  Skull_112 and four most similar retrieval skull using various methods

下载: 全尺寸图片 幻灯片

由图 9可知, 每种方法都检索到了与Skull_112号颅骨最相似的Skull_12号颅骨, 图 10将各个方法检索出的颅骨进行配准并计算其$Dist$值, 检验其相似性.

图 10  各种方法与Skull_112颅骨依次最相似的颅骨配准效果

Fig. 10  The successively most similar skull registration effect compared with Skull_112 using various methods

下载: 全尺寸图片 幻灯片

从图 10可以看出, 除了文献[5]之外, 其余方法检索出的颅骨都是按照相似性从高到低排列, 因为文献[5]的方法是将采集好的所有点分为形状特征点和一般特征点, 然后针对不同特征点计算不同含义的距离, 最后通过这些距离来得出颅骨的相似度, 计算方法较简单; CF方法虽然按照相似性从高到低排列, 但是检索出来的排名第4的相似颅骨与本文、文献[8]和PM方法检索到的颅骨不一致, 这源于CF算法主要应用于在不同类的模型中进行同类模型的检索, 而颅骨是属于同类模型的, 所以检索的效果出现了偏差; 本文方法、文献[8]算法和PM算法都很好地检索到了最相似的4套颅骨数据.

3.2.2   统计效果比较

使用本文方法与各个文献中的方法对颅骨数据库中的300套数据进行相似性比较, 对每种方法得到的排名靠前的50套相似颅骨进行误差分析, 得到的$Dist$误差分布直方图如图 11所示.

图 11  $Dist$误差直方图分布

Fig. 11  The error distribution histogram of $Dist$

下载: 全尺寸图片 幻灯片

由图 11可知, 文献[5]的方法的度量值分布规律不明显, 由于其较简单的计算颅骨特征点间距离的判断方式, 导致几乎各个区间都有一定的样本, 分布零散, 且$Dist$在[0, 0.2)区间的数据只占1/10, 统计效果不理想; 文献[8]的方法充分考虑到颅骨的提取, 统计出的[0, 0.2)区间的数据较多, 且[0, 0.4)区间的占比达3/5, 统计效果较好; 文献[11]中的CF算法可以根据不同模型的曲度特征很好地识别出物体的差别, 但是由于颅骨是同一类模型, 仅依靠曲度不容易分辨出不同颅骨之间的差异, 使得CF算法的子空间特征向量统计效果受到影响; 文献[14]的PM方法使用三维颅骨的二维主流特征进行不同颅骨之间的差异判断, 统计结果显示在[0, 0.4)区间的数据超过一半, 但是[0.4, 0.6)区间的数据也有1/5, 综合来说统计效果一般; 本文方法得到的颅骨$Dist$有3/5都在[0, 0.4)之间, 且随着$Dist$的增加, 颅骨数据越来越少, 统计效果最好.

3.2.3   检索效率比较

本文采用在50个颅骨中检索到一个与待复原颅骨最相似的颅骨来进行各方法时间的对比, 用检索耗费的总时间除以50即是检索一个颅骨的平均时间耗费.各个方法的平均检索时间如表 2所示.

表 2  算法效率分析

Table 2  The algorithm efficiency analysis

算法	平均时间耗费(s)
文献[5]	0.31
文献[8]	3.52
CF	4.53
PM	1.31
本文	1.25

下载: 导出CSV 

| 显示表格

文献[5]中针对颅骨的相似性判断只是进行特征点间距离的运算, 且特征点数量少, 大大降低了运算次数, 所以检索时间耗费最少, 但是由于其准确性不足, 故该方法不能作为寻找相似颅骨的通用方法, 可用于进行颅骨相似性的辅助判断.

文献[8]使用子空间特征向量来进行两个模型的对比, 在进行子空间分割的时候引入了递归的方法, 使得效率变低, 且子空间分割的数量不固定, 针对颅骨这种凹凸较多的模型, 分割出的子空间数量较大, 同时在提取三维特征向量时也会耗费大量的时间.

CF算法首先在模型上采点, 计算出每个点的高斯曲率和平均曲率进而求出曲度值, 接着计算每个点到模型质心的距离, 最后构建距离和曲度的分布矩阵.但是其算法规定三维模型随机采点的数量需达到5 000个, 因此不论在计算曲率曲度还是距离都会耗费过多的时间.且最终进行比较时, 为了更好地区分矩阵间的差异, 使用了曼哈顿距离, 在一定程度上也增加了算法运行的时间.

PM算法对三维模型采用二维主流形的方法进行比较, 算法运行时间主要用在二维主流形构造上, 且构造时将经典的PCA算法进行优化并结合网格逼近等思想.所以运行效率显著提高.

本文方法将三维颅骨进行二维投影, 且对颅骨这种单一模型进行针对性的特征提取.所以本文算法相对上述4种方法体现出更高的运行效率.

4.   结束语

本文在从颅面数据库中选择与待复原颅骨最相似的参考颅骨时, 将三维问题转化为二维问题, 提出一种结合夹角距离信息和剖面信息的二维投影轮廓特征提取算法.在保证检索查准率的同时, 显著提高了轮廓特征提取方法的计算效率.本文在进行投影时, 不仅考虑到了轮廓内部点, 而且考虑到了三维模型的空域信息.同时将模拟切割模型的思想应用于投影面, 融合剖面特征与轮廓特征, 提出了基于多尺度二维投影的三维颅骨检索算法.

最后使用ICP和TPS相结合的方法对检索出的颅骨相似程度进行验证, 使得验证的结果更加准确.本文算法整个过程完全避免了人工干预, 自动实现与待复原颅骨最相似的参考颅骨的选择, 提高了度量效率.