基于网格抽取的牙齿模型快速分割算法

马 勇 柯永振 杨 帅

(天津工业大学计算机科学与软件学院 天津 300387)

0 引 言

随着社会的快速发展，人们的生活水平得到了很大的提高。与此同时，很多人开始重视自身的外表，而牙齿作为外表的一个重要组成部分，直接影响着个人的仪容和生活质量。由于一般的金属牙套存在形状突出、难清洗等缺点，所以近年来兴起的隐形牙套受到了人们的广泛欢迎。制作隐形牙套主要有牙齿数据获取以及三维重建、牙齿分割、绘制牙弓线、虚拟排牙以及隐形牙套生成等过程，而牙齿分割是非常重要的一个过程。

近十年，在牙齿模型分割领域中，专家和学者提出了许多牙齿分割的算法。文献[1]针对牙龈上的轮廓线，一种基于能量方程的牙齿分割算法被提出，通过求解能量方程，在条件收敛之前，该算法对牙龈上的轮廓线不断地进行迭代并进行移动，这种算法虽然没有进行人工交互，但是对于很多排列不整齐的牙齿并不能进行精确的分割。文献[2]采用洪泛法提取牙齿上的信息，通过这些信息识别牙齿之间的凹区域，即牙齿间的边界，该算法对于排列比较整齐的牙齿有较精确分割结果，但对于排列不规则的牙齿没有良好的分割效果，而且分割效率非常低。文献[3]提出了交互标记控制算法，该算法交互选取种子面片，然后根据相邻两个面的曲率进行区域生长，该算法分割速度较慢，而且存在过度分割现象。文献[4]提出一种基于改进后的交互标记控制算法，该方法在分割过程中采用人为的设定阈值作为控制参数来控制队列的操作，这种算法使牙齿的分割速度得到了提高，但对于一些不规则牙齿并不能进行精确的分割。文献[5]提出一种基于调和场的牙齿分割交互算法，该算法在分割过程中采用调和场的属性提取等值线，进而提取牙齿边界，然后再交互进行牙齿分割，这种方法提高了分割的精度，但对于一些不规则牙齿并不能得到良好的分割效果。文献[6]通过计算牙齿三角网格的曲率以及进行一系列的形态学的开闭操作，提出了一种对于三角网格进行形态学操作的分割算法，这种算法对牙齿可以进行比较精确的分割，但是对于牙齿之间的融合区域，并没有提出有效的解决方法，而且分割的速度较慢。文献[7]针对牙齿之间的融合区域，找到了一种解决的方法，该方法首先剔除牙齿之间的融合区域然后再对没有融合区域的牙齿重新建模，那样做破坏了原有的牙齿网格模型上的拓扑结构。文献[8]对提取的牙齿特征区域进行膨胀和腐蚀操作，然后利用形态学开闭操作对牙齿的特征区域进行细化，将牙齿分离，该方法对排列规则的牙齿有良好的分割结果，而对排列不规则的牙齿，分割结果较差，而且分割效率较低。文献[9]提出了一种基于调和场的自动分割算法，该算法能有效的一次性获取牙齿网格模型上的所有牙齿区域，但该算法并不能进行单个牙齿分离。

本文在交互标记控制算法的基础上，提出了一种基于网格抽取的牙齿模型快速分割算法。该算法主要改进了三个方面：第一，进行网格抽取，可以适当减少网格的数量，加快网格的分割速度；第二，进行网格平滑，消除网格中的噪声，提高网格分割的精度；第三，改变数据结构，取消排序算法，减少算法的运行时间。

1 交互标记控制算法

交互标记控制算法[4]是在分水岭算法[10]的基础上提出来的。该算法遵循Hoffman等[11]提出的极小值法则，并按照Fast Marching Watershed[12]理论进行区域扩展从而确定分割区域的位置和范围[13-14]。

该算法的思想描述如下：

交互标记控制算法首先在每个牙齿上人工拾取一个三角网格，并对其进行标记。接着利用堆栈作为存放三角网格的数据结构，对堆栈进行初始化操作，即把每个牙齿上拾取并被标记的三角网格的相邻的三角网格全部压入堆栈中。然后把堆栈里的三角网格按照相对弯曲程度进行从大到小的排序，每次从堆栈中取出位于栈顶的三角网格，即相对弯曲程度最大的三角网格，将该网格进行标记，并且将与该面片相邻的且没有被标记的三角面片按照相对弯曲程序的大小顺序进入堆栈中。此过程循环进行，直到堆栈为空。至此，牙齿上所有的三角网格都被标记了一次，通过标记，可以确定每个网格属于某个区域，即属于某个牙齿上。

该算法的具体流程为：

(1) 堆栈初始化，即把所有被人工交互标记的三角网格的相邻网格按照曲率的大小压入堆栈中。

(2) 区域生长，即每次从堆栈中取出位于栈顶的三角网格，记为，则将网格face′标记为flag′，并将与该网格face′相邻的并且没有被标记的网格gi、flag′和相应的曲率H(face′,gi)构成的结构体按曲率大小排序，然后顺序压入堆栈中。

2 基于网格抽取的牙齿模型快速分割算法

由于之前的交互标记控制算法分割时间长且存在欠分割或过分割的现象，为了提高算法的分割速度和精度，本文提出了一种运行速度较快，分割精度较高的牙齿模型快速分割算法。

基于网格抽取的牙齿模型快速分割算法思想描述如下：

(1) 网格抽取 对于输入的牙齿三维网格模型，进行网格抽取，适当缩减网格的数量。

(2) 网格平滑 利用拉普拉斯平滑算法对三角网格进行平滑处理，以消除网格的噪声。

(3) vector容器初始化 交互操作，使已经标记的三角网格相邻的网格进入vector容器中。

(4) 区域生长 访问牙齿网格所有面片，根据极小值法则[11]，对网格面片进行分组，进而实现牙齿分离。

基于网格抽取的牙齿模型快速分割算法的流程图如图1所示。

图1 基于网格抽取的牙齿模型快速分割算法流程图

2.1 网格抽取

一般经过激光扫描得到的STL格式的牙齿模型的网格数量非常多，为了能够加快算法的运行速度，本文采用了文献[15]提出的网格缩减算法进行牙齿模型网格抽取，即适当减少牙齿模型网格的数量。多次实验表明，网格抽取10%，即网格数量减少10%时，使用本文提出的算法，既不会影响牙齿网格分割的精度，又提高了网格的分割速度。该算法首先将网格G=中的每个顶点Vi分类，并插入到优先队列S1中。优先级队列是基于误差来删除顶点并重新三角化，无法删除或三角化的顶点将被跳过。接着处理优先级队列S中每个顶点Vi(即使用边折叠来删除顶点，然后进行三角剖分)，继续迭代，直到优先级队列S1为空。最后处理所有剩余的顶点，沿着锐利的边缘将网格划分成单独的片段，并且将其插入到优先级队列S2中。再次处理优先级队列S2，直到队列S2为空。算法的伪代码如下：

网格抽取算法伪代码：

Input: G

//输入网格的几何和拓扑信息

classify Viand Push S1

//识别网格顶点并放入队列S1

While(S1!=empty) do

If(delete Viand re-triangulated )

//删除不符合条件的顶点并且重新三角化

Vitriangulated

End if

Else ViPush S2

//不能被删除的顶点进入队列S2

End while

While(S2!=empty)

Vitriangulated

//剩余的顶点三角化

End while

Return G1

2.2 网格平滑

经过激光扫描得到的STL文件，往往存在着噪声，为了消除噪声，提高分割的精度，本文把拉普拉斯平滑算法用在了对牙齿三角网格的平滑处理上。该算法首先找到网格G=的每个顶点Vi的1邻域顶点集合，其1邻域定义为与点Vi所有共边的点的合集：

L1(Vi)={Vi|∃edge(Vi,Vj)}

(1)

然后为每个顶点Vi构建1邻域数组Array[n]，n是每个顶点Vi的1邻域顶点的个数。最后迭代所有顶点，对于每个顶点Vi，根据1邻域数组Array[n]里所有顶点坐标的平均值，即采用伞状算子：

(2)

来修改Vi的坐标。重复多次迭代(根据多次实验的结果，迭代次数为20次时，效果最好)。算法的伪代码如下：

网格平滑算法伪代码：

Input: G

//输入网格的几何和拓扑信息

While(count < 20) do

//迭代20次平滑

For i < vertex number of G do

n = L1(Vi)

//每个顶点执行方程(2)

new array[n]

modify coordinate of Viby U(Vi)

//根据方程(3)修改每个顶点的坐标

End for

End while

Return G1

2.3 替换堆栈结构

传统的交互标记控制算法，利用堆栈存储牙齿模型的面片信息，并根据直接插入排序算法对堆栈里的三角面片进行排序。鉴于堆栈元素后进先出的特点，本文采用vector容器代替堆栈存储牙齿模型的面片信息。由于直接插入排序算法比较耗时，本文不再使用排序算法进行排序，充分利用vector容器自身的特点，可以直接在容器内找到曲率最大的元素，这样节省了大量时间。

2.4 vector容器初始化

在每个牙齿上人工拾取一个三角网格并对其进行标记f，找出每个拾取的三角网格f相邻的三个网格fi，把每个网格f的ID，标记m，以及曲率c(f，fi)存入对应的结构体数组中，并放入vector容器中；而曲率c(f，fi)的计算方法如下：

设两个三角网格f1、f2，它们的单位法向量分别为n1和n2，这里设f1、f2的共享边界为AC，AB为其中一个三角面片的非公共边(如图2所示)，则标记的两个三角网格f1,f2，定义式(3)作为它们的相对弯曲程度的函数。

H(f1,f2)=H(f2,f1)=

(3)

图2 相邻两个网格之间相对弯曲程度的计算

2.5 区域生长算法

如果vector容器不为空，每次从vector中取出曲率最大的三角面片fmax，判断fmax的访问标记Visited是否为真，若Visited为真，将fmax直接从vector容器中删除，继续迭代；若Visited为假，设置fmax的访问标记Visited为TRUE，将fmax插入到对应的数组Array[m]中(每个牙齿对应着一个数组，存放访问标记Visited为真的三角面片信息)，把fmax相邻的且访问标记Visited为假的面片fi放入vector容器中，继续迭代。其伪代码如下所示：

区域生长算法的伪代码：

Input:G

//输入网格的几何和拓扑信息

While(vector!=empty)

Select fmax

//选择曲率最大的三角面片

If(Visited==true)

Delete fmaxfrom vector

//删除已经被标记的三角面片

End if

ElseVisited=true

//设置访问标记为真

fmaxPop Array[m]

Neighboring face of fmaxand Visited is false Pop vector

//把相邻的且访问标记为假的面片放入vector中

End while

Return Array[m]

3 算法实验及分析

为验证本文提出的基于网格抽取的牙齿模型快速分割算法的精确性和快速性，利用激光扫描得到的STL文件格式的牙齿模型进行实验，所有实验均在处理器为Intel(R) Core(TM) i5-3470 CPU @ 3.20 GHz，内存4 GB，64位Windows7系统的电脑上运行的。开发软件为Visual Studio 2012和Qt5.5.1，编程语言为C++，所用的工具包为VTK可视化工具包。

图3是两个激光扫描得到的STL文件格式的牙齿网格模型(这两个网格模型均来自天津某医院)。模型一的三角网格数量为365 092，顶点数量为182 544，大小为18 290 KB，体积为49 464.945 m3，面积为11 487.523 m2；模型二的三角网格数量为506 702，顶点数量为253 353，大小为24 742 KB，体积为78 645.768 m3，面积为20 356.143 m2。图4是两个STL文件格式的网格牙齿模型经过交互标记控制算法分割后，不同角度的分割结果图。图5是两个STL文件格式的牙齿网格模型经过本文的算法分割后，不同角度的分割结果图。从图中可以看出，交互标记控制算法分割牙齿存在明显的欠分割(单颗牙齿没有被完全分割出来)现象，特别是在双尖牙处(如图4圆圈围住的牙齿)。而本文的算法分割牙齿，基本上解决了欠分割现象的存在(如图5圆圈围住的牙齿)。其主要原因是交互标记控制算法的区域生长是根据直接插入排序算法从堆栈中每次取曲率最大的三角面片，而双尖牙上的凹区域比较明显，直接插入排序算法则会把双尖牙上凹区域的三角面片排到堆栈的底部，这样区域生长时无法取出该区域的三角面片，所以极容易产生欠分割的现象。而本文的算法利用vector代替堆栈结构，并且没有利用直接插入排序算法，而是利用了vector容器自带的取最大值的方法，这样区域生长时可以非常准确地取出双尖牙凹区域的三角面片，所以本文的算法解决了交互标记控制算法存在的欠分割的现象。

图3 两个牙颌模型

图4 交互标记控制算法分割结果

图5 本文算法分割结果

表1列出了两个牙齿模型在两种算法中的分割时间，可以看出基于网格抽取的牙齿网格快速分割算法的分割时间有了明显的提高，和原来的交互标记控制算法相比，分割时间的加速比提高了6到9倍左右，而且牙齿模型的三角面片数量越少，分割的时间越快。其主要原因是：第一，本文提出的算法包含了网格抽取算法，在牙齿网格分割前，进行了网格抽取，适当地减少了牙齿网格的数量，从而加快了网格的分割速度；第二，交互标记控制算法中包含了一个插入排序算法，而插入排序算法的时间复杂度是O(n2)，所以交互标记控制算法的时间复杂度是O(n2)，而本文提出的算法使用了vector容器自带的寻找最大值的方法，没有使用插入排序算法，因为vector容器自带的寻找最大值的方法的时间复杂度为O(n)，所以本文提出的算法的时间复杂度是O(n)，因此本文提出的算法的分割时间明显快于交互标记控制算法的分割时间。从实验中也可以看出，网格模型的顶点数越多，体积和面积越大，算法执行的速度越慢，但分割的精度越高，对后续操作的影响也越小。

表1 两种算法分割时间比较

4 结 语

本文提出了一种基于网格抽取的牙齿模型快速分割算法。多次实验的结果表明，该算法主要包含三个优点：1) 分割速度快。和传统的交互标记控制算法相比较，本文提出的算法牙齿分割速度有了明显的提高。2) 分割精度高。传统的交互标记算法存在欠分割，本文提出的算法一般不会出现欠分割。3) 适应性强。对于大部分的三维牙齿网格模型，本文提出的算法都可以进行精确的分割。同时，本文的算法除了适用于牙齿网格分割外，对于其他网格模型的分割，也可以用本文提出的算法。当然，本文提出的算法也有一定的局限性，对于一些牙齿网格模型，体积和面积太大，本文的算法不能进行精确的分割，而且对于网格数量超过80万的三维牙齿模型，分割速度也有点慢。

参 考 文 献

[1] Kronfeld T,Brunner D,Brunnett G.Snake-Based Segmentation of Teeth from Virtual Dental Casts[J].Computer-Aided Design and Applications,2010,7(2):221-233.

[2] Kumar Y,Janardan R,Larson B,et al.Improved Segmentation of Teeth in Dental Models[J].Computer-Aided Design and Applications,2013,8(2):211-224.

[3] 李成军,张弛,汪国平.交互标记控制的快速网格分割[J].北京大学学报(自然科学版),2006,42(5):662-667.

[4] 吴松和.牙齿隐形矫正CAD软件开发[D].石家庄:河北科技大学,2013.

[5] Li Z,Wang H.Interactive Tooth Separation from Dental Model Using Segmentation Field[J].Plos One,2016,11(8):e0161159.

[6] Wu K,Chen L,Li J,et al.Special Section on CAD/Graphics 2013:Tooth segmentation on dental meshes using morphologic skeleton[J].Computers & Graphics,2014,38(1):199-211.

[7] Yuan T,Liao W,Dai N,et al.Single-tooth modeling for 3D dental model[J].International Journal of Biomedical Imaging,2010,2010(9):1029-1034.

[8] Mouritsen D A.Automatic segmentation of teeth in digital dental models[D].Birmingham: the university of alabama at birmingham,2013.

[9] Liao S H,Liu S J,Zou B J,et al.Automatic Tooth Segmentation of Dental Mesh Based on Harmonic Fields[J].Biomed Research International,2015,2015(3):1-10.

[10] Mangan A P,Whitaker R T.Partitioning 3D surface meshes using watershed segmentation[J].IEEE Educational Activities Department,1999,5(4):308-321.

[11] Hoffman D D,Singh M.Salience of visual parts[J].Cognition,1997,63(1):29-78.

[12] Koschan A F.Perception-based 3D triangle mesh segmentation using fast marching watersheds[C]//Computer Vision and Pattern Recognition,2003.Proceedings.2003 IEEE Computer Society Conference on.IEEE,2003,2:27-32.

[13] Liu R,Zhang H.Segmentation of 3D Meshes through Spectral Clustering[C]//Computer Graphics and Applications,Pacific Conference.IEEE Computer Society,2004:298-305.

[14] Katz S,Tal A.Hierarchical mesh decomposition using fuzzy clustering and cuts[J].Acm Transactions on Graphics,2003,22(3):954-961.

[15] Schroeder W J.Decimation of triangle meshes[J].Acm Siggraph Computer Graphics,1992,26(2):65-70.