基于条件生成对抗网络的缺失牙体功能性咬合面形态设计方法

2021-02-22田素坤袁福来孙玉春程筱胜

中国机械工程 2021年3期

田素坤戴宁袁福来孙玉春俞青程筱胜

1.南京航空航天大学机电学院，南京，2100162.北京大学口腔医学院口腔修复科，北京，100081 3.南京大学医学院附属口腔医院修复科，南京，210008

0 引言

龋齿是一种引起牙列缺损、丢失以及严重影响人类口颌健康的常见慢性口腔疾病，成年人(>20岁)患病率约为87%[1]。人类自身特有的咬合状态是牙齿存在和执行功能的基础，并且正确的咬合关系决定牙齿之间的接触形式以及影响牙齿磨耗形态特征的分布。牙齿磨耗是伴随人一生的生理变化过程，能真实反映下颌运动的生理状态，并且与口颌系统结构及功能的健康协调密切相关[2-3]。生理性磨损使上下牙颌面接触广泛，有助于消除早期咬合接触点，减少侧向力，重建咬合平衡。因此，恢复正确的咬合接触关系不仅能满足患者的生理需要，而且能维持修复体的长期稳定性和口颌健康。同时，建立合适的咬合接触关系也是评价修复治疗成功的标准。

随着计算机图形学的发展，牙科CAD/CAM修复技术凭借其高效性和智能化特征，逐步取代了传统口腔修复体的手工设计与制造[4-5]。目前大多数商业牙科CAD/CAM修复系统，如CERCE3D(德国)、3 Shape(丹麦)均以标准冠数据库作为修复软件中的重要组成部分[6-8]，借助与当前牙位一致的标准牙模型的形态信息，并采用一系列变形操作算法使其形态尽可能与缺失牙残留区域相互吻合[4,9]。尽管牙科CAD/CAM修复系统可以高效、准确地重建出缺失牙的咬合面形态，但是在设计过程中需要大量的手动交互操作，并且修复后的咬合面质量完全取决于医生操作的熟练程度及其主观设计思想，因此仍无法设计出最适合患者的个性化全冠修复体。

重建出正确的咬合接触关系以及具有天然牙齿的形态特征是牙齿修复成功的关键。因此，STEINBRECHER等[10]采用迭代Laplacian曲面编辑技术对标准牙表面进行变形操作，实现缺失牙齿咬合面的重建。在此基础上，JIANG等[11]为不同变形区域分配不同的权重因子来调整牙齿表面形态，但是牙齿网格模型质量不均匀导致迭代变形后的咬合面出现细节失真。依据同名牙齿间具有形态相似的咬合面，张长东等[12]采用基于颌位空间约束的统计方法构建缺损牙齿的平均牙冠模型，并通过分析缺损牙体与平均牙冠模型的特征差异重建出缺失区域，但该方法依赖于构建大样本同名牙齿数据库，并且需要手动选取大量特征点用于建立样本牙齿间点与点的对应关系。范然等[9]采用基于特征线约束的移动最小二乘变形算法实现标准牙与预备体的匹配。郑淑贤等[4]采用基于径向基函数的三维变形设计方法重建缺失嵌体的咬合面形态，该方法通过构建预备体与标准牙之间的特征映射关系实现缺损牙齿表面重建，但义齿表面形态完全取决于所选特征点的数量，并且特征点数量越多,设计过程越复杂。ZHANG等[13]提出一种双因子约束变形的牙齿嵌体建模框架，但该框架在嵌体设计过程中需要更多的用户交互。以上方法仅考虑在静止咬合状态下缺损牙齿咬合面的形态设计。针对这一问题，TIAN等[14]构建了一种基于相对轨迹驱动的虚拟动态调颌系统，重建出具有咀嚼功能的缺失牙齿咬合面形态，但是该系统的调颌效率较低。

以上方法主要通过手动控制变形点来调整标准冠的表面形态，虽然可以获得理想的咬合面形态，但是需要依靠医生熟练的操作技能以及鲁棒的变形算法，并且难以打破人为的主观设计。鉴于此，随着深度学习的快速发展，以卷积神经网络(convolutional neural network, CNN)为代表的典型网络结构具有较强的自学习能力和鲁棒性，并广泛应用于口腔医学各个领域[15-16]。LEE等[17]开发出一种基于GoogLeNet网络的计算机辅助检测系统，用于牙周炎和龋齿的诊断和预测。MIKI等[18]利用AlexNet网络从CBCT图像中成功标记出牙齿类别。XU等[19]采用提取的20×30图像作为2D CNN网络的输入对牙颌模型进行分割，并利用模糊聚类算法对分割边界进行优化, 有效地提高了牙颌模型的分割准确率。TIAN等[20]尝试使用基于3D CNN的多级层次学习网络框架从三维牙颌模型中提取高维特征，实现对牙齿类型的分割及分类，并且得到比较理想的结果。以静态咬合干涉作为约束条件，袁福来等[21]提出基于生成模型的咬合面形态设计方法，均未考虑牙齿磨耗的影响，以至于生成的咬合面不具备正常咀嚼功能。

本文以缺失牙齿咬合面形态设计作为研究对象，提出一种基于条件生成对抗网络(conditional generative adversarial network, CGAN)[22]的功能性咬合面形态生成网络。通过计算预备体、条件模型及其磨耗面的深度信息构建咬合面形态生成网络的样本数据集。利用颌位空间位置关系以及牙齿磨耗作为网络模型的约束条件，生成同时满足正确咬合接触关系和正常咀嚼功能的功能性咬合面。

1 全冠修复体咬合面形态设计方法

本文提出一种基于颌位空间约束以及牙齿磨耗引导的功能性咬合面形态自动设计方法，首先采用启发式搜索策略和正交投影方法提取牙齿磨耗面并生成对应样本数据集的深度图像；然后联合使用颌位空间位置关系和牙齿磨耗面训练两阶段式生成网络模型合成功能性咬合面。最后采用区域生长法和蒙皮法设计出具有功能性咬合面的全冠修复体。具体设计流程如图1所示。

图1 全冠修复体设计技术路线Fig.1 Overview of full crown restoration design workflow

1.1 数据预处理

1.1.1牙齿磨耗面提取

牙列中牙齿的形态、位置以及咬合运动的特异性，都直接影响牙齿接触区域上磨耗的分布形式。磨耗面近似为一个表面光滑且边界清晰的倾斜平面[2-3]，通常在磨牙上主要分布于咬合面、切牙边缘及邻接面上，如图2所示，图中颜色表示不同咬合运动方向形成的磨耗面形状[3]。磨耗面边缘轮廓属于特征线的范畴，并且描述了咬合面的功能特征区。

图2 磨牙牙齿磨耗分布[3]Fig.2 Tooth wear facet distribution on molar[3]

生物特征线提取是牙科CAD系统中的关键环节，启发式搜索技术因其高效性和复杂度低而广泛应用于特征线提取中。实验采用基于启发式搜索策略的交互半自动提取方法获得磨耗面的光滑边缘轮廓。启发式搜索算法描述从网格模型的一个顶点到其相邻顶点的有效可遍历路径过程[13]。从起始点到路径终点的局部代价和最小，则为搜索的最低成本路径[4]。图3a为启发式搜索示意图。

(a)启发式搜索示意图

(b)边缘轮廓优化方法图3 牙齿磨耗提取方法Fig.3 Tooth wear facet extraction

在磨耗面边缘轮廓提取中，最优路径被定义为从给定的起始点VS到终点VE的最小累积成本路径。局部启发式代价函数f(v)定义如下:

f(v)=g(v)+h(v)

(1)

g(v)=fdir+fD1+fCh(v)=fD2

式中，g(v)为从起始点到当前搜索点的实际代价函数；h(v)为当前搜索点到目标点的启发代价函数；fdir为方向代价函数，可以维持特征线具有较好的平滑效果；fD1和fD2为距离代价函数，用于限制路径搜索的终点；fC为特征代价函数，确保搜索到的目标点最大程度描述磨耗面边缘轮廓的显著特征。

直接采用启发式搜索方法提取的磨耗面边缘轮廓线呈锯齿状，因此采用一种简单高效的路径优化算法来平滑磨耗面的边缘轮廓，如图3b所示。假定Vi-1、Vi+1为磨耗面边缘轮廓上两点，Vi为这两点的中间点。优化算法具体步骤如下。

(1)计算向量Vi-1Vi与ViVi+1的夹角φ，当φ>170°，光顺该段路径；否则，执行步骤(2)。

(2)构造Vi一环邻域三角面片集合T，计算Vi-1与Vi+1连线的中点Vm，并将距离Vm最近的网格顶点标记为Vn，通过Vi-1、Vi+1和Vn构建平面。

(3)计算平面与集合T中三角面片的交点，用计算出的交点(V′i和V″i)替换Vi(图3b)。

(4)重复执行步骤(1)～步骤(3)，直到完成搜索路径的优化。

1.1.2深度图像构建

特征提取技术是基于内容的三维模型检索的核心，其结果的优劣直接影响生成网络的性能。考虑三维模型高自由度及不同模型间坐标系不一致会对检索结果产生较大影响，本文采用包围盒技术和虚拟旋转球技术对三维模型进行规范化预处理，用于提高检索的效率和准确性，并且使处于不同坐标系的相似模型具有相同的相似度。

深度图像是带有距离信息的视觉图像[23]，为使投影得到的深度图像最大程度地保留牙冠的空间形状信息，采用自适应视距的正交投影方法，即在正视方向上，将投影板自适应地调整到牙冠的包围盒处，并将其划分为256×256的网格面板，如图4所示，使投影得到的深度图像细节特征更丰富，区分度更高。

通过计算每个网格中心点到牙冠表面的最短垂直距离d，构建网格中心点P(xm,yn,dmn)与对应像素点Pp(m,n)之间的转换关系：

Pp(m,n)=255max(0,h-‖P(xm,yn,dmn),Po‖)/h

(2)

式中，h=7为距离阈值；‖·‖表示从网格中心点P(xm,yn,dmn)到咬合面的最短垂直距离。

图4 深度图像生成示意图Fig.4 Depth maps generation

1.2 咬合面生成网络模型

随着大规模训练数据集的出现，以及pix2pix作为一种典型的生成网络结构在图像到图像转换任务中取得的突破性进展[22]，本文在此基础上构建颌位空间位置约束以及牙齿磨耗引导的功能性咬合面生成网络模型，如图5所示。该网络由两个相似的生成模型和一个判别模型构成。直接以牙齿磨耗面为条件生成功能性咬合面较为困难，因此本文采用两阶段式(two-stage)方法来解决这一问题，并且每个阶段只关注一个方面(图1)。对第一阶段(Stage-Ⅰ)，使用pix2pix的变体来生成满足颌位空间位置关系的咬合面基本形状；而对第二阶段(Stage-Ⅱ)，通过在Stage-Ⅰ结果中添加牙齿磨耗面约束条件来增添咬合面的细节特征，进而生成具有功能特征的咬合面形态。

图5 咬合面生成网络模型Fig.5 Generative network model of occlusal surface

1.2.1条件生成对抗网络

生成对抗网络(generative adversarial network, GAN)[24]的出现为计算机视觉的应用提供了新方法和新框架。受博弈论中的二人零和博弈的启发，GAN由一个生成器G和一个鉴别器D组成，并且两者都受到对抗性学习机制的训练。GAN的目的是估计现有数据的潜在分布并从同一分布中生成新的数据样本。条件生成对抗网络(CGAN)是在GAN的基础上添加条件扩展而成的，其目标函数V(D,G)为

Ez～Pz(z)(log(1-D(G(z,c))))

(3)

式中，E(·)为计算期望值；c为条件变量；x为采样于真实数据分布Pdata(x)；z为采样于先验分布Pz(z)；D(x,c)为判别器鉴别输入数据来自目标数据分布的概率；G(z,c)为生成器生成的数据。

1.2.2Stage-Ⅰ颌位空间约束的咬合面生成

鉴于条件生成对抗网络具有可控的训练过程、稳定的生成模型以及条件可扩展性，本节以预备体作为输入数据x1，无磨耗牙冠作为目标数据z1，生成满足颌位空间位置关系的咬合面基本形态。对抗损失促使生成图像的分布接近真实图像，但是，该损失不直接惩罚实例到实例的映射[21]。因此，采用L1距离函数作为Stage-Ⅰ的生成损失：

LL1(G)=Ex,c,z(|z-G(x,c)|)

(4)

式中，x为输入数据；z为目标数据；G(x,c)为生成器生成的数据；|·|表示绝对偏差。

牙齿修复治疗中，只关注咬合面的解剖学特征，则无法完全恢复缺失牙齿的功能特征。因此，必须考虑咀嚼运动过程中上下颌牙齿间的接触情况，即理想的修复体咬合面与其对颌牙齿间存在一定的干涉强度[14]。本节采用对颌牙齿c1及与对颌牙间隙d作为颌位空间约束条件，引导生成器G1合成具有正确咬合关系及咬合接触分布的初始咬合面形态。其目标函数为

LcGAN(G,D)=Ex1,c1,d,z1(logD(x1,c1,d,z1))+
Ex1,c1,d,c(log(1-D(G(x1,c1,d,c))))

(5)

式中，x1为预备体数据；z1为含无磨耗的目标冠数据。

将L1损失重新表示为

LL1(G)=Ex1,c1,d,z1(|z1-G(x1,c1,d)|)

(6)

结合L1损失约束后，生成器G1和判别器D1的损失函数分别表示为

(7)

(8)

式中，λL1为平衡权重。

1.2.3Stage-Ⅱ磨耗约束的功能性咬合面生成

基于Stage-Ⅰ GAN生成的具有与目标牙冠相似姿态的初始咬合面，Stage-Ⅱ中的生成器G2旨在生成咬合面特征的差异映射，使初始牙冠更接近目标牙冠，并且以Stage-Ⅰ的结果和咬合磨耗面作为输入,用于丰富咬合面的功能特征。以磨耗面c2作为功能特征约束条件，引导生成器G2合成具有颌位空间位置关系的功能性咬合面。为计算具有目标冠z2的生成器G(x2)，同样采用L1距离函数作为Stage-Ⅱ的生成损失，其目标函数表示为

LcGAN(G,D)=Ex2,c2,z2(logD(x2,c2,z2))+
Ex2,c2,c(log(1-D(G(x2,c2,c))))

(9)

式中，x2为无磨耗牙冠数据；z2为含有磨耗的目标冠数据。

最终生成器G1和判别器D1的损失函数与式(7)、式(8)相同。

1.3 三维牙颌网格模型重建

散乱点云数据的网格曲面重建在CAD/CAM中起着重要的作用。本文采用Delaunay三角剖分优化准则，实现基于区域增长法[13]的牙颌网格模型重建，并取得了良好效果。区域增长法是从种子三角形出发，不断将其邻域中的最优标记点纳入新生成的三角网格中。定义点云数据集合P={Pi|i= 1,2,…，n}，该方法的基本步骤如下。

(1)选取点云数据集合中的点P1作为种子点。

(2)计算邻域内点到P1的欧氏距离，并以距离最小的外部点作为下一个点P2，连接这两点形成边界边P1P2，并将其添加到边界边链表Edge_List中。

(3)根据三角剖分优化准则，寻找能与P1P2形成最佳三角形的最佳点P3。

(4)连接边P1P3和P2P3，判断其是否包含在链表Edge_List中，否则将其作为新边界边插入到该链表中。

(5)选取链表Edge_List中的队首元素，重复步骤(3)和步骤(4)，直到链表中的边界边全部配对为止。

区域增长过程中，边结构存储在链表Edge_List中，该链表存储三角形面片的三个顶点，并且首尾相连形成封闭空间。三个顶点按照“右手准则”存储，便于计算三角形面的法向量。

1.4 全冠修复体设计

图6 上颌第一磨牙解剖特征Fig.6 Detailed features of the first right mandibular molar

设计前，首先将预备体偏置一定距离，用于模拟黏结层的使用(图7)。其次通过使用生成的咬合面(SecL1)和黏结层(SecL2)的边界曲线作为基线来设计连接体。蒙皮剖分作为CAD系统中常用的曲面定义方法[25]，实验为增加连接体形状的可控性，采用基于B样条曲面插值的蒙皮剖分方法来自动设计中间连接体。

图7 全冠设计结构示意图Fig.7 Sketch map of full crown design

具体设计步骤如下。

(1)提取边界曲线SecL2上的所有点，即SQ= {qi|i= 1, 2，…，n}。

(2)以SQ为参考点序列，采用平面相交法计算边界曲线SecL1上对应的匹配点序列SP，即SP= {pi|i= 1, 2,…，n}。平面相交法原理如图8a所示，通过两个中心点(c1，c2)和参考点qi构成的平面计算SecL1上的交点。由于边界曲线封闭，所以相交于两点(pi,p′i)，如果piqic2>p′iqic2，则pi与qi为对应的匹配点，否则p′i为计算的交点。

(3)计算两组参考点序列的中点序列SK：SK={ki|ki=(pi+qi) /2,i=1,2，…，n}。

(4)以(pi,ki,qi)作为一组控制点，计算B样条曲线作为脊线，并对其均匀离散化(图8b)。

(5)依次连接相邻两脊线之间的交点(图8b)，得到连接体的网格曲面。

(a)平面相交法原理图 (b)网格划分图8 连接体设计示意图Fig.8 Sketch map of the connector design

2 实验结果分析

2.1 实验数据及网络参数设置

本项目与北京大学口腔医院和南京市口腔医院合作，通过3Shape扫描设备为实验收集患者受损及修复后的三维数字牙颌数据，并采用自主研发的功能模块交互提取牙齿磨耗面。考虑下颌第一磨牙患龋率最高[1-2,4]，因此选择缺失36号或46号牙的患者作为研究对象。依据不同个体同名牙齿间具有形态特征相似的咬合面[4,9]，实验选取500个样本用于训练，100个用于测试。对于每个样本，都对应5组数据：预备体、牙齿磨耗面、对颌牙齿、无磨耗冠(牙科专家设计的未经过打磨的牙冠)、含磨耗冠(咬合运动过程中反复测试打磨的牙冠)。

为训练两阶段式生成模型，本文使用Adam方法对网络进行优化[26]，同时设置参数：λL1= 100，动量参数β1= 0.5，动量参数β2= 0.999，ReLU斜率为0.2，初始学习率为0.0002，批次为1，迭代次数设置为50 000。

2.2 结果分析

为验证提出的咬合面形态生成网络的有效性和正确性，实验采用基于TensorFlow的深度学习框架[26]在硬件配置为Intel(R) Xeon(R) Platinum 8168CPU(2.70 GHz)和GeForce GTX 1080Ti GPU 的计算机上进行网络训练。

2.2.1训练收敛性及测试结果

训练过程中，通过训练损失调整网络参数，以监视网络的收敛速度。图9显示咬合面生成网络模型的训练损失与迭代次数的关系。随着训练的进行，损失值持续下降，表明该网络在开始时成功克服了梯度膨胀/消失的问题，实现了损失值的稳定下降，同时说明所提出的网络模型可实现有效训练。

图9 训练损失曲线Fig.9 Training loss curves

如上所述，验证了提出的功能性咬合面形态生成网络模型的3个子任务，测试结果如图10所示。对于任务一，直接训练Stage-Ⅰ GAN生成256×256的牙齿磨耗图像，用以验证牙齿磨耗作为指导条件的有效性。由图10a可以看出，在预备体数据上可以生成与目标磨耗分布非常接近的结果，表明牙齿磨耗作为约束条件能够起到引导作用。如图10b所示，可以发现Stage-Ⅰ GAN能够生成比较合理的缺失牙咬合面图像，并且可以捕获牙冠的整体结构，如牙尖和牙沟(图6)。但是，Stage-Ⅰ的结果比较粗糙，并且丢失许多细节特征。如图10c所示，Stage-Ⅱ GAN生成了具有更加明显尖窝分布的咬合面图像，能更好地反映咬合面的功能特征。

(a)磨耗面 (b)Stage-Ⅰ(c)Stage-Ⅱ图10 实验测试结果Fig.10 Experimental test results

2.2.2生成图像质量评估

依据上述给出的可视化结果(图10b和图10c)，分别采用结构相似性(structure similarity index measure，SSIM)和特征相似性(feature similarity index measure，FSIM)定量评估咬合面生成的质量。SSIM用于评估两幅图像之间的整体结构相似性;FSIM从不同角度测量两幅图像上局部位置像素的差异。表1给出不同阶段生成咬合面图像与目标冠图像之间的评价结果。与Stage-Ⅰ GAN相比，Stage-Ⅱ GAN的SSIM和FSIM值相对增加了0.059和0.032，证明了提出的网络模型具有生成高质量咬合面形态图像的能力。

表1 生成图像评价指标

2.2.3咬合面质量评估

依据距离-灰度关系，采用网格重建算法得到三维牙齿网格模型，结果如图11a所示。为评估咬合面的质量，在相同邻牙约束的条件下，计算出每个阶段生成的咬合面与目标牙冠之间的偏差(图11b)，并且标记出主要功能区的偏差位置。从图中可以看出，主要功能区域中，Stage-Ⅰ结果的偏差相对较大(最大约为0.92 mm)。如图11a的第一列(Stage-Ⅰ)所示，在大多数情况下，Stage-Ⅰ GAN能够生成缺失牙齿的基本表面形态。以牙齿磨耗面和Stage-Ⅰ的结果作为条件，Stage-Ⅱ GAN可为Stage-Ⅰ的结果增加功能特征，并且标记区域偏差值(最大约为0.78 mm)均在允许的干涉阈值范围内[14]。更重要的是，该结果表明Stage-Ⅱ GAN模型能够生成真实的图像结果，并且包含比单一阶段模型更令人信服的细节特征。

(a)三维牙颌重建结果(b)偏差分析图11 三维牙颌模型质量评估Fig.11 Quality evaluation of 3D dental model

2.3 真实患者全冠修复体设计过程

为验证本文方法的临床适用性，图12给出了一名冠缺损患者(牙齿编号#46)的全冠设计完整过程。其中包括咬合状态下缺失牙齿的邻牙及对颌牙数据(图12a)。首先采用启发式搜索算法和正交投影法提取牙齿磨耗并计算样本数据的深度图。然后利用已训练好的Stage-Ⅰ GAN和Stage-Ⅱ GAN模型得到的最佳参数进行实验测试，得到具有患者正常咀嚼功能的缺失牙咬合面(图12b)。同时，采用网格重建算法将生成的缺失牙齿的咬合面图像重建为三维网格模型(图12c、图12d)。其次结合计算出的黏结层和中间连接体，设计出缺失牙具有功能性咬合面的全冠修复体(图12e)，并且与邻牙保持良好的接触关系(图12f)。最后从图12c可以看出，生成的牙冠外表面具有天然牙齿的个性化解剖特征(图6)，同时符合文献[27]中得到的基于特征线的磨牙设计结果，其中，Curve BC、Curve CF和Curve LC分别表示连接相邻牙齿颊侧尖点、牙齿中央窝和牙齿舌侧尖点形成的曲线。经牙科专家鉴定，该方法设计出的全冠修复体具有较高的临床适用性。

图12 全冠设计过程Fig.12 Full crown design process

3 方法对比

为进一步验证两阶段式生成模型的有效性，本文将其与CAD方法[10-11]、统计学方法[12]以及基于生成模型的方法[21]进行比较。结果如图13所示。

图13a为在咬合运动过程中反复测试和打磨的牙冠，即为目标冠，标记的主功能区为牙齿磨耗分布的主要区域。图13b为采用Laplacian曲面变形技术设计的磨牙冠，该方法首先从标准冠数据库中选取同名牙齿模型，然后通过定位和缩放操作获取缺失牙冠的初始近似位置形态，最后通过控制牙冠表面上的一组变形点实现咬合面形态调整，但是由于操作过程中缺乏对咬合面合理有效的控制，导致咬合面局部变形区域出现畸形或功能特性丢失。图13c为基于统计学方法计算的平均牙冠，该方法采用200个同名样本冠计算出一个具有生物遗传学形态的牙冠，但是该方法需要在每个样本牙冠表面手动选取50个特征点进行样本匹配，并且计算得到的牙冠细节特征比较少。图13d是使用生成模型得到的磨牙冠。该网络模型以静态咬合状态下对颌牙齿间的间隙距离为条件生成牙齿咬合面，并未考虑牙齿磨耗对全冠修复体的影响，因此该结果不能真实地反映咀嚼运动过程中咬合面的功能特征。通过对比图13a、图13d和图13e发现，本文提出的两阶段式方法可以在咬合面上生成更多合适的功能特征区域。