一种基于牙体长轴的质点弹簧模型的牙龈变形仿真方法与流程

本发明属于图像处理领域，涉及一种改进的简化质点弹簧模型在牙龈变形中的仿真。

背景技术：

近年来，随着3D打印技术和数字几何处理技术的不断发展，三角网格模型成为三维模型的主流表示方式，继而如何对三角网格模型进行交互式地实时仿真，一直是研究人员关注的一个研究热点。在变形技术发展的几十年间，为使得变形仿真系统更加稳定、真实、实时、精确，学者们提出和延伸出数十种变形仿真方法，许多学者根据不同的角度对其进行分类，主要分为几何方法和物理方法。

其中几何方法中的自由形变技术由Sederberg和Parry于1986年提出，其主要思想是直接通过控制点来改变物体的轮廓形状(Sederberg T.W，Parry S.R.Free-form deformation of solid geometric models[J].Computer Graphics，1986，20(4)：151-160.)，是一种经典且强有力的算法，从某种程度上被认为是建立在传统曲线和曲面几何造型理论基础上的。这种方法有明显的缺点：交互过程不直观，灵活性高不易控制，变形结果不易达到预期结果，而且当控制点的数目较多时，会加大时间开销。蒙皮法即基于骨骼的蒙皮变形技术或骨架驱动的变形技术，最初由Thalmann等提出，其基本思想是将三维网格的变形与人体的运动相类比，将面模型的点即皮肤关联到通过编辑得到的骨架结构上，网格的变形由骨架的运动来带动(Magnenat-Thalmann N，Laperrire R，Thalmann D.Joint-dependent local deformations for hand animation and object grasping[C]//In Proceedings on Graphics interface’88.1988.)。该方法操作简单直观，适用于运动变形仿真，并且编辑的骨架模型，一次建立之后就能被反复使用，因此被广泛应用于人体动画、游戏等领域。但由于几何方法不考虑组织的力学本构关系，缺乏生物特性，且其需构造矩阵等式及求解造成实时性差等原因，以及计算机计算能力的提高，使得基于物理的方法成为近年来研究的热点。

现在有多种较为成熟的用于模拟软组织的变形的物理方法，它们大致可以分为两类：一种是基于网格结构的变形方法，主要有质点弹簧模型、有限元模型和边界元模型；另一种是不依靠网格的方法，主要有粒子系统，有限球方法，光滑粒子流体动力学方法以及无网格方法(郭煜，秦璟.虚拟手术中软组织形变仿真研究综述[J].集成技术，2013(02)：52-61.)。其中，质点弹簧模型和有限元模型较为常用，这两种方法具有简单、灵活拓扑结构的质点弹簧模型最显著的优势就是速度很快，视为可满足实时性的要求，结构简单可根据需要在模型上进行改进，加入几何修正也不是难事，因此质点弹簧模型在变形仿真中有广泛地应用，但其稳定性不佳，另外，在标准方程中不能保持变形模型的体积不变，因为软组织几乎都是不可压缩的，所以这将会产生坍塌不真实的变形结果。高精度的有限元模型的优势在于可以解决各种不同的复杂情况，因为它相对于网格的几何结构以及外部载荷都是独立的，在描述变形系统时，只需要较少的材料参数，而且它还有坚实的理论基础。但其计算量庞大，时间较长，难以达到实时性，尤其是对于非线性有限元。

技术实现要素：

在经典的质点弹簧模型中，一般用三类弹簧，结构弹簧、剪切弹簧和弯曲弹簧，分别表现结构应力，剪应力和弯曲力。但已知牙龈的变形区域是不规则三角网格，不适合将这三类弹簧全部使用，因此有使用简化的质点弹簧模型或者只有剪切力的弹簧，但这种模型也不适用不规则三角网格。

针对现有技术的不足，本发明提出了一种基于牙体长轴的质点弹簧模型的牙龈变形仿真方法。

本发明所采用的技术方案为：

一种基于牙体长轴的质点弹簧模型的牙龈变形仿真方法，包括如下步骤：

1)确定牙龈模型；

2)计算出牙体长轴的位置；

3)存储质点的初始位置、邻接点以及它对应于牙体长轴的映射点等信息；

4)确定每一个时间步长内即每次迭代中质点的位置、速度和受力；

5)当拖动牙冠时，根据移动矩阵计算牙体长轴的新位置；

6)判断停止条件，输出仿真结果。

所述步骤(1)中的牙龈模型是根据过渡区域的网格曲率差，生成一条牙齿与牙龈的分割线，作为变形模型的主动点，向外扩展N层(优选地，N＝8)的面片组成变形牙龈区域为形变区域，最外层的小球组成约束点，从而形成牙龈模型。

所述步骤(2)中的牙体长轴是一条垂直穿过牙体即牙冠和牙根中心的假想直线。

所述步骤(3)中将网格模型中三角面片的顶点作为质点，边作为连接质点的弹簧；再建立变形区域网格模型顶点与牙体长轴之间的映射，且该映射关系在变形过程中保持不变；将质点向牙体长轴垂直投影的点作为弹簧的另一端，从而建立牙龈与牙体长轴之间的体弹簧。

所述步骤(4)计算出质点的位置、速度和受力；根据隐式欧拉求精确值的质点弹簧模型快速模拟方法，结合已知质点在t时刻在位置x_t、速度v_t和受力f_t，计算出质点在t+h时刻的位置x_t+h、速度v_t+h和受力f_t+h，其中h为时间步长；则隐式欧拉的表达式为：

x_t+h＝x_t+hv_t+h

同理由隐式欧拉的表达式可得：

v_t+h＝v_t+hm^-1f_t+h

两边同乘以h得hv_t+h-hv_t＝h²m^-1f_t+h，和上式对比得：x_t+h-2x_t+x_t-h＝h²m^-1f_t+h则受力为：

所述步骤(6)中的停止条件是质点的速度或位移量小于给定值；当停止条件不满足时，则计算体弹簧的形变量。

相对于现有技术，本发明的优点在于：

根据牙龈模型的特点选用简化的质点弹簧模型，即模型中只有结构弹簧。在简化的质点弹簧模型上进行改进，在牙龈和牙根之间增添了一种弹簧——“牙根弹簧”，用牙体长轴代替牙根与牙龈面模型建立联系，牙体长轴作为弹簧的一端，牙龈面模型上的点为弹簧的另一端，由此构成体弹簧。由牙根的移动来驱动弹簧的拉伸和压缩，进而模拟牙槽骨对牙龈变形的影响，来达到含有体信息的仿真效果，构建出一个基于牙体长轴的质点弹簧模型。对待弹性参数时，本发明对按弹簧原始长度的大小进行分级设置，弹簧原始长度越大，特性系数越小，反之，系数越大。

为了证明本发明的可行性，进行了试验仿真，通过试验结果表明，添加了体弹簧的简化质点弹簧模型，相较于一般的简化质点弹簧模型在变形效果上更加饱满，同时系统能够更快达到迭代停止条件，总耗时减少。

附图说明

图1为本发明的流程图

图2为单牙分割线示例图

图3为建立标架示例图

图4为牙龈变形区域图

图5为牙龈变形区域放大并线框显示图

图6为三种类型牙根图

图7为牙体长轴图

图8为牙体长轴矢量示例图

图9为牙体长轴示例图

图10为P_i的Laplacian坐标图

图11为网格顶点与牙体长轴的映射图

图12为简化的质点弹簧模型方法和本发明方法的对比图

图13为简化的质点弹簧模型方法和本发明方法的时间消耗对比图

具体实施方式

系统读入牙颌模型的STL文件后，计算并得到牙龈分割线，将单颗牙齿和牙龈进行分离。本发明是在分割后的模型上进行牙龈变形方法的研究。如图1所示，具体建模步骤如下：

步骤(1)确定牙龈模型。

研究牙龈变形方法一般都是对单颗牙齿进行操作，因此本发明简化牙齿分割操作部分，只对单颗牙齿进行分割，先选中需矫正的牙齿，再选中与它相邻的牙齿，如图2所示，这里识别出相邻牙齿是为了在下一步扩展生成牙龈变形区域时，不会把相邻牙齿的面片也包括进去。之后生成分割线，调整分割线到理想的位置，分割牙齿、牙龈，并对先选中的牙齿即矫正牙齿建立标架，用于之后对牙齿的旋转操作，它是由3条相互垂直的长轴和分布在长轴两端的6个小球组成，如图3所示。

最后，变形区域如图4-5所示。图4中牙齿与牙龈边界的小球是由上一步骤中的分割线获取的，作为变形模型的主动点，向外扩展N层的面片组成变形牙龈区域为下文的形变区域，本发明中N设为8，最外层的小球组成约束点。从图5可以看出这是不规则三角面片组成的，有的顶点有8个邻接面片，有的只有3个。

步骤(2)计算出牙体长轴的位置。

牙根的情况比较复杂，有三种：单根、双根和三个根，如图6所示，虽然牙根的外观并不规则，但在口腔医学上，为了方便描述，认为它们都有一条牙体长轴。牙体长轴是一条垂直穿过牙体即牙冠和牙根中心的假想直线，如图7所示。

本发明用牙体长轴代替牙根与牙龈面模型建立联系，牙体长轴作为弹簧的一端，牙龈面模型上的点为弹簧的另一端，由此构成体弹簧。在求解中，若用牙冠中心和牙龈曲线的中心之间的连线表示，容易造成牙根倾斜，因为这两个点之间的距离比较近。因此，可以用牙冠中心点和牙龈中心点的连线表示，或是从构造牙根的算法中获取。在实验部分，将牙体长轴定义为从牙冠中心指向牙龈变形部分中心的矢量，如图8所示。图9中红色的线段就是计算得到的牙体长轴。

步骤(3)存储质点的初始位置、邻接点以及它对应于牙体长轴的映射点等信息。

Laplacian坐标可以很好的存储网格的拓扑信息和几何信息，在三角网格变形前，先对网格顶点的Laplacian坐标进行存储。给定三角网格，其中P代表点的集合。如图10所示。设P₁...P_j是点的邻接点，则微分坐标定义为：

其中，N(i)表示P_i点的邻接点的集合，d_i表示N(i)中点的个数。

根据步骤(1)知，本发明的质点弹簧模型网格就是变形的牙龈区域，将三角面片的顶点作为质点，边作为连接质点的弹簧。当简化的质点弹簧模型建立后，需建立变形区域网格模型顶点与牙体长轴之间的映射，且该映射关系在变形过程中保持不变。将质点向牙体长轴垂直投影的点作为弹簧的另一端，就建立牙龈与牙体长轴之间的体弹簧。这时，映射的垂直距离就是弹簧的原长，需保存下来用于计算变形时弹簧的形变量。

其中需要设置弹簧的弹性参数。

观察图8-9，可以看出牙龈变形区域可以近似看作是一个漏斗样的模型，越靠近牙冠位置的牙龈距离牙根越近，越远离牙冠位置的牙龈距离牙根越远，也可以说在移动牙齿时，远离的牙龈受牙根的影响越小，考虑到这个情况，本文在设置添加弹簧的弹性系数时，按弹簧原始长度的大小进行分级设置，弹簧原始长度越大，特性系数越小，反之，系数越大。

步骤(4)确定每一个时间步长内即每次迭代中质点的位置、速度和受力。

根据胡克定律可以将质点弹簧的运动规律模拟牙龈的变形过程，则每个质点的运动都满足拉格朗日运动方程：

其中x表示质点位置，m表示质点质量，γ表示弹簧的阻尼系数，k_x表示弹簧的刚度系数，其中ε＝l_org-l_cur，l_org表示弹簧的初始长度，l_cur表示弹簧的当前长度，等式的右边f表示质点受到的外力。

由于本发明不考虑重点的因素，故整个质点运动系统满足微分方程组：

已知质点的初始位置，就可以求得X_t+Δt。更新求出每个质点的新位置后，反复迭代，就可得到任意时刻的质点的位置。

一阶常微分方程初值问题的数学模型一般形式为：

根据隐式欧拉求精确值的质点弹簧模型快速模拟方法，结合已知质点在t时刻在位置x_t、速度v_t和受力f_t，计算出质点在t+h时刻的位置x_t+h、速度v_t+h和受力f_t+h，其中h为时间步长。则隐式欧拉的表达式为：

x_t+h＝x_t+hv_t+h

根据隐式欧拉的表达式可得hv_t+h＝x_t+h-x_t和hv_t＝x_t-x_t-h，两式相减得：

hv_t+h-hv_t＝x_t+h-2x_t+x_t-h

同理由隐式欧拉的表达式可得：v_t+h＝v_t+hm^-1f_t+h，两边同乘以h得hv_t+h-hv_t＝h²m^-1f_t+h，和上式对比得：

x_t+h-2x_t+x_t-h＝h²m^-1f_t+h

为了简化表达式，令x_t+h为x，2x_t-x_t-h＝y，同时质点质量用矩阵表达，则

M(x-y)-h²f(x)＝0

其解符合下式的最小值：

而f_l代表体弹簧弹力的势能，L是网格刚度系数Laplacian算子矩阵，是由弹簧两端的质点决定的，J是一个弹簧的系数矩阵，是由系统中所有弹簧所决定的。

求g(x)的最小值：

式中M+h²L是一个稀疏矩阵并对称半正定，令求得最小值的解为：

(M+h²L)x＝My+h²(Jd-f_b-f_l)、

等式中d是通过t时刻的位置x(t)计算得出的，x_i，x_j代表弹簧的两端：

对于以上的数值求解过程中，在迭代前只需对M+h²L进行一次LLT快速分解，之后的迭代求解看作是一个局部线性方程，因此尽管系统中的质点和弹簧比较多，这种算法却能大大减小求解的时间开销。

步骤(5)当拖动牙冠时，根据移动矩阵计算牙体长轴的新位置。

移动牙齿时，牙体长轴的位置也跟着改变。具体计算过程如下：

(1)将变形区域的网格顶点垂直映射到牙体长轴上，如图11所示，牙体长轴上的点S_i即为网格顶点P_i相对应的映射点。S_i与P_i为弹簧的两端，它们之间的距离为弹簧的原始长度。

(2)移动牙齿后，牙体长轴上的映射点均左乘牙齿的移动矩阵，就可以得到移动后的牙体长轴上映射点S_i的位置。

步骤(6)判断停止条件，输出仿真结果。

本发明中的停止条件是质点的速度或位移量小于给定的一个值。当停止条件不满足时，则计算体弹簧的形变量，将体弹簧产生的弹力作为质点下一次迭代所受到的外力的一部分，循环步骤(4)。

用简化的质点弹簧模型与添加体弹簧的改进模型在牙龈模型上进行挤压操作和恢复操作，移动量为1.2，对变形效果和时间消耗进行对比，如图12和图13所示。

在对数值积分部分即每次迭代都需重复一下几个步骤：

(1)判断是否对M+h²L进行了LLT分解；

(2)获得当前时刻各个质点的位置；

(3)根据质点位置计算更新矩阵d以及质点的Laplacian坐标；

(4)计算质点的弹簧弹力势能f_l；

(5)更新步骤(4)中各个变量。

(6)判断是否满足停止条件，是的话停止迭代，否的话重复1～6。

试验分析得出，改进的质点弹簧模型变形效果相较简化的质点弹簧模型比较饱满，但也丢失了原模型的几何细节，这是由这种变形模型的本质所决定的；在时间消耗方面，单次迭代的时间增加0.1～0.2ms，但减少了总的迭代次数，添加的弹簧使系统能够更快达到迭代停止条件，因此总耗时减少。