基于骨架驱动的牙齿形态设计

黄香，程筱胜，戴宁，闫国栋，董光雷

(南京航空航天大学机电学院，江苏南京 210016)

随着科技的飞速发展，计算机辅助设计(CAD)与计算机辅助制造(CAM)技术进入了口腔修复学领域，这不仅提高了制作口腔修复体的效率，而且也减轻了患者的痛苦。目前已有以下著名的口腔CAD/CAM修复系统［1］:德国 Sirona公司的 CEREC 3D系统、德国Kavo公司的Everest系统和丹麦3Shape A/S公司的3SHAPE系统等。

将空间几何变形应用于牙齿形态的设计是这些系统中的一个关键技术。国内的口腔CAD相关研究一般集中在北京大学口腔医学院［2］、第四军医大学［3］等口腔医学院，研究的平台还停留在利用逆向工程软件作一些面向应用的二次开发，缺少核心技术。

1984年Barr［4］首次提出了整体和局部变形的思想。该变形手段只适用于特定的几类几何形状，要产生任意的、比较复杂的形状变形则有很大的难度，不适合牙齿表面形态的设计。1986年Sederberg等［5］提出了自由变形(FFD)的概念。该方法适用于柔性物体动画，变形范围广，但工作量较为繁重。Celniker等［6］提出了基于物理模型的曲面造型方法，Lazarus等［7］则提出了轴变形(AxDf)的方法。依据所使用变形工具的不同，可将上述变形技术分为4类［8］:基于体的变形、基于曲面的变形、基于曲线的变形、基于点的变形。其中，曲线驱动变形技术具有多方面优势，既操作简单(点驱动变形优势)、易于实现，同时也可以实现比较复杂的变形(体、面驱动变形优势)。本研究中的骨架线由参数曲线来表示。

Lazarus等［7］在构建标架时，若曲线某一点处的曲率不存在，则该点处标架构建失败。建立模型与曲线的映射关系时，可以选择模型上的点与线上的最近点建立关系。这些最近点可以用曲线的凸包性和Oslo算法［9］确定，也可以使用离散样条曲线来确定。基于上述研究，作者提出了一种快速实现牙齿修复体表面形态设计的方法。

1 骨架驱动模型变形原理

骨架线驱动模型变形如图1所示，变形原理为:先为模型M构造近似骨架C，建立近似骨架C与模型M之间的映射关系;移动曲线上的控制点使曲线从C变形到C';依据曲线在离散点处的新局部坐标，更新模型上的点坐标。这一变形过程可以用一个简单的关系式:F:Ω3→C→Ω3来表示［10］。

图1 轴线驱动变形Fig 1 Axial curve drive deformation

整个变形的流程如图2所示:

图2 变形流程图Fig 2 Flowchart of deformatio n

在建立映射关系时，模型需要被划分成一个个子空间，这些子空间的连续性由曲线的连续性来保证。由于曲线是连续且光滑的，并且模型的子空间被划分的很密，因此，可以近似地认为这些子空间也是连续且光滑的。在变形时，可以保证模型的曲面具有较好的光顺性。

2 “骨架线”的构建

参数曲线的表现形式有:贝塞尔(Bezier)曲线、均匀B样条曲线和非均匀有理B样条(NURBS)曲线等。其中Bezier曲线缺乏局部特性;均匀B样条曲线虽然局部性较好，但是没有保留Bezier曲线端点处的几何性质;NURBS曲线同时可以克服上述的缺点，但由于引入了权因子，一旦权因子选取的不合适，可能会导致很坏的参数化，从而使问题变得复杂化［11］。综合多方面的因素，我们采用3次准均匀B样条构造参数曲线，它具备良好的端点性质和局部性质，支持曲线局部编辑，便于对模型的局部区域实施变形。

勾画曲线的交互手段类似于文献［12］中所使用的方法。给定一个原始牙齿模型，用户依据牙齿的表面形态，对该模型进行勾画，选定兴趣区域，这个兴趣区域被参数曲线首末两控制点构成的切平面所限定。勾画时，用户可以自由地调节屏幕视点，使操作更加便利。勾画完成后所得到的曲线位于模型的表面。用户可以通过平移、旋转等操作调整曲线的位置，使其控制点落于模型之内。把控制点位于模型内部的三维(3D)参数曲线定义为“骨架线”。

3 网格变形

3． 1 定义局部标架

为明确映射关系，首先需要在参数曲线(“骨架线”)上建立局部坐标系。在建立坐标系之前，需要将曲线均匀离散化，并把离散点数据存储到容器中。构造曲线上的局部活动标架有多种方法。其中弗朗内特(Frenet)标架是一种最常见的活动标架，但该标架在曲线曲率消失的地方无法进行定位，而且，在曲线的拐点处，该标架的主法矢和副法矢有可能会反向或产生不必要的旋转，进而可能会导致模型的暴力扭曲［13］。为避免上述的问题发生，采用如下方法来建立坐标标架。

将活动标架F表示为一个四元组(t，T，N，B)，其中t为标架原点对应曲线C的差值参数，T=C'(t)/‖C'(t)‖，为C上参数t点处的单位切向量，‖·‖为欧式范数，N、B也是标架原点处的单位向量，T、N、B两两正交构成F的三个坐标轴。

首先构造一个t=0处的起始标架。令T0为该处的单位切矢T0=C'(0)/‖C'(0)‖，然后将T0投影到XY平面得到T0p，将T0p旋转90°后再正规化即为N0，如图3所示。B0则由T0和N0的叉积计算得到，即B0=T0×N0。如果T0的方向与世界坐标系下的向量(0，0，1)的方向一致，则该计算过程将会失效。为了解决这个问题，我们检验T0及其投影向量的合理性，选择世界坐标系下的任意一个平面(XY、YZ或者XZ)，保证T0投影到该平面之后的向量不是一个零向量。

图3 起始标架计算Fig 3 The creation of the first local frame

沿曲线在t∈［0，1］的定义域内构建其余标架。已知当前t=ti处曲线上点的坐标为Pi、单位切矢Ti，位于该点处的平面 πi。将向量 Ni－1投影到平面 πi上，并将其正规化即为Ni，则第三个向量Bi为Bi=Ti×Ni。在计算过程中要保证 α ＜90°，α 为 Bi与 Bi－1之间的夹角大小即这种算法简单高效，避免了曲线暴力扭曲的发生。

3． 2 划分区域

在变形之前，需要构建曲线与原始网格上的所有点之间的映射关系。我们采用文献［14］的方法来确定网格与曲线之间的关系。如图4所示，〈·|·〉表示两个向量之间的点积，vm表示网格上的点，Pi表示曲线上点的坐标即局部坐标原点，Ti即单位切矢，Pi和Ti确定平面πi。划分区域时，所用的判定式定义如下:

图4 模型点、标架原点与单位切矢位置关系Fig 4 The relationship among the model vertex，the origin of frame and the unit tangent

依据上述判别公式，程序每循环1次，需要遍历的网格顶点数就减少一部分，这就有效地提高了计算速度。我们为局部标架及其对应的网格点索引值建立一个结构体，用C++语言定义如下:

struct Pt_Mesh

{

std::vector＜HPoint＞base_pt;/*表示离散后的曲线上的点的坐标即 Pi，HPoint 是一个三维坐标点类*/

std::vector＜int＞ index; /*Pi所对应的那一块网格数据点的索引值*/

std::vector＜float＞ind_rm;/*网格上的每个点所对应的一个轴向参数，在变形过程中保持不变，会在接下来的章节中进行描述*/

};

3． 3 模型变形

定义 πk－1和 πk分别是由向量(Nk－1，Bk－1)和(Nk，Bk)所张成的空间平面。位于这两个平面之间的点集vj可以由 πk－1平面上的标架来确定它的局部坐标。

设OXYZ为全局坐标系，位于平面πk－1和πk之间的网格上的点为 vm，Pk－1为标架 Fk－1的原点，而 Tk－1、Nk－1、Bk－1构成了该标架的 3个坐标轴，则 vm在局部标架Fk－1下的坐标为(¯x，¯y，¯z)，它们之间的变换关系如下:

这样就建立了全局坐标系与局部坐标系之间变换公式。令vm沿着Tk－1的方向投影到平面πk－1和πk上的点分别为vmp和vmn，vm和轴向参数rm通过等式

关联起来，如图5所示。dist是两个点之间的欧式距离，这个参数在变形过程中保持不变。

图5 网格点vm分别投影到面πk－1和面πk，得到vmp和vmnFig 5 vmand its projections vmpand vmnon πk －1and πk，respectively

当曲线C变形后，曲线上的离散点处的标架将进行重新计算。设变化后的曲线为C'。由于在计算曲线C的起始标架时，T0投影到世界坐标平面时具有随机性，若采用相同的方法来计算C'的起始标架，则有可能会出现两次标架差别很大的情况。因此，我们采用如下的方法来计算变化后的起始标架［9］。曲线C在 t=0 处的标架为(T0，N0，B0)，T'0为变化后的曲线在t=0的单位切矢。令ω为T0和T'0之间的夹角，向量R=T0×T'0，N0和B0沿着R的方向旋转ω角度得到N'0和B'0。起始标架计算完成后，其余标架的计算方法同原始曲线在t=tk处计算标架的方法。从而vm新的坐标值为

但是这一公式并没有考虑到控制点在移动曲线上时，曲线产生轴向拉伸的情况。为了把这种拉伸效果传递到网格，将v'm的公式变形为

其中，v'mp为变化后的网格点在 π'k－1上的投影，v'mp与向量T'k－1所确定的直线与平面π'k的交点为v'mn，如图5所示。

曲线上标架的数量决定了模型变形结果的平滑程度。如果标架的数量较少，则变形效果有可能令人不十分满意。但如果标架数量过多的话，又会延长计算时间。因此，我们需要在两者之间找到一个平衡点。

4 实验与分析

开发实现平台:Windows XP系统，Visual C++．net，Hoops图形显示包。作者通过切牙、磨牙的整体变形和双尖牙的牙尖变形来验证算法的有效性与可行性。由于曲线采用了3次准均匀B样条，因此，曲线上的某个点发生变化时，其附近区域也会有相应的改变。图6～8分别表示3类牙的变形效果。

4． 1 切牙整体变形

图6 切牙整体变形Fig 6 The global deformation of anterior tooth

图6显示了切牙整体旋转的过程:a图是切牙原始模型;b图是调整后的参数曲线，位于原始模型内部;c图是变形后的切牙牙齿形态，其中虚线表示原始“骨架线”，实线表示变化后的“骨架线”，箭头方向表示变形方向。

4． 2 磨牙整体变形

图7 磨牙整体变形Fig 7 The global deformation of molar tooth

图7显示了磨牙整体拉伸的过程:a图是磨牙原始模型;b图是磨牙调整后的“骨架线”，位于原始模型内部;c图是拉伸以后的磨牙牙齿形态，箭头方向表示变形方向。

4． 3 双尖牙局部变形

图8 双尖牙局部变形Fig 8 The local deformation of bicuspid tooth

图8显示了双尖牙某一牙尖局部拉伸的过程:a图是双尖牙原始牙尖形态;b图是对变形区域进行曲线勾画，以及调整参数曲线后的形态，曲线位于模型内部;c图是变形后的双尖牙牙尖形态，其中虚线表示原始“骨架线”，实线表示变化后的“骨架线”，箭头方向表示变形方向。

选取切牙整体绕X轴顺时针旋转10°作为考察对象，如表1所示。

表1 切牙变形结果的比对(绕X轴旋转10°)Tab 1 The comparison of anterior tooth's deformation results

由表1我们可以看出:当标架数量增多时，划分区域的时间、局部标架时间和变形时间都会相应变长，但变形效果越来越好。当标架数量为200和500时，牙齿表面形态的效果相差无几。考虑到程序的计算速度以及变形效果的好坏程度，选取每条曲线标架数量为200。

5 结束语

现今的口腔CAD/CAM修复，一般是调用数据库中的标准牙冠作为修复体，该修复体是一个统计模型，尽管具备了牙齿形态的大部分特征，但每一个患者的牙齿形态特征还是存在一定差异性的。我们依据邻牙、颌牙牙面的约束关系，可以对修复体进行整体变形，也可以对修复体的局部特征(例如牙尖)进行修改。我们提出的基于骨架线驱动的牙齿形态设计其优点如下:(A)在牙齿表面勾画曲线时，可以随意地删除、增加曲线上的控制点，若对整体都不满意则可以进行重绘。这提高了操作的交互性与简便性。(B)在建立骨架上的局部标架时，我们所使用的方法克服了Frenet标架的缺点，骨架线既可以是直线也可以是曲线。也就是说，即使线上的曲率不存在也可以在该点处定义标架;同时也克服了因曲线出现暴力扭曲而导致的不满意的牙齿形态变形。(C)牙齿模型区域划分时所用的算法计算简单，实用性较好，可以有效地提高计算效率。(D)牙齿模型上的点从局部坐标映射回全局坐标时，我们考虑到曲线变形会出现轴向拉伸的情况，因此，为了使牙齿形态的变形具有更好的光顺性，对等式(5)做了一些改进，进而将这种拉伸传递到网格。

该方法还有一些问题亟待解决，如多条参数曲线如何实现控制变形，以及在变形过程中如何通过相应算法消除局部自交。这些将是我们今后研究的重点。

航空航天大学，2006．

［2］吕培军，李彦生，王勇，等．国产口腔修复CAD-CAM系统的研究与开发［J］．中华口腔医学杂志，2002，37(5):367-370．

［3］金树人，姚月玲，高勃，等．应用快速成型法制作磨牙树脂全冠［J］．第四军医大学学报，2003，24(8):700-702．

［4］BARR A H．Global and local deformations of solid primitives［J］．Computers Graphics，1984，18(3):21-34．

［5］SEDERBERG T W，PARRY R．Free-form deformation of solid geometric model［J］．Computer Graphics，1986，20(4):151-160．

［6］CELNIKER G，GOSSARD D．Deformable curve and surface finite-elements for free-form shape design［J］．Computer Graphics，1991，25(4):257-266．

［7］LAZARUS F，COQUILLART S，JANCENE P．Axial deformation:an intuitive technique［J］．Computer Aided Design，1994，26(8):607-613．

［8］徐岗，汪国昭，陈小雕．自由变形技术及其应用［J］．计算机研究与发展，2010，47(2):344-352．

［9］魏斌，袁修干．基于NURBS曲面的轴变形方法［J］．北京航空航天大学学报，1997，23(5):546-550．

［10］JAMES G，DOMINIQUE B．A survey of spatial deformation from a user-centered perspective［J］．ACM Transactions on Graphics，2008，27(4):1-21．

［11］施法中．计算机辅助几何设计与非均匀有理B样条［M］．北京:高等教育出版社，2001:309．

［12］KHO Y，GARLAND M．Sketching mesh deformations［C］//Proceedings of the 2005 Symposium on Interactive 3D Graphics．Washington D C:ACM，2005:147-154．

［13］ BLOOMENTHAL J．Graphics Gems［M］．San Diego:Academic Press Professional，1990:567-571．

［14］FORSTMANN S，OHYA J，KROHN-GRIMBERGHE A，et al．Deformation styles for spline-based skeletal animation［C］//Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on computer animation，San Diego:Eurographics Association，2007:141-150．