不规则复杂牙冠模型的建立及优化分析

闫胜昝， 章 健， 苏本跃

(安庆师范学院计算机与信息学院，安徽 安庆 246011)

不规则复杂牙冠模型的建立及优化分析

闫胜昝， 章 健， 苏本跃

(安庆师范学院计算机与信息学院，安徽 安庆 246011)

从牙齿扫描的点云数据出发，提出优化的快速行进距离场构建算法用于实现牙冠模型的高精度快速生成。借助计算机辅助工程中的应力有限元分析方法，获得网格单元上的受力情况，建立有限元实体网格单元与三角面片网格模型面片间应力关系模型。利用仿真引导设计的方法实现常规临床治疗手工加工过程的数字化和虚拟化，从而有效加速建模过程，并将模型修整部分提前到设计阶段。通过创建并优化后的模型可简化治疗程序、缩短治疗周期、减轻病人痛苦。

不规则网格曲面；点云数据；有限元分析；距离场；力矢量分解

口腔修复学作为口腔医学的一个重要组成部分，是医学与现代科学技术相结合的产物。近年来，随着计算机技术的不断发展，特别是计算机图形图像处理技术、虚拟现实技术和现代医学的飞速发展与交叉融合、相互渗透，利用计算机进行口腔修复并设计出实用的医学仿真系统成为计算机应用研究的一大热点，其研究基础是牙齿建模研究。CAD源自工业领域，广泛应用于制造业。工业上使用的CAD系统，设计加工对象通常是标准的几何体，具有一定的数学表达，而牙齿模型形态各异，多为复杂的不规则曲面，因此在实现个性化设计、满足美观和精度等多方面的要求上有一定的难度。

近年来，制作义齿的相关技术取得了一定进展，Williams等[1-3]采用具有力回馈功能的人机交互设备，进行了可摘义齿的CAD与金属支架快速成形，并在患者口内初步验证了适合性。国内学者采用两步法进行全口义齿的CAD软件开发以及全口义齿修复体的快速成形辅助制作，初步完成了全口义齿辅助制作系统的基础理论和基本方法的探索工作，实现了CAD技术在全口义齿修复领域的突破[4-6]。随着研究和应用的不断深入，在各种CAD系统中，相关数据的收集、分析、归纳和整理在逐步深入，与设计相关的口腔修复数据库也日益完善，特别是在嵌体、贴面、固定桥、个性化种植体的形态和功能设计上都已具有较高的水准和精度。然而，CAD阶段的问题也并未完全解决，口腔修复CAD系统的研究方向主要集中于传统修复领域的拓展，制作完整修复体时必须考虑咬合、邻面接触关系，需要模拟个性化下颌咬合运动关系进行检查、调整。而在设计制作牙齿全冠模型时，需要一种合理有效的等距算法，对分割后的牙冠模型进行等距放大，然后利用得到的等距面和原牙冠模型相结合得到牙齿全冠模型。在设计时，应避免饰瓷层过厚、过薄以及由于基底冠厚度不均匀导致的材料应力，采用适当的曲面等距算法，可以将饰瓷崩裂的风险降至最低。传统的对网格模型的等距方法主要包括基于面的等距[7]和基于点的等距[8]。其中，基于面的等距方法是先将网格模型中的所有面片沿各自法向量方向移动一个定值，然后，通过裁剪和延伸等操作，使其构成闭合网格，该方法在操作过程中需要处理大量的面片相交问题。而基于点的等距则是将网格中的各个控制点沿其法向量方向移动一个定值，然后，根据原模型的拓扑关系来生成新模型。这两种方法的优点是运算速度快，在等距距离较小时比较实用。然而，当等距距离较大时，自交与断裂现象的出现将不可避免。此外，使用隐式曲面的等距[9]也是一种思路，其主要原理是通过选取原模型中部分控制点，由这部分控制点及其法向量信息构造出一个基于径向基函数(radial basis function，RBF)的隐式曲面，适当调整隐式曲面中的参数便可得到等距曲面。这种算法可以在等距距离较小的情况下不产生自交，缺点是生成的等距模型精度不够，主要原因在于构造隐式曲面时，需要大量的控制点作为RBF的插值点，当插值点增多时，运算的复杂度将会急剧增高，而如果减少插值点数量则必然会影响到生成的等距曲面的精度。对于控制点数目较大的牙齿模型来说，使用隐式曲面的等距在精度问题上不能满意。为了避免出现自交现象，同时在保证精度的前提条件下提高算法的时效性，Lorensen等[10-11]采用移动立方体算法(marching cubes，MC)或双重轮廓线算法(dual contouring，DC)来处理等距问题。

本文利用基于 DC算法构造的优化等距算法处理等距问题，采用优化的快速行进算法(fastMarching)构建距离场，每次循环均计算采样点到最近三角面片的距离值，而非用从前次循环递推出的间接距离表示，在保证牙体模型精度的前提下，获得较好的提速效果；引入CAE分析方法计算牙冠模型的应力/应变；利用有限元分析的应力/应变结果指导牙冠模型调整优化，提高模型一次制作精度，减少对其进行手工精加工和调磨所需的时间，减轻患者痛苦。

1 模型重建

1.1 优化的fastMarching距离场构建算法

Sethian[12]提出的fastMarching算法，是一种用于处理运动界面演化问题的算法，现已广泛应用于图像处理，曲面重建的领域。fastMarching算法主要分3步：

(1) 将空间中所有采样点分为3类：① alive：已经被访问过的点；② active：正在或即将被访问的点；③ faraway：远离的点，是否需要访问未知。

算法先将原始数据作为种子，在二维平面中是一条曲线，在三维空间中就是一个曲面。对于牙齿模型来说，可以将原模型的点云数据或三角网格的控制点作为种子。构建距离场时，由种子出发，种子标记为active，均匀栅格中的所有采样点标记为faraway。

(2) 访问所有标记为active的采样点，将访问到的采样点的标记改为 alive。然后，访问此采样点 1邻域范围内的其他采样点，对于均匀栅格来说，采样点的 1邻域包括了与其在各方向的距离都不超过1个采样间距的共26个采样点，访问到的采样点，若标记为 faraway，则将其标记改为active，并计算该点与网格模型间的距离值，若标记不为faraway，则不做处理。

(3) 对第(2)步进行迭代，直到所有标记为active的采样点与网格模型间的距离值都大于一定阈值时截止。在处理曲面重建问题时，阈值设定为均匀栅格的采样间距即可，而在处理曲面等距问题时，阈值则要设定为等距距离加上均匀栅格的采样间距。

采用fastMarching算法来构建距离场，在整个的运算过程之中，被访问到的采样点除了最后一个层次外，均为后续步骤中需要用到的，没有对更多的采样点进行访问，因此，采用fastMarching算法来构建距离场在运算时间上必然比传统算法更快速。

当fastMarching算法加快了运算速度后，精确地计算被访问到的采样点与网格模型间的距离值就成为了另一个关键问题。为了更好地发挥其运算速度快的优点，通常会利用较前的波层中的采样点的位置及对应法向量信息，估算出当前遍历波层中的采样点与模型间的距离值，这样的做法在精度要求不高或迭代层数较少时尚可接受，对于精度要求很高的牙齿模型处理将不能适用。

为了解决采用 fastMarching算法来构建距离场时的精度问题，并在传统的距离场构建算法的基础上加快运算速度，本文以fastMarching算法为基础，在计算采样点到网格模型间的距离值时，结合传统距离场算法中的计算方法，提出了优化的fastMarching距离场构建算法，该算法不但能保证精度要求，而且具有较高的运算速度，并有利于模型优化时网格变形的调整。

1.2 算法实现及性能分析

本文优化的 fastMarching距离场构建算法实现的伪代码如下：

BEGIN(算法开始)

//更新 active堆栈(原 active堆栈已加入 alive堆栈，现active堆栈为遍历过程中寻找到的)

14 while(当 active堆栈中所有顶点与模型距离超过某一阈值时终止)

}

END(算法结束)

其中，计算采样点与网格模型间的距离值时，采用了一种快速获得距离值的方法[13]。对于某个采样点v，将三角面片表示成 T(s,t)=B+sE0+tE1的形式，其中 B为三角形的任意一个顶点，设定为基本点，E0、E1分别为另外两个顶点与基本点之间的向量。当满足s,t∈[0,1],S+t＜1时，T(s,t)表示的点处于该三角形之内。通过计算v与T(s,t)之间的距离，可得到一个关于s,t的函数，对该函数取偏导数并求极值，即可算出三角形T(s,t)所在平面上与点v最近的点(即垂足)所对应的s,t的值。需要注意的是，这里计算出的S,t的值不一定满足s,t∈[0,1],S+t＜1，即点v并不处于三角形T(s,t)的正上方，此时，还需根据s,t值的情况进行分类并调整，将s,t的值按情况增减至[0,1]范围以内，然后，将求得的s,t的值带入到式T(s,t)中，就能得到三角形内与采样点v最近的点坐标。进而，根据空间两点间距离公式求得距离值。

针对A、B、C三种牙齿模型进行了算法比较，如表1所示。

表1 使用不同方法的耗时对比

从表 1可知，对于同一种模型，本文算法较传统算法耗时降低 24%以上，速度平均提高25.92%。同时，对于同一模型，两种方法使用的控制点数和三角面片数相同，在计算中由于本文算法第一步即对采样点进行了active和faraway的过滤，因此在计算中实际调入内存中的控制点数和三角面片数较传统算法要少，所以本文算法的内存空间占用也较传统算法更少。

模型A～C分别为采用两种算法获得的牙齿模型，如图 1～3所示。可见本文算法在降低内存占用、提高速度的基础上仍能保证一定的模型精度(优化的fastMarching算法使距离场更精确)。

图1 模型A等距模型效果

图2 模型B等距模型效果

2 CAE分析

2.1 CAE分析原理及结果

采用以上建模算法得到的牙齿CAD模型，由于原模型的准确性以及等距算法的精度问题，使得牙齿模型未必能与原牙完全一致，并与周围牙齿之间保持很好的匹配度和密合度。在治疗中，对于初步制作的牙冠模型，一般通过将其装入原齿位置后的人为咬合运动，利用咬合纸检查修复体与周围牙之间的密合度问题，对于留下咬合纸印记较重的区域再通过手工修磨进行校正，其间患者承受很大痛苦。随着CAE技术的广泛应用，特别是在制造业中的应用，用仿真引导设计的思想已经成熟。

为在义齿实体模型制作之前设计出满足个体需求和精度要求的三维CAD模型，将CAE技术引入到牙齿的模型设计中来，对采用以上算法建立的牙齿模型或牙冠模型，进行有限元分析，利用分析得到的应力、应变分布云图指引设计的修正，即将手工修磨提前到设计阶段，在计算机中实现模型修整，从而可避免或减少患者痛苦。图4给出了采用本文设计思想后的换牙流程与传统换牙流程的对比，可知大量的调整工作在计算机中实现，既提高了速度，又减少了病人痛苦。

图4 传统换牙方法与基于CAE换牙方法的流程对比

应力有限元分析的原理是通过将连续体离散成有限单元的结构离散化，由于不同单元间节点相互连接，位移相同，选择单元位移模式，以节点位移值来计算单元内所有点处位移值的插值函数形成形函数，再通过多个单元形函数的联系，逼近真实解。利用最小势能原理、变分原理、几何方程和物理方程建立起单元刚度矩阵，获得单元节点力与节点位移之间的关系。再组合形成整体刚度矩阵，计算等效节点力，求未知节点位移。

鉴于应力有限元分析在本研究中的作用，为减少计算工作量，仅考虑应力、应变两个物理参数建立本构方程，且认为两者呈线性关系，遵循虎克定律。假设牙齿模型的材料符合各向同性线弹性材料的特性，即仅有两个独立的材料参数，分别为弹性模量E和泊松比v。弹性区域内，应力和应变之间满足以下虎克定律。

将应力放在等号左端，应变放在等号右端可改写成：

式中，εm称为体应变；eij为应变分量；σij为应力分量；下标 i,j分别取值1,2,3，依次代表 x,y,z。

根据总体平衡条件及总能量关系建立物理方程，再通过离散后在每个单元内采用分片插值求得各个单元和节点的应力和应变。根据应变和位移的关系：可积分求得 x,y,z方向上的位移u、v和w。

具体分析方法：在前述方法建立的牙齿CAD几何模型上，定义约束条件及接触面等物理特性，从而建立有限元分析模型，计算得到牙齿或牙冠接触表面的彩色应力分布云图以及各单元节点的应力、应变值，进而利用其指导牙齿或牙冠设计模型的修整。建立满口牙齿模型，目标是使分析所得各个接触位置的应力值趋于均匀，对较高区域进行不断修正从而达到这一目标。

目前满口牙齿模型正在构建之中，本文仅对单齿模型进行应力有限元分析，旨在介绍利用有限元分析结果指导牙齿设计模型修整的方法。图5给出了在牙根部施加固定约束，在牙冠3个高点施加一定载荷时的应力有限元分析结果。由此可见，不同位置的不同应力对应不同的颜色且有值。因此，可根据应力、应变关系建立整体均匀化的约束条件，通过对局部网格模型的调整实现模型的修整优化。

图5 应力分析结果及有限元分析模型

对于牙齿模型进行应力有限元分析时采用 10节点六面体单元，经分析计算可得到各个有限元单元节点上的应力，应力矢量有方向和大小，每个节点有坐标位置，因此，可将有限元网格模型上计算得到的各个节点应力映射到建模得到的三角面片网格模型上，再通过力的分解获得每个三角面片上的力，指导三角面片网格模型的调整。由于牙齿模型表面的复杂性，凹曲面及面片间的联动问题还有待进一步研究。图 6(a)给出了单个10节点六面体单元上10个节点的应力示意，并给出了相邻2个体单元构成模型表面的2个面上节点的应力情况。图 6(b)给出了 2个实体单元面上节点应力对应在三角面片网格模型上的情况。根据力的合成可由实体单元面上节点应力计算出该单元表面上的应力情况，然后采用按面积加权的方法计算出对应三角网格模型各三角面片上应力的大小和方向，用作三角面片网格模型调整优化的指导参数。设三角网格模型某个三角面片 P的面积为 AP，与该三角面片相关的单元表面数目为n，相关单元表面的面积分别为 Am，则该三角面片分别与相关单元表面相交的面积为APAm，第m个相关单元上的应力为σijm，则该三角面片上的各应力分量σijP可按以下公式求得：

图6 力的示意

式中，i,j分别取值为1,2,3，依次代表 x,y,z方向。

最后，由应力分量求出合成应力及其方向，即得三角面片上应力的大小和方向。

2.2 网格模型优化

在制作假牙或全冠修复体时，有可能超出患者牙齿原来的大小，也就是说比患者原先健康的牙齿所占的空间更大，这样患者佩戴后就感觉不适或无法正常咬合。因此，牙医在制作牙修复体前，最好能在计算机中，对修复体模型和装入的对颌牙模型进行观察，分析出双方是否有相交的部分，若存在相交，则应对相交部位的修复体模型进行修整。

判断修复体模型与对颌牙模型之间是否存在咬合相交及其部位的方法：①计算对颌牙模型中的各控制点与全冠模型之间的距离，若修复体模型与对颌牙模型之间没有相互咬合的部分，则对颌牙模型的各控制点均应处于修复体模型的外部；如果在对颌牙模型之中存在某些处于修复体模型内部的控制点，则控制点就处于两者咬合的部位。②基于满口牙齿模型应力有限元分析的结果，应力、应变应比较均匀，反之出现应力、应变峰值的区域即为两者的咬合相交部位。

对于存在咬合相交的情况，要对修复体模型进行调整。首先要对咬合的程度进行分析，若咬合区域很小，所涉及修复体模型的控制点不多，则可以仅对这些控制点做微调处理，即沿着其法向量的逆方向位移其与对颌牙之间的距离即可；而如果咬合区域较大，所涉及修复体模型的控制点多，再使用以上的方法将出现大量的翻转和自交现象，上述方法也就不再适用，还需要进行更好地处理。为此本文提出基于应力有限元分析所得节点应力，及前述三角面片网格模型控制点受力的计算方法，计算得到相交区域面片控制点的受力情况，根据相关面片控制点的受力大小及方向，决定各控制点的调整距离和方向，从而实现非全凸曲面、不规则曲面的动态调整。

3 结 束 语

提出了基于扫描点云数据进行快速建模的优化fastMarching算法，并引入CAE技术对牙齿模型进行应力有限元分析，设计提出了有限元分析实体网格模型应力结果与三角面片网格模型面片及控制点受力情况的映射模型，用于准确指导不规则复杂三角网格模型的调整优化。以牙齿模型为例，分析了本文提出优化算法的优越性，以及由三角网格模型生成实体模型并进行应力有限元分析的可行性和指导意义。下一步将继续研究基于应力有限元分析结果的三角网格模型局部调整的优化算法问题，进一步完善有限元分析模型，建立更具现实意义的分析模型并验证。

[1] Williams R J, Bibb R, Eggbeer D, et al. Use of CAD/CAM technology to fabricate a removable partial denture framework [J]. The Journal of Prosthetic Dentisty, 2006, 96(2)： 96-99.

[2] Bibb R, Eggbeer D, Williams R. RapidManufacture of removable partial denture frameworks [J]. Rapid Prototyping Journal, 2006, 12(2)： 95-99.

[3] Eggbeer D, Bibb R, Williams R. The computer-aided design and rapid prototyping fabrication of removable partial denture frameworks [J]. Proceedings of the Institute ofMechanical Engineers, Part H： Journal of Engineering inMedicine, 2005, 219(3)： 195-202.

[4] 吕培军, 李彦生, 王 勇, 等. 国产口腔修复CAD-CAM系统的研究与开发[J]. 中华口腔医学杂志, 2002, 37(5)： 367-370.

[5] 孙玉春, 吕培军, 王 勇. 用于全口义齿计算机辅助设计的虚拟半可调(牙合)架[J]. 北京大学学报： 医学版, 2008,40(1)： 92-96.

[6]Sun Yuchun, Lv Peijun, Wang Yong.Study on CAD&RP for removable complete denture [J]. ComputerMethods and Programs in Biomedicine, 2009, 93(3)： 266-272.

[7]Maekawa T. An overview of offset curves andSurfaces [J]. Computer-Aided Design, 1999, 31(3)：165-173.

[8] Chen Yong, Wang Hongqing, Rosen D W, et al. A point-based offsetting method of polygonalMeshes [EB/OL]. [2014-06-01]. hppt：//www-bcf.usc. edu/～yongchen/.

[9] 钮叶新, 戴 宁, 袁天然, 等. 基于隐式曲面的三角网格模型等距算法[J]. 中国制造业信息化, 2007, 36(3)： 57-61.

[10] Lorensen W E, Cline H E.Marching cubes： a high resolution 3DSurface construction algorithm [J]. Computer and Graphics, 1987, 21(4)： 163-169.

[11] Ju Tao, Losasso F,SchaeferS, et al. Dual contouring of Hermite data [J]. Computer and Graphics, 2002, 21(3)：339-346.

[12]Sethian J A. LevelSetMethods and fastMarchingMethods [M]. Cambridge： Cambridge University Press, 1999： 86-99.

[13] Eberly D. Distance between point and triangle in 3D [EB/OL]. [2013-06-01]. http：//www.geometrictools. com/Documentation/DistancePoint3Triangle3.pdf, 2008.

Modeling and Optimization Analysis of Irregular Complex Dental Crown

YanShengzan, Zhang Jian, Su Benyue
(School of Computer & Information, Anqing Teachers College, Anqing Anhui 246011, China)

From points cloud data ofScanned tooth, the optimized fastMarching distance field construction algorithm is used to achieve rapid generation of high-precision dental crownModel. And by using CAEStress finite element analysisMethod the force of each element and node is gained, the force relationship between finite elementSolidMesh and triangularMeshModel patch is established. UsingSimulation-guiding-designMethod theManual processing in normal clinical treatment is digitalized and virtualized. Thus theModelingSpeed is effectively increased and theModel recondition is advanced to the designing period. Through thisMethod theModel is created and optimized based on CAE, which canSimplify the treatment process,Shorten the treatment period and reduce patients′Suffering.

irregularMeshSurface; points cloud data; finite element analysis; distance field; decomposition of force vector

TP 391.7

A

2095-302X(2015)02-0198-07

2014-10-08；定稿日期：2014-10-24

国家自然科学基金资助项目(61340016)；安庆师范学院青年科研基金资助项目(KJ201220)

闫胜昝(1978–)，女，河北辛集人，讲师，博士。主要研究方向为计算机技术应用、CAD/CAE等。E-mail：hebeiysz@163.com通讯作者：苏本跃(1971–)，男，安徽芜湖人，教授，博士。主要研究方向为CAD、CAGD等。E-mail：bysu@aqtc.edu.cn