基于赋权形态学分析的三维面部对称参考平面构建方法

朱玉佳，赵一姣，郑盛文，温奥楠，傅湘玲△，王 勇△

(1. 北京大学口腔医学院·口腔医院，口腔医学数字化研究中心 国家口腔疾病临床医学研究中心 口腔数字化医疗技术和材料国家工程实验室 口腔数字医学北京市重点实验室，北京 100081； 2. 北京大学口腔医学院·口腔医院口腔修复教研室，北京 100081； 3. 北京邮电大学计算机学院(国家示范性软件学院)，北京 100876； 4. 北京邮电大学可信分布式计算与服务教育部重点实验室，北京 100876)

面部对称性是人类学面部吸引力及口腔美学的关键因素，由于发育过程中的生物因素和环境因素影响，几乎没有绝对对称的人脸[1-4]。面部不对称畸形的临床表现具有多样化和复杂性的特性，且多见于面下1/3[5-6]。随着数字化诊疗技术的发展，三维面部对称性分析成为口腔正颌、正畸矫治制定相应治疗计划的基础环节，也是口腔美学修复设计的重要环节[7-9]，对称参考平面的正确与否对诊疗结果有非常大的影响。三维面部对称参考平面(symmetry reference plane，SRP，也称正中矢状平面)的构建是三维面部对称性分析的前提和关键。

传统方法常基于面部三维数字模型上重要解剖标志点或结合自然头位构建对称参考平面，该方法的核心在于面部解剖标志点的选择[10-11]。既往研究对标志点的选择标准不一，较难实现适用于各类型面部畸形的共性方法，因此，近年来基于三维面部模型本体及其镜像模型重叠关联构建对称参考平面的方法(简称本体-镜像关联法)得到越来越多的关注[12]。本体-镜像关联法的原理是对三维面部模型(本体模型)及其镜像模型的几何形态进行最优三维重叠匹配，重叠后的联合模型称为关联模型，该关联模型理论上为一对称模型，通过分析关联模型的对称平面，从而确定原始三维面部模型(本体模型)的对称参考平面。

本体-镜像关联法的核心环节是三维面部本体与其镜像间的最优重叠算法，现有研究主要采用迭代最近点(iterative closest point，ICP)算法和普氏分析(Procrustes analysis，PA)算法[13-15]。ICP算法在本领域应用较早，它不依赖于解剖标志点，软件可自动化完成本体和镜像三维面部模型间的全数据最优匹配重叠，但对复杂面部畸形的情况，临床还需辅助人工对面部非畸形区域数据进行筛选，称为“区域ICP算法”[16-17]。PA算法则基于面部重要解剖标志点，在保证三维面部本体和镜像模型上同名解剖标志点一一对应关系的基础上，实现本体和镜像标志点集间平均距离最小的重叠，从而获得本体和镜像模型的最优重叠位置。PA算法对面部重要解剖标志点的关注较符合口腔临床的诊断经验和习惯，并被证实具有较好的可靠性[14]，但也同样存在对复杂面部畸形的适合性欠佳问题[18]，分析其原因在于现有PA算法不具备对面部个性化畸形特征(不同程度的不对称)进行量化评估和权重分配的机制。

本研究在既往PA算法研究的基础上，拟通过形态学研究方法中的欧氏距离矩阵分析(Euclidean distance matrix analysis，EDMA)方法，对面部重要解剖标志点进行量化评估，并基于该量化评估结果对面部解剖标志点进行量化赋权，从而建立一种基于赋权普氏分析算法(weighted Procrustes analysis，WPA)的三维面部对称参考平面自动构建方法，并应用该方法初步对口腔临床较为常见的下颌偏斜类面部不对称畸形病例进行算法的临床适用性研究。

1 资料与方法

1.1 研究对象

选取就诊于北京大学口腔医院正畸科、颌面外科、修复科的下颌偏斜面部不对称畸形患者30例，纳入标准为：患者处于自然头位时，颏前点与经鼻根点垂直于双侧瞳孔连线的面中线偏差大于3 mm。本研究获得北京大学口腔医院生物医学伦理委员会批准(PKUSSIRB-202054042)，患者均对本研究内容、目的充分知情并同意。

1.2 实验设备及软件

1.3 三维面部数据的获取和处理

采集三维面相之前进行设备校准以保证能够准确取像，患者在距离照相设备135 cm处，临床医师指导患者至自然头位状态，患者采取端坐位，双眼平视前方，保持眶耳平面(Frankfort plane)与地平面平行，于面部表情自然放松时获取数据。三维面部数据能够采纳的标准为：有效显示面部轮廓、分辨率高、没有明显的移动和闭口。使用逆向工程软件Geomagic Studio 2013对采集的患者面部三维数据进行必要的处理，包括删除多余数据(如裁减头发、颈部)和修补孔洞。调整原始三维面部模型空间位姿至自然头位，使自然头位坐标系的眶耳平面与全局坐标系XZ平面重合，矢状面与全局坐标系YZ平面重合。由口腔临床高年资医师根据临床经验完成本体面部模型(Model_Org)重要解剖标志点集的提取，在全面部区域选取发际点、额点、鼻根点、鼻尖点、鼻底点等共计32个解剖标志点(中线10个、双侧各11个)，间隔1周提取3次，求坐标均值作为本体标志点集(LMK_Org)，标志点集如图1所示。计算本体标志点集的重心坐标，将本体模型平移至标志点集重心与全局坐标系原点重合，保存为OBJ文件。

Midline landmarks: glabella (Gb), trichion (Tri), pronasale (Prn), nasion (N), subnasale (Sn), labiale superius (Ls), labiale inferius (Li), sublabiale (Sl), pogonion (Pg), gnathion (Gn); bilateral landmarks: superciliary ridge (Su), exocanthion (Ex), endocanthion (En), pupil (Pu), alare (Ala), subalare (Sal), zygion (Zg), tragion (Tr), crista philtre (Cph), cheilion (Ch), gonion (Go).

1.4 标志点不对称性的量化评估

采用EDMA法量化评估三维面部解剖标志点的不对称性，并依据评估结果进行不对称性排序[19-20]。EDMA法是一种基于一系列标志点间的距离组成的矩阵来比较两个个体几何形态差异的方法。通过分别计算单一个体所有标志点之间的欧氏距离，建立欧氏距离形态矩阵，两个个体矩阵中相对应的距离线段的比值组成一个矩阵，即为两个体的形态差异矩阵(form difference matrix，FDM)。

为研究单一标志点对左、右脸形态差异的贡献大小(不对称性程度)，通过逐一删除成对的单侧标志点或中线单个标志点(i=1、2…21)计算剩余标志点集的FDM，并定义T(-i)为删除第i个标志点后FDM元素的平均值。T(-i)值越接近1，代表删除该标志点后剩余标志点集的不对称程度有所降低，即该点在标志点集中不对称性的贡献较大，可认为T(-i)最小值对应的标志点不对称程度最高，因而其不对称性排序最高。将各点的不对称性排序定义为EDMA值，按T(-i)值从大到小，EDMA值分别为1～21。

Examples were the exocanthion (Ex), zygion (Zg), cheilion (Ch), nasion (N), pronasale (Prn). The yellow line segment was the distance between the bilateral landmarks, the green line segment was the distance between the midline and bilateral landmarks, and red points inside the blue dotted line were right landmark set and left landmark set, respectively.

1.5 面部对称参考平面的构建

1.5.1经典PA算法对称参考平面的构建(实验组一) 将30个病例模型的本体标志点集(LMK_Org)中32个标志点的三维空间坐标输入基于python语言编制的经典PA算法程序，该程序首先需要对本体标志点集基于YZ平面镜像获得镜像标志点集(LMK_Mir)，计算基于经典PA算法实现本体和镜像标志点集间配对标志点的重叠效果，并对坐标变换后的镜像标志点集与本体标志点集中的对应标志点(32对)通过求解方程来拟合对称参考平面，定义为PA对称参考平面(SRP_PA)。

1.5.3区域ICP算法(参考组) 同样基于上述Model_Org模型，在Geomagic Studio软件中由高年资医师手工选取本体模型上面部对称性良好的区域，对Model_Org基于YZ平面获得镜像模型Model_Mir，根据专家选取区域进行ICP配准(保证本体模型固定、镜像模型浮动调整)，本体与镜像模型的专家定义区域重叠后，获得本体-镜像联合模型并计算对称参考平面，定义为“真值平面”(SRP_Ref)， 最终构建出30例患者面部数据的“真值平面”。同一三维面部数据基于上述PA、WPA及区域ICP算法构建的对称参考平面效果如图3所示。

1.6 数据分析

计算30个下颌偏斜三维面部模型通过PA算法和WPA算法构建的对称参考平面(SRP_PA和SRP_WPA)与真值平面(SRP_Ref)的夹角，分别记做Ang_PA和Ang_WPA，计算各算法平面角度误差的平均值和标准差。

使用SPSS 21.0软件对实验组一(PA算法)和实验组二(WPA算法)的平面角度误差进行K-S正态性检验，采用配对t检验(双侧检验)的方法对WPA算法、PA算法的平面角度误差进行统计学推断，检验水准α均取双侧0.05。

2 结果

PA算法和WPA算法的对称参考平面的角度误差(Ang_PA和Ang_WPA)均服从正态分布，PA算法组的平面角度误差为2.06°±0.86°，WPA算法组为1.53°±0.84°，配对t检验显示两组间差异有统计学意义(P<0.05)，WPA算法组的平面角度误差较小，说明30例患者的面部数据经WPA算法构建的对称参考平面与“真值平面”更为接近(表1)。

采用颏部至专家定义的“真值平面”的距离评估患者的面下部不对称畸形程度，30例不同偏斜程度患者的平面角度误差分布规律见表2。各分组数据显示，WPA算法的平面角度误差均小于PA算法，且WPA算法的角度误差均未超过2°，随患者面下部不对称畸形程度增大，两算法的角度误差也有所增大，但WPA算法的误差增幅明显小于PA算法，特别是在明显下颌偏斜组中两者差距最大，说明WPA算法比经典PA算法对下颌显著偏斜患者的适合性更好。

3 讨论

3.1 WPA算法标志点的赋权策略较符合口腔临床专家的诊疗思维

回顾三维面部对称参考平面本体-镜像关联法的既往研究，本体和镜像模型的重叠算法主要以ICP算法和PA算法为主[13-15]。ICP算法是一种完全不参考解剖标志点的算法，虽然有研究验证了ICP关联法在构建正常面型面部数据的对称参考平面时的可靠性和可重复性，但对于面部严重不对称畸形的患者，对称性不良数据会对算法造成干扰，使对称参考平面的构建效果欠佳[16]。有研究对全局ICP关联法进行了改进，通过人工选取对称性良好的面部区域进行本体和镜像重叠运算，一定程度上改进了ICP关联法的临床适合性，称之为区域ICP关联法[16-17,21-23]。区域ICP关联法引入了人为干预，虽然提高了准确性，但也降低了算法的自动化程度。

本研究采用专家筛选区域的区域ICP关联法作为“真值组”，评价本研究WPA算法的准确性，该算法目前具有较好的临床认可度[17,23]。

PA算法与ICP算法的显著区别在于PA算法更为关注面部解剖标志点的参考价值，这种算法的核心思想是基于与对称性判断密切相关的解剖标志点集的重叠，引导构建本体-镜像关联模型，从而获得对称参考平面。这种算法的原理比较符合口腔医生临床诊疗的习惯和经验，近年来也得到广泛关注。有研究证实了PA关联法在面部无明显不对称的患者中适用性良好[18]。但是，与ICP算法一样，对称性不良的PA标志点会对PA算法产生匹诺曹效应(Pinocchio effect)[24]，同样存在着对复杂面部畸形对称参考平面构建效果欠佳的问题。PA算法的改进方向之一是对PA标志点进行筛选，Gateno等[22]通过“递归PA算法”对标志点进行排序，删除明显不对称的标志点(离群值)，使用剩余的标志点进行PA运算，从而避免不良标志点的干扰，但这种算法对于复杂面部畸形的患者(绝大部分标志点对称性均不理想)，可能会剔除过多的标志点，出现局部过优化的结果倾向。PA算法的另一改进方向即为对不同标志点进行赋权，有研究使用稠密序列类标志点作为面部蒙版，将其映射到患者模型上，通过迭代筛选离群点，使用PA算法进行本体与镜像蒙版重叠，并依据每个标志点的匹配质量进行赋权[25-27]。

本研究提出的WPA算法是基于三维几何形态学分析编制自动化算法对标志点的对称性进行量化评估，并将量化评估结果作为标志点赋权参数，实现具有个性化特征权值分配的对称参考平面构建。通过计算WPA算法面部对称参考平面与专家定义对称参考平面间的夹角，显示30例下颌偏斜面部畸形患者的平均角度误差均小于2°，与传统无加权PA算法相比，WPA算法构建的结果与临床专家定义的“真值平面”更为接近(图3)，且对于严重不对称畸形患者(下颌偏斜程度大于12 mm组)的优势更为明显。Wu等[28]的研究表明，两平面角度差异大于6° 时易被人感知，WPA算法与真值平面的角度误差小于2°，表明WPA算法构建的对称参考平面的准确性几乎等同于专家参考平面。

本研究WPA算法的实现完全基于开放的python软件平台实现，从人脸标志点信息读取到对称参考平面的构建，不依赖于任何第三方软件，实现了完全自动化，与基于商业软件平台ICP算法的专家组算法相比，WPA算法具有更好的二次开发性，有利于数字化诊疗的基层推广。

3.2 EDMA法用于面部标志点不对称性量化评估具有可行性

本研究建立的WPA算法是在本体-镜像关联法的框架下，基于三维几何形态学分析中的EDMA法对三维面部重要解剖标志点的不对称程度进行量化评估，并以赋予标志点相应权重的形式体现其为对称参考平面构建的贡献差异，该赋权算法是本研究的一个创新点。以往文献报道，EDMA法已被应用于颅面形态测量[29-31]、性别二态性研究[32-33]、口腔正畸学牙弓形态学分析[34-36]等领域，其原理为通过个体几何形态上标志点间的欧氏距离组成的矩阵反映个体的形状和大小，并通过个体间对应矩阵元素的比值反映个体间的形态差异。聂琼等[37]进行错牙合畸形牙弓形态分析时，通过删除标志点找出对形态差异贡献相对较大的标志点。本研究借鉴了EDMA法个体间差异分析的原理，将人脸左侧、右侧标志点集分别与中线标志点集组合构建为左脸形态点集和右脸形态点集，将左脸和右脸作为互为镜像的个体，借助EDMA法进行左、右脸形态差异的量化分析，从而获得标志点不对称性的量化分析结果。通过删除某一标志点计算其余标志点的T(-i)值，分析该点对左、右脸三维形态不一致程度上的贡献，进而通过自动化算法对每个标志点的T(-i)值进行逐一计算和排序，不对称性排序最低者权值为1，最高者权值为0，左右对应标志点权值相同，权值根据排序呈线性单调递减。30例面下偏斜患者三维面部数据分析显示，鼻根点、鼻尖点、内眦点、外眦点等面上及面中部标志点的算法赋权总体较高，颏前点、颏下点、下颌角点等面下部标志点的算法赋权总体偏低，该算法结果较符合口腔临床针对面下偏斜病例的常规诊断经验，EDMA法用于面部标志点不对称性量化评估具有较好的临床适合性。

本研究基于EDMA值的WPA算法构建三维面部对称参考平面，突出优点是可以在不降低算法自动化程度的前提下(无需人为干预)，尽可能模拟面部对称参考平面构建过程中专家经验对解剖标志点参考价值权重的表达，初步验证了WPA算法对下颌偏斜面部不对称畸形患者具有较好的适合性，后续还需进一步扩大样本类型及提高算法的自动化程度(如三维面部标志点自动标定)，对不同面部分区畸形和不同程度面部畸形的适合性进行进一步的统计学和测量学分析，为临床应用提供指导。