[发明专利]一种基于现代人软组织分布的古人类头骨颅面交互复原方法

权利要求书

1.一种基于现代人软组织分布的古人类头骨颅面交互复原方法，其特征在于，包括以下步骤：

步骤1，古人类头骨的几何修复：

步骤1.1，古人类头骨外表面点云的提取和建模：

步骤1.1.1，在古人类头骨模型上手工标定左耳孔、右耳孔、左眼眶下缘点和眉心点，从而建立法兰克福坐标系，并将古人类头骨模型调整至法兰克福坐标系下，实现模型姿态校正；

步骤1.1.2，计算古人类头骨三维模型的包围盒box＝{x_min，y_min，z_min，x_max，y_max，z_max}，然后以法兰克福坐标系中Z轴作为圆柱体长轴进行等间距圆柱采样，则采样点p_i＝{x_i，y_i，z_i}可表示为以下公式(1)：

其中：r表示颅骨模型包围盒的对角线的长度，θ表示在坐标系XOY平面中的采样角度，dt表示沿Z轴的采样层间距离，layer表示沿Z轴的层数；

步骤1.1.3，对于Z轴上任意点o_i＝{0，0，oz_i}，oz_i为Z轴坐标，定义以p_i点为起点且为方向的射线计算所述射线与古人类头骨模型的交点，并将所述交点作为外表面顶点，当所述射线与颅骨模型产生多个交点时，则选取t最小的交点作为外表面顶点；

步骤1.1.4，针对步骤1.1.3计算获得的古人类头骨的外表面顶点，采用基于最短对角线的三角剖分算法实现相邻层轮廓线上对应顶点的连接，完成颅骨外表面的三维建模，建模过程中依次计算古人类头骨外表面相邻两层轮廓线上邻接四个顶点对角线的长度，选择最短对角线作为顶点连接的方式，从而实现古人类头骨外表面点云的三角剖分；

步骤1.1.5，通过在相邻层中较短的轮廓线上插入顶点的方法实现模型修复，判断各个非锐角三角形中最长边所在轮廓线的位置，采用B样条函数在该轮廓线的最长边中插入新顶点，连接该三角形中最长边对应的顶点和新插入的顶点，实现三角剖分；

步骤1.2，古人类头骨对称平面的计算：

步骤1.2.1，在古人类头骨模型的正中矢状面上手工标定不少于3个特征点，记为landmark_i＝{lx_i，ly_i，lz_i}，定义对称平面方程为z＝a₀x+a₁y+a₂，则采用最小二乘方法拟合对称平面的平面方程，可表示为以下公式(2)：

上式应满足：

即：

求解上述线性方程组，可得a₀，a₁，a₂；步骤1.2.2，利用拟合对称平面将古人类头骨模型分为原始模型和待修复区域两部分，并依据镜像原理计算原始模型的对称数据，所述数据将作为残缺区域几何修复的依据；

步骤1.2.3，分别提取对称数据和待修复区域的边界顶点集合，分别记为V＝{v_c}和U＝{u_d}，定义边界对应点能量项，RX＝{RX_c}为对称数据中全部顶点对应的仿射变换，(v_c，u_d)表示迭代过程中搜寻最近点确定的顶点间的对应关系，为权值；定义局部保刚性能量项为其中F表示由待修复模型的顶点-边组成的邻接矩阵，为克罗内克乘积算子，G＝diag(1，1，1，1)表示对角矩阵；则通过建立基于边界约束的保刚性能量函数E(RX)实现对称数据和待修复模型的非刚性配准；所述保刚性能量函数由边界对应点能量项E_boundary(RX)和局部保刚性能量项E_s(RX)两部分构成，即E(RX)＝E_boundary(RX)+αE_s(RX)，其中α为权值；

步骤1.3，古人类头骨残缺区域的虚拟修复：

步骤1.3.1，基于步骤1.2.3提取的古人类头骨模型的边界点集，确定待修复区域的边界点集；

步骤1.3.2，融合配准后的对称数据和待修复区域，实现残缺古人类头骨几何模型的虚拟修复；

步骤2，古人类头骨与现代人的颅骨的配准：

步骤2.1，现代人颅面数据库中参考颅骨模型的选择：

步骤2.1.1，计算古人类头骨模型和现代人颅面数据库中每个颅骨各顶点的高斯曲率，基于曲率特征和颅骨的几何坐标分割眉弓区域；

步骤2.1.2，计算古人类头骨模型和现代人颅面数据库中每个颅骨眉弓区域各顶点高斯曲率的均值，并将与古人类头骨高斯曲率均值最接近的现代人样本作为参考模型，参考模型包括颅骨模型和面貌模型；

步骤2.1.3，计算步骤2.1.2中参考模型的面部软组织分布，面部软组织分布将作为古人类头骨初步颅面复原结果的依据；

步骤2.2，古人类头骨与现生人参考颅骨的三维模型配准：

步骤2.2.1，在古人类头骨模型和现代人参考颅骨模型上手工标定特征点集，包括头面部区域的眉弓、眼眶、颧骨、鼻骨、脸颊、牙齿和正中矢状面；进一步针对已标定的特征点集，采用曲面参数拟合方法在古人头骨和现代人参考颅骨模型上自动生成对应点集合；

步骤2.2.2，定义特征点能量项E_l(X)＝∑||X_blp_b-lq_b||²，其中lp_b和lq_b分别表示步骤2.2.1中参考颅骨模型和古人类头骨模型的特征点，(lp_b，lq_b)表示参考颅骨模型和古人类头骨模型中具有相同序号的特征点对；定义自动生成的对应点项E_d(X)＝∑||X_gvp_g-ug_g||²；定义配准迭代过程中模型间最近距离对应点项E_c(X)＝∑δ_l||X_lvz_l-uz_h||²，δ₁为权值；定义局部保刚性能量项为其中F表示由待修复模型的顶点-边组成的邻接矩阵，为克罗内克乘积算子，G＝diag(1，1，1，1)表示对角矩阵；则可以建立特征点能量项E_l(X)和自动生成的对应点项E_d(X)、以及配准迭代过程中模型间最近距离对应点项E_c(X)和局部保刚性能量项E_s(X)四部分构成的保刚性能量函数其中γ，ρ和为权值，实现古人类头骨和现代人参考颅骨的非刚性配准；

步骤3，基于参考模型的软组织分布的颅面复原：

步骤3.1，将用于颅面复原的古人类头骨和现代人参考颅骨的对应点集合分别记为QSkullMarks和RSkullMarks，对应点集包括以下两部分：

(1)步骤2.2.1中手工标定的特征点和自动生成的对应点；

(2)采样点，将对步骤2.2.2中建立的古人类头骨和现代人参考颅骨的对应点集合进行等距离采样，为了提高颅面复原结果的准确性，采样过程中将使得点云均匀覆盖整个古人类头骨模型，以使相同序号的采样点具有近似相同的位置；

步骤3.2，依据已经建立的对应点集QSkullMarks和RSkullMarks，首先采用最近点迭代算法实现现代人参考颅骨向古人类头骨的刚性变换，计算其刚性变换X_R；实验中定义刚性变换能量函数为其中j表示步骤3.1中古人类头骨和现代人参考颅骨对应点集的序号，n表示对应点的数量；采用奇异值分解方法计算旋转变换R和平移变换T；然后，采用薄板样条函数实现经过刚性变换后的现代人参考颅骨向古人类头骨的非刚性变形，计算其仿射变换X_NR；实验中定义薄板样条函数为通过求解线性方程组计算薄板样条函数的系数；最后，将变换矩阵X_R和X_NR应用于现代人的参考面貌模型，实现古人类头骨面貌的初步三维复原，复原结果记为CoarseFace；

步骤3.3，将古人类头骨和该颅面复原结果变换在法兰克福坐标系下，并计算古人类头骨各顶点的软组织厚度；

步骤4，颅面复原结果候选数据集的构建：

步骤4.1，选择步骤2.1.2获得的现代人参考面貌模型为参考模型；

步骤4.2，利用非刚性配准算法实现现代人颅面数据库中三维面貌模型点云间的稠密对应，使得每个三维面貌模型由相同数量的顶点构成且相同序号的顶点具有近似的解剖学位置：

步骤4.2.1，依据颅面数据库中已经标定的面貌特征点集，建立特征点约束的保刚性能量函数E(FX)，实现步骤4.1中的参考面貌模型向颅面数据库中其他三维面貌模型的非刚性配准和点云稠密对应，所述保刚性能量函数由三部分组成，记为E(FX)＝E_l(FX)+χE_c(FX)+φE_s(FX)，其中χ和φ表示权值；FX为参考面貌模型顶点的仿射变换，保刚性能量函数项包括：特征点能量项E_l(FX)、模型间最近距离对应点项E_c(FX)和局部保刚性能量项E_s(FX)；

步骤4.2.2，利用步骤4.2.1的非刚性配准结果，计算现代人颅面数据库的平均面貌模型，并将其作为新的参考模型；

步骤4.2.3，重复步骤4.2.1的过程，实现步骤4.2.2中的新的参考模型向现代人颅面数据库中每个三维面貌模型的非刚性配准和点云稠密对应，点云稠密对应后的三维面貌模型顶点数量与步骤2.1.2获得的参考面貌模型的顶点数量相同，相同序号的顶点具有近似的位置；

步骤4.3，将颅面数据库中每个样本的软组织分布作为依据，实现古人类头骨的颅面复原：

步骤4.3.1，依据现代人颅面数据库中已经标定的颅骨特征点集，建立特征点约束的保刚性能量函数，实现步骤2.1中的参考颅骨模型向现代人颅面数据库中其他颅骨模型的非刚性配准和点云稠密对应；

步骤4.3.2，依据步骤2.2.2获得的古人类头骨和现代人参考颅骨顶点间的对应关系，建立现代人颅面数据库中其他颅骨模型与古人类头骨顶点间的对应关系；

步骤4.3.3，依据步骤3.1中基于颅骨配准的颅面复原方法，将颅面数据库中每个样本的软组织分布作为依据实现古人类头骨的面貌复原，建立颅面复原结果候选数据集，记为EstimateFace＝{eface_k}，k＝0，1，2，...，t，其中eface_k表示利用颅面数据库中第t个样本的软组织分布实现颅面复原的结果；

步骤5，基于统计形状模型的颅面交互复原：

步骤5.1，针对步骤4.3.3建立的颅面复原结果候选数据集EstimateFace，采用主成分分析方法计算主成分及eface_k对应的主成分系数，建立颅面复原结果的统计形状模型，则步骤3.2获得的初步三维复原模型可表示为其中表示颅面复原结果候选数据集EstimateFace的平均面貌模型，wθ＝[wθ₁，wθ₂，...，wθ_s]∈R^3×m×s表示EstimateFace构成的协方差矩阵对应的特征向量，Coarsefaceb_s表示初步颅面复原面貌CoarseFace对应的主成分系数，s表示主成分个数，其值将由贡献率决定；

步骤5.2，针对颅面复原结果候选数据集EstimateFace的主成分，通过设置各主成分系数的权值，发现主成分与三维面貌几何形状的关系；

步骤5.3，通过交互地编辑Coarsefaceb_s的数值，实现古人类头骨面貌复原模型的交互编辑，新的颅面复原结果表示为：其中μ_s表示对每个主成分系数进行编辑的值，l表示三维面貌模型归一化时的系数；

步骤5.4，将古人类头骨三维模型和最终复原三维面貌模型叠加并透明显示，实现颅面复原结果的评价。