CN106780591B - 一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法 - Google Patents

Abstract

本发明涉及一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，包括如下步骤：1颅骨稠密对应点云；2面貌稠密对应点云；3颅面形态关系可视分析；4基于软组织分区的颅面形态关系表示；5未知身源颅骨的面貌复原。本发明利用基于主成分系数的可视分析和最小二乘回归的定量表示方法，实现了颅面形态关系的分析，解决了颅面点云数据量大、主成分几何意义不易确定、颅面形态关系不易定量表示的问题。利用基于软组织厚度的颅面分区方法，克服了颅面形态关系复杂、不同区域间颅面形态关系不一致的问题，提高了颅面形态关系的准确性。最终利用基于分区的颅面形态关系，实现了未知身源颅骨的面貌复原。

Description

一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法

技术领域

本发明涉及一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，属于计算机图形学、数字几何处理、人体测量学技术领域。

背景技术

人体头面部几何形态复杂，颅骨形态决定了面貌的基本形态，颅骨和面貌形态间存在相关性，对其关系进行研究已经成为法医人类学领域研究的热点。传统研究方法大多针对尸源进行软组织厚度测量、头面部测量指标统计分析，然而该方法存在研究对象样本有限，测量指标数量少的局限。随着医学影像设备的发展，CT、MRI等影像设备已经被广泛应用于活体样本的数据采集和三维建模，为新技术和新方法的开展提供了数据基础。

目前，针对颅面形态关系的研究方法主要包括三大类：1，软组织厚度测量和统计。法医人类学家利用计算机软件在影像数据或者三维模型上标定特征点，然后统计不同年龄、性别、种族人群特征点处的软组织厚度。目前，各国学者已经测量了澳大利亚、葡萄牙、埃及、高加索、捷克、土耳其、匈牙利、斯洛伐克、印度、美国、法国、日本、巴西、韩国、巴基斯坦、中国多个国家或民族的软组织厚度，并通过t-test、多元回归分析方法研究面部软组织的分布，借此发现颅骨和面貌之间的相关性。为了保证软组织测量的一致性，Stephan基于已经发布的软组织厚度数据集，统计分析了6700余个样本的软组织厚度。2，头面部几何测量指标的统计。法医人类学家利用计算机软件在影像数据或者三维模型上标定颅骨和对应的面貌特征点，然后分别测量颅骨几何测量项和对应的面貌几何测量项，包括距离、角度、面积的指标，进而分析颅骨和面貌形态之间的关系，如颅骨形状和脸型的关系、颅骨五官和面貌五官之间的关系。3，基于稠密点云的颅面形态表示。Berar等在2006年针对颅面稠密对应点云首次提出基于线性计算模型的颅面形态关系方法，该方法将颅骨点云和对应的面貌点云看作一个向量，通过最小二乘拟合方法实现颅面形态关系的表示。考虑到颅骨五官与面貌五官间更为复杂的关系，Zhang在Berar方法的基础上将颅骨和面貌模型按照五官分区，将每部分的颅骨点云和对应的面貌点云看作一个向量。为了显示地表示颅面形态关系，针对颅骨和面貌模型，利用偏最小二乘回归方法建立颅骨点云坐标和面貌点云坐标间的相互关系，然而由于颅骨和面貌点云数量大，给表示模型计算造成了极大困难。因此，依据统计形变模型(statistic morphable model)的思想，现有方法大多将颅骨和面貌模型分别表示为稠密对应点云，进而利用主成分分析方法实现颅骨和面貌模型的降维表示，最后通过机器学习方法学习颅骨样本主成分系数和面貌样本主成分系数间的关系，实现颅面形态关系的定量表示。该类方法的核心是利用主成分系数表示颅骨和面貌的几何形态。目前，支持向量机、特征根回归、偏最小二乘回归等机器学习方法已经被用于颅面形态关系的定量表示。段福庆将颅骨和面貌模型按照五官分区，然后利用偏最小二乘回归方法实现颅面形态关系的学习。邓擎琼建立了考虑年龄、性别、身体质量指数的回归模型，提高了表示的准确性。

但上述方法仍不能有效解释颅骨主成分和面貌主成分的几何意义，不能有效解释基于主成分和线性回归的颅面形态分析和表示方法的可行性。另一方面，由于颅面形态关系复杂且面部软组织在各区域的分布并不一致，在头部整体数据集上采用主成分分析技术进行表示将降低颅面形态关系表示的准确性。因此需要依靠面部软组织分布，建立一种更有效的颅面分区方法，使得各分区内形态关系尽可能一致，从而提高颅面形态关系表示的准确性。此外，颅骨稠密对应点云和面貌稠密对应点云的建立方法还亟待加强研究。

发明内容

本发明的目的是针对现有技术存在的问题，提供一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，本发明通过特征自动提取和基于隐式函数和顶点能量约束最优结合的非刚性配准算法，实现颅骨点云间和面貌点云间的稠密对应。通过可视和定量分析结合的方法，解决基于主成分分析的颅面形态关系几何意义难以解释的问题。通过基于软组织分区的颅面形态关系线性学习方法，解决复杂颅面形态关系难以准确表示的问题，提高无身源颅骨面貌复原的准确性。本发明应用于刑事案件、考古学中未知身源颅骨的面貌复原，并能推广到应用于医学领域中人体器官、肿瘤等形状的几何分析，本发明包括以下步骤：

步骤1：颅骨稠密对应点云；

步骤2：面貌稠密对应点云；

步骤3：颅面形态关系可视分析；

步骤4：基于软组织分区的颅面形态关系表示；

步骤5：未知身源颅骨的面貌复原。

进一步地，所述步骤1包括以下具体步骤：

步骤1.1：针对颅骨三维网格模型，定义内点和边界点，实现颅骨孔洞边缘的自动识别；并通过孔洞边界长度、孔洞形状和位置，分析确定孔洞之间的对应关系；

步骤1.2：基于高斯映射和动态区域增长算法，实现颅骨几何形状变化较大区域包括上颌、下颌、颧骨等的自动分割；

步骤1.3：从颅骨模型数据集中选择选择两个模型，一个作为参考颅骨，另一个作为目标颅骨，依据步骤1.1和步骤1.2产生的两个颅骨模型特征，利用最近点迭代算法实现参考颅骨向目标颅骨的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个颅骨模型的非刚性配准；

步骤1.4：计算目标模型的每个顶点与变形后的参考模型的最近点，记录其顶点序号，建立目标模型与参考模型顶点的对应关系；

步骤1.5：依据步骤1.4计算获得的顶点序号，针对原始参考颅骨模型，生成与目标颅骨点云对应的参考颅骨顶点坐标；

步骤1.6：从颅骨模型数据集中选择其他模型作为参考颅骨，重复步骤1.3-步骤1.5，直到遍历所有模型为止，从而建立颅骨模型间的顶点对应关系，即所有颅骨模型具有相同的顶点个数且对应顶点具有近似的解剖学位置。

所述步骤1.1进一步包括：

步骤1.1.1：定义颅骨三维网格模型为skull＝{P,E}，其中P＝{p₁,p₂,...,p_n},p_i＝(x_i,y_i,z_i)∈R³表示n个颅骨顶点，E＝{e_k＝(p_i,p_j),k＝1,2,3,...,m}表示m条边。孔洞识别前首选对模型进行检验，判断其是否满足流行，并记录不满足条件的边序号。孔洞识别过程中，首先遍历模型中各个顶点p_i，确定其对应的邻接顶点集合Adj_pi，如果Adj_pi中的点通过模型的边e_k直接连成封闭的多边形，则p_i为内点，否则为边界点。然后，以任意边界点p_j为起点，根据e_k寻找下一个边界点，直到遍历完所有边界点。最终则会获得c条封闭边界，记为boundary_i＝{p_j},i＝1,2,3,...,c。

步骤1.1.2：计算每条边界的长度、中心，通过分析边界长度和中心位置即能识别左眼眶、右眼眶和鼻骨三个轮廓。通过统计训练集中样本的左眼眶、右眼眶、鼻骨边界轮廓的长度和中心坐标确定阈值。

所述步骤1.2进一步包括：

步骤1.2.1：计算颅骨模型各顶点p_i的法矢n_i＝(x_i,y_i,z_i)和高斯曲率gauss_i，通过高斯映射将颅骨顶点映射到单位球上，映射后的顶点坐标为

步骤1.2.2：颅骨模型的上颌、下颌、颧骨区域的几何形状复杂，曲面几何形状变化大，通过分析顶点的曲率和法矢法向完成上述区域的分割。选择高斯曲率gauss_i最大的顶点作为种子点，采用动态区域增长算法合并邻接顶点，直到遍历所有顶点后停止，区域增长的条件决定了特征提取的准确性，其区域合并条件定义如下：①||p'_i-p'_j||＜δ，高斯球上两个相邻顶点p'_i和p'_j间的距离描述了两个法矢方向的差异，其阈值为δ；②||gauss_i-gauss_j||＜ε，两个相邻顶点的高斯曲率gauss_i和gauss_j的差异，其阈值为ε；阈值均由统计反馈动态生成。

所述步骤1.3进一步包括：

步骤1.3.1：针对由步骤1.1和步骤1.2获得的参考颅骨模型的特征点集，其包括颅骨上颌、下颌和颧骨，以及眼眶和鼻骨的孔洞边缘点集两部分，和步骤1.1和步骤1.2获得的目标颅骨模型特征点集，采用最近点迭代算法实现两个颅骨模型的刚性配准，配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性，参考颅骨模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³；

步骤1.3.2：为了快速建立参考颅骨模型和目标颅骨模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考颅骨模型与目标颅骨模型的非刚性配准，进而将最近点作为对应点。参考颅骨模型点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应将作为步骤1.3.3中对应点误差能量项中的初始对应关系；

步骤1.3.3：定义具有保刚性的能量函数，计算参考颅骨模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考颅骨和目标颅骨的非刚性配准。该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值。具体而言，对应点误差能量项为w_i为权值，参考模型与目标模型的初始对应关系由步骤1.3.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1，1，1，1)为对角矩阵；特征误差能量项为即由步骤1.3.1产生的对应点集。

进一步地，所述步骤2包括以下具体步骤：

步骤2.1：面貌模型表面特征线提取，计算面貌三维模型每个顶点的法线和高斯曲率，统计高斯曲率值局部最大且邻接顶点法矢夹角差异较大的顶点，该顶点作为特征点；面貌五官模型耳朵、鼻子、嘴、眼睛自动分割，以顶点高斯曲率局部最大的顶点作为种子顶点，基于高斯映射和动态区域增长算法实现分割；

步骤2.2：从面貌模型数据集中选择两个模型，一个作为参考面貌，另一个作为目标面貌，针对步骤2.1生成的两个面貌模型特征，利用最近点迭代算法实现参考面貌向目标面貌的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个面貌模型的非刚性配准；

步骤2.3：变形后的参考模型的每个顶点与目标模型的最近点即为当前顶点的对应点，记录其顶点序号。针对原始参考面貌模型，生成与目标面貌点云对应的参考面貌顶点坐标；

步骤2.4：从面貌数据集中选择其他模型作为参考面貌，重复步骤2.1-步骤2.3，直到所有面貌均已建立对应关系时停止。

所述步骤2.2进一步包括：

步骤2.2.1：针对由步骤2.1获得的参考面貌模型特征点集和目标面貌模型特征点集，采用最近点迭代算法实现两个面貌模型的刚性配准，配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性，参考面貌模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标面貌对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³；

步骤2.2.2：为了快速建立参考面貌模型和目标面貌模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考面貌模型与目标面貌模型的非刚性配准，进而将最近点作为对应点，参考面貌模型的点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标面貌的对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应作为步骤2.2.3中对应点误差能量项的初始对应关系；

步骤2.2.3：定义具有保刚性的能量函数，计算参考面貌模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考面貌和目标面貌的非刚性配准。该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值。具体而言，对应点误差能量项为w_i为权值，参考模型与目前模型的初始对应关系由步骤2.2.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵；特征误差能量项为即由步骤2.2.1产生的对应点集。

进一步地，所述步骤3包括以下具体步骤：

步骤3.1：利用主成分分析方法，降维表示颅骨稠密点云，计算特征值和特征向量。利用主成分分析方法，降维表示面貌稠密点云，计算特征值和特征向量。

步骤3.2：为了观察各主成分对模型几何形状的影响，针对颅骨模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₁δ₁，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₁为该主成分系数的方差，同理，针对面貌模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₂δ₂，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₂为该主成分系数的方差；

步骤3.3：为了观察基于主成分表示的颅骨和面貌模型间的相互关系。每次分别取相同贡献率的颅骨主成分和面貌主成分，将当前选择的颅骨主成分系数和对应的面貌主成分系数分别设置为给定值如-1.0、-0.8、-0.6、-0.4、-0.2、0.2、0.4、0.6、0.8、1.0，其他主成分系数值为0，显示颅骨模型和对应的面貌模型；

步骤3.4：计算颅骨主成分系数和对应的面貌主成分系数的相关性，判断其是否近似满足线性关系；如果满足线性相关，则采用步骤4.4中的最小二乘法进行形态关系学习。

进一步地，所述步骤4包括以下具体步骤：

步骤4.1：针对数据集中的每个样本，计算每个顶点的软组织厚度值；针对所有样本，计算每个顶点的软组织厚度均值和方差，利用改进的K均值聚类算法，按软组织厚度进行聚类如分为四类；聚类过程中首先均匀采样设定聚类中心和较小粒度的聚类条件，然后完成初始分类并建立各分类间的邻接关系无向图，最后以含有顶点数较多的分类为中心通过合并邻接分类完成指定数量的聚类；

步骤4.2：针对目标颅骨和目标面貌模型，依据每个顶点对应的软组织厚度分类，将颅骨顶点和面貌顶点进行分区，实现基于软组织厚度的颅面分区；

步骤4.3：针对颅骨各个分区点云，计算主成分系数和特征向量；针对面貌各个分区点云，计算主成分系数和特征向量；

步骤4.4：针对每个分区数据集，设Skull_l×p＝[α_1,p,α_2,p,...,α_l,p]和Face_l×q＝[b_1,q,b_2,q,...,b_l,q]分别为该分区中每个样本颅骨的主成分和对应面貌的主成分，则颅骨和面貌之间的形态关系M＝{M_i,i＝1,2,L,k}能表示为M_i＝argmin||Skull×M_i-Face||²+λ²||M_i||²。利用最小二乘法求解M得M_i＝(Skull^T·Skull+λI)^-1·Skull^T·Face，其中λ为权值，I为单位矩阵，则颅骨和面貌间的形态关系表示为M＝{(M₁,M₂,L,M_k)}。

进一步地，所述步骤5包括以下具体步骤：

步骤5.1：针对待复原的未知身源颅骨，利用步骤1.1-步骤1.5实现未知身源颅骨模型与目标颅骨模型的非刚性配准，建立顶点间的对应，实现未知身源颅骨的分区和顶点对应；

步骤5.2：针对未知身源颅骨的每个分区，计算每个分区对应的主成分(skull₁,skull₂,L,skull_k)。依据颅面形态关系M＝{(M₁,M₂,L,M_k)}，计算每个颅骨分区主成分对应的面貌分区主成分(Face₁,Face₂,L,Face_k)。根据计算获得每个面貌分区主成分及其对应的特征向量，计算面貌分区复原结果。

步骤5.3：建立边缘约束的能量方程，求解各分区中每个顶点对应的仿射变换，实现面貌分区复原结果的光滑融合；

所述步骤5.3进一步包括：

步骤5.3.1：对每个分区点云进行三角剖分，建立顶点间的邻接关系。统计每个分区的顶点数量，并将顶点数量最多的分区作为目标模型，其他分区作为参考模型并向该分区变形；

步骤5.3.2：提取各分区复原结果的边缘轮廓点集。计算目标模型与每个参考模型边缘轮廓间的对应关系，对应的参考模型的边缘点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标模型对应的点集记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³；

步骤5.3.3：定义具有保刚性的能量函数，计算参考模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现各分区融合。该能量函数由边缘误差项E_edge(X)和局部刚性能量项E_s(X)两部分构成：E(X)＝E_edge(X)+αE_s(X)，其中α为权值。具体而言，边缘误差能量项为w_i为权值，其对应关系在步骤5.3.2中生成；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵。

本发明相对于现有技术具有如下的优点和积极效果：

1本发明通过基于高斯映射和区域增长结合的方法，自动提取颅骨形状变化大的区域作为特征，将自动提取的对应轮廓边缘和特征作为特征点集，建立了隐式函数和顶点能量约束最优结合的非刚性配准方法，提高了颅骨模型配准的准确性。本发明由于能够自动提取特征点集，能有效克服配准过程中需要人工参与标定特征点的不足。通过建立非刚性能量优化和顶点对应关系自动寻找的迭代机制，在非刚性配准过程中实现顶点对应关系的建立，提高了颅骨模型顶点间对应关系寻找的准确性。

2本发明能够自动提取面貌特征点集，通过建立非刚性能量优化和顶点对应关系自动确定的迭代机制，提高了面貌顶点间对应关系寻找的准确性，有效克服配准过程中需要人工参与标定特征点的不足。

3本发明通过合理选择和设置颅骨主成分和面貌主成分，能够可视地显示颅骨几何形态、面貌几何形态与主成分之间的关系，发现了颅骨主成分和面貌主成分的几何意义。通过主成分之间的相关性检测，进一步解释了基于主成分和线性回归的颅面形态关系分析方法的几何意义。

4本发明通过统计面部稠密顶点的软组织厚度，实现了基于软组织分布的颅面分区。相比较于基于位置的五官分区方法，使得各分区内的颅面形态关系更为一致，提高了利用最小二乘回归方法进行颅面形态关系定量表示的准确性，克服了过拟合和拟合不足的问题。

5本发明通过建立边缘约束的能量方程，提高了分区复原结果融合的准确性。

本发明提供一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，是对现有颅面形态分析方法进行了创新和扩展。通过自动确定对应特征和隐式函数和顶点能量约束最优结合的非刚性配准方法，提高了颅骨点云和面貌点云稠密对应的准确性；通过分析面部稠密软组织分布，建立了基于软组织分布的颅面分区模型，提高了颅面形态关系表示的准确性。通过建立边缘约束的能量方程，提高了分区复原结果融合的准确性。最终应用于无身源颅骨的面貌复原。

附图说明

图1本发明所述一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法的流程图。

具体实施方式

以下结合附图对本发明的实施例进行详细描述，但本发明并不仅仅限于这些实施例。本发明涵盖任何在本发明的精髓和范围上做的替代、修改、等效方法以及方案。

如图1所示，本发明一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法包括如下步骤：

1、颅骨稠密对应点云；

2、面貌稠密对应点云；

3、颅面形态关系可视分析；

4、基于软组织分区的颅面形态关系表示；

5、未知身源颅骨的面貌复原。

如图1所示，在所述步骤1颅骨稠密对应点云和步骤2面貌稠密点云完成后，根据颅面形态分析的方式进行步骤3颅面形态关系可视分析和步骤4基于软组织分区的颅面形态关系表示。最后依据颅面形态关系完成步骤5未知身源颅骨的面貌复原。

所述步骤1包括以下具体步骤：

步骤1.1：针对颅骨三维网格模型，定义内点和边界点，实现颅骨孔洞边缘的自动识别；并通过孔洞边界长度、孔洞形状和位置，分析确定孔洞之间的对应关系；

步骤1.2：基于高斯映射和动态区域增长算法，实现颅骨几何形状变化较大区域包括上颌、下颌、颧骨等的自动分割；

步骤1.3：从颅骨模型数据集中选择选择两个模型，一个作为参考颅骨，另一个作为目标颅骨，依据步骤1.1和步骤1.2产生的两个颅骨模型特征，利用最近点迭代算法实现参考颅骨向目标颅骨的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个颅骨模型的非刚性配准；

步骤1.4：计算目标模型的每个顶点与变形后的参考模型的最近点，记录其顶点序号，建立目标模型与参考模型顶点的对应关系；

步骤1.5：依据步骤1.4计算获得的顶点序号，针对原始参考颅骨模型，生成与目标颅骨点云对应的参考颅骨顶点坐标；

步骤1.6：从颅骨模型数据集中选择其他模型作为参考颅骨，重复步骤1.3-步骤1.5，直到遍历所有模型为止，从而建立颅骨模型间的顶点对应关系，即所有颅骨模型具有相同的顶点个数且对应顶点具有近似的解剖学位置。

所述步骤1.1进一步包括：

步骤1.1.1：定义颅骨三维网格模型为skull＝{P,E}，其中P＝{p₁,p₂,...,p_n},p_i＝(x_i,y_i,z_i)∈R³表示n个颅骨顶点，E＝{e_k＝(p_i,p_j),k＝1,2,3,...,m}表示m条边。孔洞识别前首选对模型进行检验，判断其是否满足流行，并记录不满足条件的边序号。孔洞识别过程中，首先遍历模型中各个顶点pi，确定其对应的邻接顶点集合Adj_pi，如果Adj_pi中的点通过模型的边e_k直接连成封闭的多边形，则p_i为内点，否则为边界点。然后，以任意边界点p_j为起点，根据e_k寻找下一个边界点，直到遍历完所有边界点。最终则会获得c条封闭边界，记为boundary_i＝{p_j},i＝1,2,3,...,c。

步骤1.1.2：计算每条边界的长度、中心，通过分析边界长度和中心位置即能识别左眼眶、右眼眶和鼻骨三个轮廓。能通过统计训练集中样本的左眼眶、右眼眶、鼻骨边界轮廓的长度和中心坐标确定阈值。

所述步骤1.2进一步包括：

步骤1.2.1：计算颅骨模型各顶点p_i的法矢n_i＝(x_i,y_i,z_i)和高斯曲率gauss_i，通过高斯映射将颅骨顶点映射到单位球上，映射后的顶点坐标为

步骤1.2.2：颅骨模型的上颌、下颌、颧骨区域的几何形状复杂，曲面几何形状变化大，通过分析顶点的曲率和法矢法向完成上述区域的分割。选择高斯曲率gauss_i最大的顶点作为种子点，采用动态区域增长算法合并邻接顶点，直到遍历所有顶点后停止，区域增长的条件决定了特征提取的准确性，其区域合并条件定义如下：①||p'_i-p'_j||＜δ，高斯球上两个相邻顶点p'_i和p'_j间的距离描述了两个法矢方向的差异，其阈值为δ；②||gauss_i-gauss_j||＜ε，两个相邻顶点的高斯曲率gauss_i和gauss_j的差异，其阈值为ε；阈值均由统计反馈动态生成。

所述步骤1.3进一步包括：

步骤1.3.1：针对由步骤1.1和步骤1.2获得的参考颅骨模型的特征点集，其包括颅骨上颌、下颌和颧骨，以及眼眶和鼻骨的孔洞边缘点集两部分，和步骤1.1和步骤1.2获得的目标颅骨模型特征点集，采用最近点迭代算法实现两个颅骨模型的刚性配准，配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性，参考颅骨模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³；

步骤1.3.2：为了快速建立参考颅骨模型和目标颅骨模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考颅骨模型与目标颅骨模型的非刚性配准，进而将最近点作为对应点，参考颅骨模型点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应将作为步骤1.3.3中对应点误差能量项中的初始对应关系；

步骤1.3.3：定义具有保刚性的能量函数，计算参考颅骨模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考颅骨和目标颅骨的非刚性配准。该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值。具体而言，对应点误差能量项为w_i为权值，参考模型与目标模型的初始对应关系由步骤1.3.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵；特征误差能量项为即由步骤1.3.1产生的对应点集。

进一步地，所述步骤2包括以下具体步骤：

步骤2.1：面貌模型表面特征线提取，计算面貌三维模型每个顶点的法线和高斯曲率，统计高斯曲率值局部最大且邻接顶点法矢夹角差异较大的顶点，该顶点作为特征点；面貌五官模型耳朵、鼻子、嘴、眼睛自动分割，以顶点高斯曲率局部最大的顶点作为种子顶点，基于高斯映射和动态区域增长算法实现分割；

步骤2.2：从面貌模型数据集中选择两个模型，一个作为参考面貌，另一个作为目标面貌，针对步骤2.1生成的两个面貌模型特征，利用最近点迭代算法实现参考面貌向目标面貌的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个面貌模型的非刚性配准；

步骤2.3：变形后的参考模型的每个顶点与目标模型的最近点即为当前顶点的对应点，记录其顶点序号。针对原始参考面貌模型，生成与目标面貌点云对应的参考面貌顶点坐标；

步骤2.4：从面貌数据集中选择其他模型作为参考面貌，重复步骤2.1-步骤2.3，直到所有面貌均已建立对应关系时停止。

所述步骤2.2进一步包括：

步骤2.2.1：针对由步骤2.1获得的参考面貌模型特征点集和目标面貌模型特征点集，采用最近点迭代算法实现两个面貌模型的刚性配准。配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性。参考面貌模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标面貌对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³。

步骤2.2.2：为了快速建立参考面貌模型和目标面貌模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考面貌模型与目标面貌模型的非刚性配准，进而将最近点作为对应点。参考面貌模型的点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标面貌的对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应作为步骤2.2.3中对应点误差能量项的初始对应关系。

步骤2.2.3：定义具有保刚性的能量函数，计算参考面貌模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考面貌和目标面貌的非刚性配准。该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值。具体而言，对应点误差能量项为w_i为权值，参考模型与目前模型的初始对应关系由步骤2.2.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵；特征误差能量项为即由步骤2.2.1产生的对应点集。

进一步地，所述步骤3包括以下具体步骤：

步骤3.1：利用主成分分析方法，降维表示颅骨稠密点云，计算特征值和特征向量。利用主成分分析方法，降维表示面貌稠密点云，计算特征值和特征向量。

步骤3.2：为了观察各主成分对模型几何形状的影响，针对颅骨模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₁δ₁，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₁为该主成分系数的方差，同理，针对面貌模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₂δ₂，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₂为该主成分系数的方差；

步骤3.3：为了观察基于主成分表示的颅骨和面貌模型间的相互关系。每次分别取相同贡献率的颅骨主成分和面貌主成分，将当前选择的颅骨主成分系数和对应的面貌主成分系数分别设置为给定值如-1.0、-0.8、-0.6、-0.4、-0.2、0.2、0.4、0.6、0.8、1.0，其他主成分系数值为0，显示颅骨模型和对应的面貌模型；

步骤3.4：计算颅骨主成分系数和对应的面貌主成分系数的相关性，判断其是否近似满足线性关系。如果满足线性相关，则采用步骤4.4中的最小二乘法进行形态关系学习。

进一步地，所述步骤4包括以下具体步骤：

步骤4.1：针对数据集中的每个样本，计算每个顶点的软组织厚度值；针对所有样本，计算每个顶点的软组织厚度均值和方差，利用改进的K均值聚类算法，按软组织厚度进行聚类如分为四类；聚类过程中首先均匀采样设定聚类中心和较小粒度的聚类条件，然后完成初始分类并建立各分类间的邻接关系无向图，最后以含有顶点数较多的分类为中心通过合并邻接分类完成指定数量的聚类；

步骤4.2：针对目标颅骨和目标面貌模型，依据每个顶点对应的软组织厚度分类，将颅骨顶点和面貌顶点进行分区，实现基于软组织厚度的颅面分区；

步骤4.3：针对颅骨各个分区点云，计算主成分系数和特征向量；针对面貌各个分区点云，计算主成分系数和特征向量；

步骤4.4：针对每个分区数据集，设Skull_l×p＝[α_1,p,α_2,p,...,α_l,p]和Face_l×q＝[b_1,q,b_2,q,...,b_l,q]分别为该分区中每个样本颅骨的主成分和对应面貌的主成分，则颅骨和面貌之间的形态关系M＝{M_i,i＝1,2,L,k}能表示为M_i＝argmin||Skull×M_i-Face||²+λ²||M_i||²。利用最小二乘法求解M得M_i＝(Skull^T·Skull+λI)^-1·Skull^T·Face，其中λ为权值，I为单位矩阵，则颅骨和面貌间的形态关系表示为M＝{(M₁,M₂,L,M_k)}。

进一步地，所述步骤5包括以下具体步骤：

步骤5.1：针对待复原的未知身源颅骨，利用步骤1.1-步骤1.5实现未知身源颅骨模型与目标颅骨模型的非刚性配准，建立顶点间的对应，实现未知身源颅骨的分区和顶点对应；

步骤5.2：针对未知身源颅骨的每个分区，计算每个分区对应的主成分(skull₁,skull₂,L,skull_k)。依据颅面形态关系M＝{(M₁,M₂,L,M_k)}，计算每个颅骨分区主成分对应的面貌分区主成分(Face₁,Face₂,L,Face_k)。根据计算获得每个面貌分区主成分及其对应的特征向量，计算面貌分区复原结果。

步骤5.3：建立边缘约束的能量方程，求解各分区中每个顶点对应的仿射变换，实现面貌分区复原结果的光滑融合；

所述步骤5.3进一步包括：

步骤5.3.1：对每个分区点云进行三角剖分，建立顶点间的邻接关系。统计每个分区的顶点数量，并将顶点数量最多的分区作为目标模型，其他分区作为参考模型并向该分区变形；

步骤5.3.2：提取各分区复原结果的边缘轮廓点集。计算目标模型与每个参考模型边缘轮廓间的对应关系，对应的参考模型的边缘点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标模型对应的点集记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³；

步骤5.3.3：定义具有保刚性的能量函数，计算参考模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现各分区融合。该能量函数由边缘误差项E_edge(X)和局部刚性能量项E_s(X)两部分构成：E(X)＝E_edge(X)+αE_s(X)，其中α为权值。具体而言，边缘误差能量项为w_i为权值，其对应关系在步骤5.3.2中生成；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵。

在本发明的实施例中：首先，针对人体头部CT影像数据，通过图像分割和边缘轮廓跟踪技术，提取颅骨和面貌的外表面，进而将获得的颅骨点云和面貌点云模型变换到法兰克福坐标系下。然后，通过提取颅骨点云模型的特征、孔洞边缘，基于能量优化的非刚性配准方法，建立颅骨稠密对应点云。特征自动提取过程中，为了消除特征误对应的问题，需要考虑曲率一致性、法线一致性、距离的约束条件并结合RANSAC方法，提高特征对应的准确性。面貌稠密点云对应过程中也需要确保特征对应的一致性，消除误匹配。最后，在颅面形态关系可视分析过程中，为了保证分析的有效性，需要选择具有相同贡献率的颅骨主成分和面貌主成分。可视分析过程中，用户应该同时设置颅骨主成分系数和面貌主成分系数值，比较模型的长度、宽度和角度的变化。

基于软组织分区的颅面形态关系定量分析过程中必须针对已经建立好稠密对应关系的颅骨点云和面貌点云。首先计算每个顶点的软组织厚度和整个样本集的软组织厚度均值，每个分区的软组织均值可以通过统计脸颊处、额头、下颌、颧骨的软组织厚度均值获得，从而提高分区结果的有效性。如果聚类后的分区过小，需要与相邻分区进行合并。针对每个分区利用最小二乘拟合颅面形态关系时，需要引入扰动项，克服过拟合和拟合不足的问题。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明公开的范围内，能够轻易想到的变化或替换，都应涵盖在本发明权利要求的保护范围内。

Claims (5)

1.一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，其特征在于，包括以下步骤：

1：颅骨稠密对应点云；

步骤1.1：针对颅骨三维网格模型，定义内点和边界点，实现颅骨孔洞边缘的自动识别；并通过孔洞边界长度、孔洞形状和位置，分析确定孔洞之间的对应关系；

步骤1.2：基于高斯映射和动态区域增长算法，实现颅骨几何形状变化较大区域包括上颌、下颌、颧骨等的自动分割；

步骤1.3：从颅骨模型数据集中选择选择两个模型，一个作为参考颅骨，另一个作为目标颅骨，依据步骤1.1和步骤1.2产生的两个颅骨模型特征，利用最近点迭代算法实现参考颅骨向目标颅骨的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个颅骨模型的非刚性配准；

步骤1.4：计算目标模型的每个顶点与变形后的参考模型的最近点，记录其顶点序号，建立目标模型与参考模型顶点的对应关系；

步骤1.5：依据步骤1.4计算获得的顶点序号，针对原始参考颅骨模型，生成与目标颅骨点云对应的参考颅骨顶点坐标；

步骤1.6：从颅骨模型数据集中选择其他模型作为参考颅骨，重复步骤1.3-步骤1.5，直到遍历所有模型为止，从而建立颅骨模型间的顶点对应关系，即所有颅骨模型具有相同的顶点个数且对应顶点具有近似的解剖学位置；

2：面貌稠密对应点云；

步骤2.1：面貌模型表面特征线提取，计算面貌三维模型每个顶点的法线和高斯曲率，统计高斯曲率值局部最大且邻接顶点法矢夹角差异较大的顶点，该顶点作为特征点；面貌五官模型耳朵、鼻子、嘴、眼睛自动分割，以顶点高斯曲率局部最大的顶点作为种子顶点，基于高斯映射和动态区域增长算法实现分割；

步骤2.2：从面貌模型数据集中选择两个模型，一个作为参考面貌，另一个作为目标面貌，针对步骤2.1生成的两个面貌模型特征，利用最近点迭代算法实现参考面貌向目标面貌的刚性配准，进一步提出通过基于隐式函数和顶点能量约束最优结合的方法，实现两个面貌模型的非刚性配准；

步骤2.3：变形后的参考模型的每个顶点与目标模型的最近点即为当前顶点的对应点，记录其顶点序号，针对原始参考面貌模型，生成与目标面貌点云对应的参考面貌顶点坐标；

步骤2.4：从面貌数据集中选择其他模型作为参考面貌，重复步骤2.1-步骤2.3，直到所有面貌均已建立对应关系时停止；

3：颅面形态关系可视分析；

步骤3.1：利用主成分分析方法，降维表示颅骨稠密点云，计算特征值和特征向量，利用主成分分析方法，降维表示面貌稠密点云，计算特征值和特征向量；

步骤3.2：为了观察各主成分对模型几何形状的影响，针对颅骨模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₁δ₁，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₁为该主成分系数的方差，同理，针对面貌模型，从第一个主成分开始，将其对应的主成分系数设置为给定值value＝3·λ₂δ₂，其中λ₁＝{-1.0,-0.8,-0.6,-0.4,-0.2,0.2,0.4,0.6,0.8,1.0}，δ₂为该主成分系数的方差；

步骤3.3：为了观察基于主成分表示的颅骨和面貌模型间的相互关系，每次分别取近似相同贡献率的颅骨主成分和面貌主成分，将当前选择的颅骨主成分系数和对应的面貌主成分系数分别设置为给定值，其他主成分系数值为0，显示颅骨模型和对应的面貌模型；

步骤3.4：计算颅骨主成分系数和对应的面貌主成分系数的相关性，判断其是否近似满足线性关系；如果满足线性相关，则采用步骤4.4中的最小二乘法进行形态关系学习；

4：基于软组织分区的颅面形态关系表示；

步骤4.1：针对数据集中的每个样本，计算每个顶点的软组织厚度值；针对所有样本，计算每个顶点的软组织厚度均值和方差，利用改进的K均值聚类算法，按软组织厚度进行聚类分为四类；聚类过程中首先均匀采样设定聚类中心和较小粒度的聚类条件，然后完成初始分类并建立各分类间的邻接关系无向图，最后以含有顶点数较多的分类为中心通过合并邻接分类完成指定数量的聚类；

步骤4.2：针对目标颅骨和目标面貌模型，依据每个顶点对应的软组织厚度分类，将颅骨顶点和面貌顶点进行分区，实现基于软组织厚度的颅面分区；

步骤4.3：针对颅骨各个分区点云，计算主成分系数和特征向量；针对面貌各个分区点云，计算主成分系数和特征向量；

步骤4.4：针对每个分区数据集，设Skull_l×p＝[α_1,p,α_2,p,...,α_l,p]和Face_l×q＝[b_1,q,b_2,q,...,b_l,q]分别为该分区中每个样本颅骨的主成分和对应面貌的主成分，则颅骨和面貌之间的形态关系M＝{M_i,i＝1,2,L,k}能表示为M_i＝argmin||Skull×M_i-Face||²+λ²||M_i||²，利用最小二乘法求解M得M_i＝(Skull^T·Skull+λI)^-1·Skull^T·Face，其中λ为权值，I为单位矩阵，则颅骨和面貌间的形态关系表示为M＝{(M₁,M₂,L,M_k)}；

5：未知身源颅骨的面貌复原。

2.根据权利要求1所述的一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，其特征在于，所述步骤5包括以下具体步骤：

步骤5.1：针对待复原的未知身源颅骨，利用步骤1.1-步骤1.5实现未知身源颅骨模型与目标颅骨模型的非刚性配准，建立顶点间的对应，实现未知身源颅骨的分区和顶点对应；

步骤5.2：针对未知身源颅骨的每个分区，计算每个分区对应的主成分(skull₁,skull₂,L,skull_k)；依据颅面形态关系M＝{(M₁,M₂,L,M_k)}，计算每个颅骨分区主成分对应的面貌分区主成分(Face₁,Face₂,L,Face_k)；根据计算获得每个面貌分区主成分及其对应的特征向量，计算面貌分区复原结果；

步骤5.3：建立边缘约束的能量方程，求解各分区中每个顶点对应的仿射变换，实现面貌分区复原结果的光滑融合。

3.根据权利要求1所述的一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，其特征在于，所述步骤1.1进一步包括：

步骤1.1.1：定义颅骨三维网格模型为skull＝{P,E}，其中P＝{p₁,p₂,...,p_n},p_i＝(x_i,y_i,z_i)∈R³表示n个颅骨顶点，E＝{e_k＝(p_i,p_j),k＝1,2,3,...,m}表示m条边，孔洞识别前首选对模型进行检验，判断其是否满足条件，并记录不满足条件的边序号，孔洞识别过程中，首先遍历模型中各个顶点p_i，确定其对应的邻接顶点集合Adj_pi，如果Adj_pi中的点通过模型的边e_k直接连成封闭的多边形，则p_i为内点，否则为边界点，然后，以任意边界点p_j为起点，根据e_k寻找下一个边界点，直到遍历完所有边界点，最终则会获得c条封闭边界，记为boundary_i＝{p_j},i＝1,2,3,...,c；

步骤1.1.2：计算每条边界的长度、中心，通过分析边界长度和中心位置即能识别左眼眶、右眼眶和鼻骨三个轮廓，通过统计训练集中样本的左眼眶、右眼眶、鼻骨边界轮廓的长度和中心坐标确定阈值；

所述步骤1.2进一步包括：

步骤1.2.1：计算颅骨模型各顶点p_i的法矢n_i＝(x_i,y_i,z_i)和高斯曲率gauss_i，通过高斯映射将颅骨顶点映射到单位球上，映射后的顶点坐标为

步骤1.2.2：颅骨模型的上颌、下颌、颧骨区域的几何形状复杂，曲面几何形状变化大，通过分析顶点的曲率和法矢法向完成上述区域的分割；选择高斯曲率gauss_i最大的顶点作为种子点，采用动态区域增长算法合并邻接顶点，直到遍历所有顶点后停止，区域增长的条件决定了特征提取的准确性，其区域合并条件定义如下：①||p'_i-p'_j||＜δ，高斯球上两个相邻顶点p'_i和p'_j间的距离描述了两个法矢方向的差异，其阈值为δ；②||gauss_i-gauss_j||＜ε，两个相邻顶点的高斯曲率gauss_i和gauss_j的差异，其阈值为ε；阈值均由统计反馈动态生成；

所述步骤1.3进一步包括：

步骤1.3.1：针对由步骤1.1和步骤1.2获得的参考颅骨模型的特征点集，其包括颅骨上颌、下颌和颧骨，以及眼眶和鼻骨的孔洞边缘点集两部分，和步骤1.1和步骤1.2获得的目标颅骨模型特征点集，采用最近点迭代算法实现两个颅骨模型的刚性配准，配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性，参考颅骨模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³；

步骤1.3.2：为了快速建立参考颅骨模型和目标颅骨模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考颅骨模型与目标颅骨模型的非刚性配准，进而将最近点作为对应点，参考颅骨模型点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标颅骨对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应将作为步骤1.3.3中对应点误差能量项中的初始对应关系；

步骤1.3.3：定义具有保刚性的能量函数，计算参考颅骨模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考颅骨和目标颅骨的非刚性配准，该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值，对应点误差能量项为w_i为权值，参考模型与目标模型的初始对应关系由步骤1.3.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵；特征误差能量项为即由步骤1.3.1产生的对应点集。

4.根据权利要求1所述的一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，其特征在于，所述步骤2.2进一步包括：

步骤2.2.1：针对由步骤2.1获得的参考面貌模型特征点集和目标面貌模型特征点集，采用最近点迭代算法实现两个面貌模型的刚性配准，配准过程中采用随机抽样一致算法去除特征点集的误对应，建立顶点间的准确对应关系，提高配准结果的准确性，参考面貌模型的特征点集记为S＝{s_i},s_i＝(x_i,y_i,z_i)∈R³，目标面貌对应的特征点集记为Q＝{q_i},q_i＝(x_i,y_i,z_i)∈R³；

步骤2.2.2：为了快速建立参考面貌模型和目标面貌模型顶点间的对应关系，采用径向基函数和带有紧支撑的径向基函数，实现参考面貌模型与目标面貌模型的非刚性配准，进而将最近点作为对应点，参考面貌模型的点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标面貌的对应点记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³，该对应作为步骤2.2.3中对应点误差能量项的初始对应关系；

步骤2.2.3：定义具有保刚性的能量函数，计算参考面貌模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现参考面貌和目标面貌的非刚性配准，该能量函数由对应点误差能量项E_d(X)、特征误差项E_l(X)和局部刚性能量项E_s(X)三部分构成：E(X)＝E_d(X)+αE_s(X)+βE_l(X)，其中α和β为权值，对应点误差能量项为w_i为权值，参考模型与目前模型的初始对应关系由步骤2.2.2获得，此后每次代过程中通过搜寻最近点确定对应关系(v_i,u_i)；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵；特征误差能量项为即由步骤2.2.1产生的对应点集。

5.根据权利要求2所述的一种基于颅面稠密对应点云的颅面形态分析及面貌复原方法，其特征在于，所述步骤5.3进一步包括：

步骤5.3.1：对每个分区点云进行三角剖分，建立顶点间的邻接关系，统计每个分区的顶点数量，并将顶点数量最多的分区作为目标模型，其他分区作为参考模型并向该分区变形；

步骤5.3.2：提取各分区复原结果的边缘轮廓点集；计算目标模型与每个参考模型边缘轮廓间的对应关系，对应的参考模型的边缘点集记为V＝{v_i},v_i＝(x_i,y_i,z_i)∈R³，目标模型对应的点集记为U＝{u_i},u_i＝(x_i,y_i,z_i)∈R³；

步骤5.3.3：定义具有保刚性的能量函数，计算参考模型每个顶点p_i对应的仿射变换进而根据仿射变换矩阵X＝[x₁ x₂ … x_n]^T实现各分区融合；该能量函数由边缘误差项E_edge(X)和局部刚性能量项E_s(X)两部分构成：E(X)＝E_edge(X)+αE_s(X)，其中α为权值，边缘误差能量项为w_i为权值，其对应关系在步骤5.3.2中生成；局部刚性能量项为其中F表示由目标模型的顶点-边组成的邻接矩阵,为克罗内克乘积算子，G＝diag(1,1,1,1)为对角矩阵。