基于CBCT图像的颌骨及牙列有限元建模方法研究

王爽 孙江 于雁云

三维有限元法作为生物力学的一种重要研究手段，于1973 年被Thresher等[1]首次应用于口腔医学。随后几十年中有限元法在口腔生物力学研究中的应用越来越广泛。在正畸领域中，有限元法被广泛用于模拟正畸力作用下的牙齿移动[2]，分析移动过程中牙周膜的力学响应[3-4]，预测矫治器作用下的牙周膜应力[5-6]等问题中。在口腔修复领域，有限元法用于评估牙齿在动态咬合力作用下的力学响应[7]，分析不同种植方案或者种植体对应的牙齿及牙槽骨的应力，以便于确定生物力学特性最佳的修复材料[8]或者确定最佳的修复方案[9]等。

建立合理的有限元模型是进行有限元分析的重要基础。目前，有限元建模主要包括以下3 种途径，即基于CT图像[3-4,9-10]、基于Micro CT图像[7]与基于CBCT的建模方法[2,5-6]。CT采用扇形束扫描，空间分辨率相对较低，尤其对细微结构显示欠佳，容易产生辐射状伪影， 所以基于CT建立的有限元模型精度相对较低。Micro CT扫描范围有限且辐射剂量大，其扫描对象通常仅为离体标本或者小动物。CBCT采用锥形束X射线扫描，只需旋转360°即可获取重建所需的全部原始数据，具有辐射量低、各向同性高空间分辨率等优势，在临床中广泛应用。因此，CBCT相比于CT与Micro CT更适用于口腔临床的有限元建模。现有的文献中，缺少基于CBCT影像的牙齿、牙周膜与牙槽骨有限元建模方法的系统研究。本文提出一种实用的颌骨及牙列三维有限元建模方法，可广泛用于各种口腔生物力学有限元分析问题。

1 材料与方法

1.1 材料

选择1 名牙周健康，牙齿排列整齐，牙列完整的健康男性青年志愿者，采用CBCT(KaVo 3D eXam,美国)进行扫描，采集下颌骨及牙齿三维信息获得DICOM格式数据。

1.2 颌骨及牙列三维重建

将DICOM格式数据导入Mimics17(Materialise，比利时)软件中，通过设置合理的灰度阈值并三维重建，得到由三角面片组成的颌骨及牙列三维模型。图 1A为灰度阈值取700～4 000三维重建得到的颌骨及牙列整体模型，该模型中包含了牙齿、牙槽骨和牙周膜，但是三者连为一体，无法得到牙齿与牙周膜的边界、牙周膜与牙槽骨的边界。由于牙齿、牙槽骨和牙周膜具有不同的物理属性，有限元分析中需要将三者分离并赋予对应的物理属性。所以，上述方法建立得到的颌骨及牙列整体模型，可以用于临床中牙齿及牙槽骨的形态观察、参数测量等，但不能直接用于建立有限元模型。不同牙位的牙骨质羟基磷灰石含量有一定差异，为了识别出所有牙齿，需要取较大的灰度阈值范围。同样牙槽骨的羟基磷灰石含量分布也不均匀，需要设置较大的灰度阈值范围才能识别出所有牙槽骨。同时，牙骨质和牙槽骨皮质羟基磷灰石含量接近，CBCT图像上反差相对小，因此二者的灰度阈值范围存在较大重叠，无法单纯通过设置灰度阈值范围将二者分离。如图 1B所示为灰度阈值取1 300～4 000的三维重建结果，模型中磨牙的牙根由于灰度阈值下限较高而缺损，但是模型中仍然包含大量的牙槽骨。

1.3 牙列三维建模

牙齿三维建模的核心问题在于从牙槽骨中识别出牙根，识别牙根的关键在于灰度阈值的设置。令某颗牙齿牙骨质灰度阈值范围为[LT,UT]，对应牙位牙槽窝处牙槽骨灰度阈值范围为[LA,UA]，通常情况下，LT>UA。识别牙齿时灰度阈值下限应取LT与UA中间的值，以识别出牙齿的外表面。LT与UA没有固定标准，其取值因人而异，且同一人不同牙位处的值也存在较大差异。因此建模中需要根据经验，通过数次尝试确定各类牙齿LT与UA中间的值，以确定识别牙齿时的灰度阈值下限。必要时需要针对不同牙位的牙齿设置不同的灰度阈值进行三维建模。

A: 灰度取值700～4 000； B: 灰度取值1 300～4 000图 1 灰度阈值范围取700～4 000与1 300～4 000对应的三维模型A: Gray threshold value of 700-4 000; B: Gray threshold value of 1 300-4 000Fig 1 Three dimensional models for gray threshold value of 700-4 000 and 1 300-4 000 respectively

将三维重建得到的三维模型保存为.stl格式文件。应用Geomagic Warp 2015软件(Raindrop，美国)打开并编辑.stl文件。牙齿三角面片与牙槽骨三角面片存在大量的连接关系，不能简单通过区域扩张等方式将二者分离。牙齿与牙槽骨交线处，二者三角面片的法向存在明显的差异，因此可以通过设置相邻单元法向角差值阈值，识别出牙齿三角面片，并删除其余三角面片。然后修补缺损的牙齿面片，删除尖角、凹陷等表面特征，得到牙齿整体模型。

经过上述方式建立的牙齿整体模型中相邻牙齿在牙冠处相连。虽然可认为各牙齿的物理属性相同，但是有限元模型中必须将其分离才能在有限元分析中得到正确的力学响应。单颗牙齿中相邻单元的法向量差异较小，而相邻牙齿在交线处法向量有明显差异。通过设置相邻单元法向角差值阈值，将相邻牙齿分离，然后借助于填充孔等方法修复分离牙齿导致的孔，得到完整的单颗牙齿模型。将所有牙齿分离并修复后得到牙列模型。

1.4 牙槽骨与牙周膜三维建模

牙槽骨内部骨松质存在大量的骨小梁，骨表面解剖结构复杂，凹凸不平致使骨密度不均，在三维重建过程中，难免出现大量不规则的尖角形状。含有不规则形状的三维模型在有限元网格划分时可能产生低质量单元，甚至导致网格划分失败。因此，用于有限元分析的三维模型应当在不影响颌骨及牙列力学特性的前提下进行必要的简化，去掉尖角、凹陷等特征。先通过设置相邻单元法向角差值阈值识别牙槽骨外表面，删除牙槽骨外表面的多余特征，通过填充方式将牙槽窝处的孔填平，建立不体现牙槽窝、骨松质的牙槽骨整体模型。

识别牙槽骨松质与骨皮质分界面，采用删除特征、填充孔等工具，在不对模型精度造成过大影响的前提下，将分界面修复成一张相对光顺的三角面片曲面，得到牙槽骨松质模型。利用三角面片模型的布尔减运算，从牙槽骨整体模型中扣除牙槽骨松质模型，得到不体现牙槽窝的牙槽骨皮质模型。牙槽骨松质模型与牙槽骨皮质模型构成不体现牙槽窝的牙槽骨整体模型。

受限于CBCT的分辨率，当前硬件与软件条件下不能从CBCT的DICOM格式数据中直接识别出牙周膜。因此本研究中假定牙周膜厚度均匀(取0.25 mm)。首先，将牙列模型沿各顶点的法向向外侧偏移0.25 mm，得到牙齿向外偏移模型。利用三角面片模型的布尔干涉运算，计算牙槽骨整体模型与牙齿向外偏移模型的交集，得到牙槽窝模型(其中包含牙周膜)。通过布尔减运算从牙槽窝模型中扣除牙列模型，得到牙周膜模型。利用布尔减运算从牙槽骨整体模型中扣除牙齿向外偏移模型，得到牙槽骨模型。

1.5 三维建模流程及其优点

前面所述的牙齿、牙周膜与牙槽骨的三维建模流程可由图 2更清晰表达。该建模方法可以完成上颌、下颌的牙齿、牙周膜与牙槽骨模型的建立。

图 2 牙齿、牙周膜与牙槽骨的建模流程图Fig 2 Modeling procedure chart for tooth, periodontal ligament and alveolar bone

上述流程中涉及到的布尔减运算、布尔干涉运算是指三角面片的布尔运算，其运算结果仍然为三角面片模型，其中不涉及任何三维实体相关运算。由于牙周膜、牙槽窝模型均通过布尔运算得到，所以牙齿与牙周膜、牙周膜与牙槽骨在分界面处的三角面片划分相同，可有效保证有限元模型中牙齿与牙周膜、牙周膜与牙槽骨连接关系的正确性。采用上述方法，建立图 1对应牙齿与牙周膜的模型如图 3A所示，牙槽骨模型如图 3B所示，下颌骨及牙列整体模型如图 3C所示。

A:牙齿及牙周膜； B：牙槽骨； C：整体图 3 牙齿与牙周膜、牙槽骨及整体的三维模型A: Periodontal ligament; B: Alveolar bone; C: Entire modelFig 3 Three dimensional models for tooth, periodontal ligament, alveolar bone and the entire model

2 结 果

应用图 2所示流程建立牙齿、牙周膜与牙槽骨的三角面片模型。有限元分析需要基于牙齿、牙周膜与牙槽骨的体网格完成，因此需要对牙齿、牙周膜及牙槽骨三角面片模型进行网格划分，得到体网格，建立有限元模型。

2.1 三角面片网格优化

以三角面片的最小内角为25°作为约束对三角面片模型重新网格划分，以提高三角面片模型的质量。由于牙周膜厚度取为0.25 mm，三角面片模型中网格的最大边长设置为0.5 mm，以保证体网格质量。重划网格时，应注意使牙槽骨与牙周膜接触面、牙周膜与牙齿接触面的网格划分一致，即通过共节点的方式建立牙槽骨与牙周膜、牙周膜与牙齿之间的正确连接关系。

2.2 体网格划分原则

理论上，六面体网格相比于四面体网格具有更高的计算精度。但是由于三维模型采用三角面片表达，无法在三角面片的基础上划分六面体单元，所以需要将其划分为四面体单元。选择10节点二次四面体网格，可实现通过较少的单元获取较高的计算精度。

对于一个成年人的牙槽骨与牙齿模型，如果采用等边长网格划分方案，即使边长取三角面片边长0.5 mm，体网格的规模将达到数百万量级，有限元模型的自由度总数可超过千万。上述问题对于普通的台式计算机求解十分困难。有限元分析的重点关注区域是牙周膜及其相连区域的牙齿与牙槽骨，加大非重点关注区域的网格尺寸对有限元计算结果不会造成影响。因此，在划分体网格时，不限制单元的最小边长，将四面体单元最小内角不小于25°作为约束条件，这样可保证牙周膜、牙周膜附近的牙齿与牙槽骨划分为细网格，其他区域划分为粗网格，在确保计算精度的同时有效降低模型的规模。同时，由于牙齿与牙周膜接触面、牙周膜与牙槽骨接触面的三角面片顶点网格划分一致，所以基于三角面片划分得到的四面体单元在上述两个接触面处满足网格一致性。采用上述原则，建立图 3所示牙齿、牙周膜有限元模型如图 4A、 4B所示，整体模型如图 4C所示。模型的节点总数为299286，牙齿单元数为105805，牙周膜单元数为122427，牙槽骨单元数为577529。

A: 牙齿； B： 牙周膜； C： 整体图 4 牙齿、牙周膜及整体有限元模型A: Periodontal ligament; B: Alveolar bone; C: Entire modelFig 4 Finite element model for tooth, periodontal ligament and mandible with full dentition

3 讨 论

本文提出一种基于CBCT影像资料的颌骨、牙周膜及牙列的有限元建模方法，该法可以针对临床患者建立个性化的高质量全口牙齿、牙周膜及牙槽骨有限元模型，以便于针对临床具体问题进行生物力学分析，为精准医疗提供理论基础。该法适用于口腔生物力学中的正畸、修复及牙周治疗等临床问题，也可用于上述问题的学术与科学研究，具有良好的通用性。

本文提出的颌骨及牙列有限元建模方法中，有限元网格划分前的所有三维操作均基于三角面片模型完成。现有文献中建立牙齿及牙槽骨模型通常基于边界表达法(Boundary Representation，BRep)的三维实体模型完成，即在构造牙齿、牙槽骨的表面模型之后，将其拟合为NURBS曲面，依据NURBS曲面构造BRep实体模型。对BRep实体模型进行网格划分，得到有限元模型[5-6，10]。传统方法存在明显的缺点。首先，将三角面片模型拟合为NURBS曲面模型时，拟合效率受限于曲面算法与曲面边界的选取，如果曲面较复杂或者轮廓线选取不当，会导致曲面拟合失败。其次，在对BRep实体进行布尔干涉、布尔减操作时，容易产生距离较近的点或者线。这些近距离的点、线对于三维显示没有影响，但是会给有限元网格质量造成不利影响，甚至导致网格划分失败。所以，传统的方法主要适用于单牙模型的建立[3，8，10]，或者相对简单、精度要求不高的牙列及牙槽骨建模[4-5，9]。本文提出的方法中三维模型采用三角面片模型，因为避免了NURBS曲面构建、NURBS实体缝合和BRep实体网格划分等问题，因此效率更高、容错性更强、三维模型更为精确且适用范围更广。