CN105046750B - 一种自动分割全颌牙齿三角网格模型的方法 - Google Patents

Abstract

本发明公开了一种自动分割全颌牙齿三角网格模型的方法。求出每个网格顶点的平均曲率和均方差曲率；获得边界特征点和边界特征区域；将牙颌三角网格模型分离为多个独立的网格区域并进行区分；用不同编号标记牙龈区域和牙齿区域；利用区域增长计算方法处理获得精确分割结果；通过其他牙齿分别到最后一颗磨牙的测地距离对牙齿进行排序并顺序切牙；去除毛刺，删除外翻面片，采用laplacian平滑方法处理后去除狭长三角面片，实现最终牙颌三角网格模型的分割。本发明利用三角网格模型的曲率分布对牙齿模型进行初步分割；利用区域增长法对牙齿模型进行精确的切割；对边界不光滑以及外翻现象进行了光顺化处理；达到快速自动分割牙齿的目的，分割精准、光滑。

Description

一种自动分割全颌牙齿三角网格模型的方法

技术领域

本发明涉及了数字几何处理领域中的网格分割技术，尤其是涉及了一种自动分割全颌牙齿三角网格模型的方法，适用于从牙颌三角网格曲面快速自动分割牙齿三角网格曲面。

背景技术

近年来，数字口腔技术发展迅速，牙齿隐形正畸也逐渐普及。而数字治疗辅助系统要完成的工作就是建立牙齿模型，为矫治方案的设计提供数据支持。要建立精准的牙齿模型，首先需要通过光学方法扫描石膏模型获得牙颌三维数据，再利用数字几何处理技术分离牙齿牙龈数据，并对牙齿缺失部分进行修复，建立可移动的牙齿模型，供医生确定矫治方案。

但由于牙齿的形状和排列因人而异，所以要精确地分离牙齿有一定的困难。而且，目前的一些牙齿分割技术中，很难一次性自动分割所有牙齿，且严重依赖人工交互，分割效率低下。

利用曲面内在属性来识别牙齿边界是一种主流切牙方法，具体包括：1)曲率估计：利用曲面的内在属性，包括主曲率(最小主曲率，平均主曲率等)，对牙齿特征进行识别；2)粗糙定位牙冠区域：利用主成分分析方法(Principal Component Analysis，简称PCA)计算输入模型的局部坐标系，并构建切割平面分离牙齿和牙龈部分，这样的方法只能得到初步的结果；3)阈值法：首先计算整体模型的曲率场，再利用曲率阈值来分割牙齿边界区域。一般来说，用户可以通过交互选择该阈值，也可以通过大量实验来获得较满意的阈值，但由于一个整体阈值难以同时满足所有牙齿分割条件，往往出现过分割或欠分割情况；4)边界优化法：由于牙冠区域存在无效特征区域，且曲率场容易受噪声影响，学者通过形态学操作对边界进行优化处理。但对于复杂牙齿模型来说，形态学操作获得的边界区域存在不完整现象；5)边界定位和优化：学者利用骨架操作抽取牙齿边界区域，然后利用边界优化方法产生平滑精确地边界。但利用该方法抽取的边界存在不闭合现象，因此，需要其他额外的操作以获得平滑精确地牙齿边界。

Kondo等提出一种全自动方法，该方法使用两幅距离图像识别相邻牙齿的分裂边界，将三维空间问题转化为图像问题，但该方法采用用平面去裁剪牙齿数据，对于存在严重错颌情况来说，无法精确定位牙齿牙龈边界。Kronfeld等提出基于snake的方法，该方法给定牙龈上的初始边界区域，然后利用特征吸收场以逐渐收敛获得牙齿边界区域。但模型边界的噪声往往影响特征场，使得该自动方法得到的特征线并不完全符合牙齿边界。其他交互式方法允许用户选择多个边界控制点集，然后利用测地线来连接相邻控制点。该方法直观清晰地分割牙齿模型，但用户需要多次旋转平移模型以确定模型边界点。用户交互过于繁杂，效率低下。Zou等基于调合场交互地分割模型。但该方法需要多次迭代程序以获得全部牙齿的分割。尽管Liao等人利用一致性调和场一次性地分割所有牙齿，但方法的准确性依旧依赖于人工交互。

总之，现有的从牙颌三角网格模型中分割牙齿三角网格模型的方法在应用中都有着各种缺陷，主要表现为精度低、效率底下和交互繁杂。

发明内容

为了解决背景技术中存在的问题，本发明所提供了一种自动分割全颌牙齿三角网格模型的方法，本方法利用数字几何处理中的多项技术来实现快速自动分割边界与边界精确调整两种模式的结合，使得牙齿分割快速、精确且边界光滑。

本发明主要包括六个步骤，流程如下所示：

A.由牙颌三角网格模型中各顶点的曲率值求出每个网格顶点P的平均曲率K_m和均方差曲率K_s；

曲面曲率是微分几何中的重要概念，描述了曲面的局部弯曲程度，曲率值是切割牙齿所需要的重要参数。

B.将满足平均曲率K_m<t₁和均方差曲率K_s>t₂中任一条件的网格顶点P作为边界特征点，t₁为平均曲率阈值，t₂为均方差曲率阈值，所有边界特征点组成边界特征点集合，取所有边界特征点及其相邻的所有三角面片的合集作为边界特征区域；

具体来说，对于一个网格顶点P，若其若其平均曲率K_m<t₁，则令P∈A₁；若其均方差曲率K_s>t₂，则令P∈A₂。由于选取的集合A₁和A₂在牙齿表面都不是完全闭合的曲线，不能把牙齿模型分为不同的独立部分，所以取A₃＝A₁UA₂作为边界特征点的集合。根据边界特征点来计算边界特征区域的三角面片。即对于一个三角面片，只要其任意一个顶点属于集合A₃，那么该三角面片为特征三角面片。

C.根据边界特征区域将牙颌三角网格模型分离为多个独立的网格区域，根据顶点数目将各个独立网格区域区分为牙龈区域、牙齿区域和杂质区域；

D.将分离出来的牙龈区域和牙齿区域作为每个牙齿或牙龈的预分割区域，并进行不同编号标记；

E.根据预分割结果，利用区域增长计算方法处理获取精确分割结果；

F.找出最后一颗磨牙，并通过计算其他牙齿分别到该颗磨牙的测地距离，来依次对其他牙齿进行排序，然后按照排序进行顺序切牙；

G.去除毛刺，删除外翻面片，并采用laplacian平滑方法对牙颌三角模型上的边缘区域进行平滑处理，然后去除狭长三角面片，完成后续边界平滑处理，实现最终牙颌三角网格模型的分割。

所述步骤A中的平均曲率K_m和均方差曲率K_s采用以下公式进行计算：

平均曲率均方差曲率

其中，K₁、K₂为网格顶点P的最大曲率和最小曲率。

所述步骤B中的平均曲率阈值t₁和均方差曲率阈值t₂可采用以下公式进行计算：取所有顶点平均曲率的最大值K_mmax和最小值K_mmin，求出平均曲率阈值取所有顶点均方差曲率的最大值K_smax和最小值K_smin，求出均方差曲率阈值

所述的步骤C具体如下：牙颌三角网格模型中除边界特征区域以外的其他区域被边界特征区域分离为多个独立的网格区域，将特征区域的点和面片做Delete标记，则没有被标记的面片会被分离成独立的网格区域。这些独立的网格区域大多包含一颗牙齿或牙龈的大部分数据。

然后对这些独立网格区域的顶点数目进行统计，并统计独立网格区域的顶点数目，根据独立网格区域的顶点数目确定该独立网格区域的性质，分为牙龈区域、牙齿区域和杂质区域，其中，牙龈区域的顶点数目>牙齿区域的顶点数目>杂质区域的顶点数目，即顶点数目最多的为牙龈区域，顶点数目小于整个牙颌模型总顶点数目的0.5％的独立网格区域为杂质区域，剩余的为牙齿区域。

所述的步骤D具体如下：对于牙龈区域和每个牙齿区域分别标记不同编号，其所包含的顶点全部标记为该顶点所在区域相同的编号，杂质区域与边界特征区域此处不进行标记，一起归为未标记区域。例如牙龈区域编号为R；牙齿区域按顶点数目排列分别编号为Q₁～Q_n，牙龈区域的顶点均标记为R，牙齿区域Q_n的顶点均标记为Q_n。

所述的步骤E中区域增长方法具体如下：

E.1)将所有杂质区域与边界特征区域的顶点组成集合U，为杂质区域与边界特征区域的每个顶点t建立一个三元组<t，m，d>，其中m为其标记，m的初始值为U，d为该顶点到相邻顶点的特征距离；

E.2)对于牙龈区域和牙齿区域中的每一个顶点，将其1-邻域中属于杂质区域与边界特征区域的顶点的三元组归属到集合Q′中，并按照特征距离d从小到大排列；

E.3)将集合Q′中特征距离d最小的顶点的标记改为与该顶点相邻的牙龈区域或牙齿区域的标记，并将特征距离d最小的顶点1-邻域中属于杂质区域与边界特征区域的顶点再加入到集合Q′中，按照特征距离d从小到大重新排列；

E.4)重复上述步骤E.3)，直至牙颌三角网格模型中的所有顶点均被标记，停止增长，得到精确分割结果。

所述步骤E.3)和E.3)中，对于任一网格顶点p到相邻另一网格顶点t的特征距离具体采用以下方式计算：

d(p,t)＝-K_m

其中，K_m是网格顶点t的平均曲率。

所述步骤F中最后一颗磨牙采用以下方式进行寻找：根据步骤E得到的分割结果(分离出来的牙齿面片集合)，分别计算整个牙颌三角网格模型的重心坐标r_0和每一颗牙齿的中心坐标{r_1,…,r_n}，中心坐标离重心坐标r_0的欧式距离最远所对应的牙齿为最后一颗磨牙。

所述步骤F中的测地距离是指牙颌三角网格模型上的两顶点沿着模型表面路径上距离。

虽然本方法能高效自动地分离牙齿和牙龈数据，但无法保证获得牙齿的边界线是平滑的。所以还需经所述步骤G对模型进行平滑修复处理。

由于分割出的边界会有很多“毛刺”，不利于后面的平滑处理，所以首先要对模型进行去毛刺处理。所述的毛刺是指三个顶点均为边界点的三角面片，若一个三角面片的三个顶点都为边界点，则可将该三角面片作为毛刺删除。

由于牙齿模型本身的缺陷，导致获得的边界区域存在“外翻”现象。“外翻”的边界会导致比较大的误差，往往影响后期的模型修复操作。因此，需要删除“外翻”面片，并对牙齿边界线进行平滑操作，平滑操作采用牙颌三角模型上的边界区域laplacian平滑和去除狭长三角面片相结合的方法。

由于laplacian smooth操作不改变带状边界的拓扑，但改变带状边界顶点的坐标，因此，会出现狭长三角形，影响后期的模型修复工作。通过检测并删除狭长三角形集合，获得边界平滑且形状基本逼近原始模型的牙齿三角网格模型。

所述步骤G中删除外翻面片具体采用以下方式：选取牙颌三角网格模型最外围的一层带状边界区域，边界区域内法向变化≥10°的区域作为外翻面片进行删除。

所述步骤G中狭长三角形具体为任意一个角≥160°的三角形。

本发明的有益效果是：

本发明方法可快速高效精确地自牙颌中分离出单个的牙齿模型，提高了精确度和工作效率，避免了复杂的人机交互过程，并且使分离出的模型边界光滑，便于后期处理。

附图说明

图1为实施例的牙颌三角网格模型。

图2为实施例的牙颌三角网格顶点平均曲率值的分布情况。

图3为实施例的牙颌三角网格顶点均方差曲率值的分布情况。

图4为实施例的牙颌三角网格中的特征区域分布情况。

图5为实施例的分割出来的单独的牙齿模型，存在一定的外翻现象。

图6为实施例的毛刺现象网格示意图。

图7为实施例的去毛刺之后的牙齿模型。

图8为实施例的删除外翻面片并经过laplacian smooth操作之后的牙齿模型。

图9为实施例的经过lapliacian smooth操作之后产生的狭长三角面片示意图。

图10为实施例的三角面片中区域增长法示意图之一。

图11为实施例的三角面片中区域增长法示意图之一。

具体实施方式

下面结合附图及具体实施例对本发明作进一步详细说明。本发明的实施是为了对本发明进一步说明，而非对本发明的发明范围的限制。

本发明方法的实施例及其具体过程如下：

导入牙颌三角网格模型，如图1所示，该模型包共含97502个顶点和195000个三角面片。

1)由牙颌三角网格模型中各顶点的曲率值求出每个网格顶点P的平均曲率K_m和均方差曲率K_s；平均曲率分布如图2所示，均方差曲率分布如图3所示。其中，色彩灰度反映了曲率值大小，颜色灰度越高，表明曲率值越小。

2)取所有顶点平均曲率的最大值K_mmax和最小值K_mmin，求出平均曲率阈值t₁和均方差曲率阈值t₂，本例中，t₁＝0.20t₂＝0.90。

对于一个网格顶点P，若其若其平均曲率K_m<t₁，则令P∈A₁；若其均方差曲率K_s<t₂，则令P∈A₂，取A₃＝A₁UA₂作为边界特征点的集合。取所有边界特征点及其相邻的所有三角面片的合集作为边界特征区域，如图4即为取并集之后的边界特征区域分布。

3)牙颌三角网格模型中除边界特征区域以外的其他区域被边界特征区域分离为多个独立的网格区域，然后对这些独立网格区域的顶点数目进行统计，根据独立网格区域的顶点数目确定该独立网格区域的性质：顶点数目最多的独立网格区域为牙龈区域；顶点数目小于整个牙颌模型总顶点数目的0.5％的独立网格区域为杂质区域在本例中取值500；其余独立网格区域为牙齿区域。

4)将分离出来的牙龈区域和牙齿区域作为每个牙齿或牙龈的预分割区域，并进行不同编号标记；例如，牙龈区域编号为R；牙齿区域按顶点数目排列分别编号为Q₁～Q_n，牙龈区域的顶点均标记为R，牙齿区域Q_n的顶点均标记为Q_n。

5)根据预分割结果，利用区域增长计算方法处理获取精确分割结果；

E.1)将所有杂质区域与边界特征区域的顶点组成集合U，为杂质区域与边界特征区域的每个顶点t建立一个三元组<t，m，d>，其中m为其标记，m的初始值为U，d为该顶点到相邻顶点的特征距离；

E.2)对于牙龈区域和牙齿区域中的每一个顶点，将其1-邻域中属于杂质区域与边界特征区域的顶点的三元组归属到集合Q′中，并按照特征距离d从小到大排列；

E.3)将集合Q′中特征距离d最小的顶点的标记改为与该顶点相邻的牙龈区域或牙齿区域的标记，并将特征距离d最小的顶点1-邻域中属于杂质区域与边界特征区域的顶点再加入到集合Q′中，按照特征距离d从小到大重新排列；

E.4)重复上述步骤E.3)，直至牙颌三角网格模型中的所有顶点均被标记，停止增长，得到精确分割结果。

如图10为区域增长法示意图。图中，左边三个实心点被标记为Q₁，右边四个空心点被标记为Q₂，其余点标记为U；将与标记Q₁、Q₂相邻的点，即图中的t₁～t₇均加入集合Q′，并将其按照d值从小到大排列为：

<t₃，U，-3.81>

<t₂，U，-3.52>

<t₇，U，-2.76>

<t₄，U，-1.01>

<t₁，U，0.19>

<t₆，U，0.27>

<t₅，U，-0.39>

如图11所示，将t₃标记为Q₂；从集合Q′中删除点t₃并将与之相邻的点t₈加入集合Q′中，重新按照d值排列，重复上述步骤。

6)找出最后一颗磨牙，并通过计算其他牙齿分别到该颗磨牙的测地距离，来依次对其他牙齿进行排序，然后按照排序进行顺序切牙；将每颗牙齿分离为独立的图层，如图5为分离之后的牙齿图层。

7)对切割好的牙齿进行去毛刺处理，若一个三角面片的三个顶点都为边界点，则可将其作为为毛刺删除，如图6所示，标记为Tr₁和Tr₂的两个三角形即为毛刺。经删除毛刺处理的牙齿模型如图7所示。

8)由于牙齿模型数据存在一定的缺失，所以切割好的牙齿边缘会出现“外翻”现象，如图5所示。对于存在外翻现象的牙齿模型，选择边缘的一层三角面片，对法向变化较大的点进行删除操作，如此重复两次，可基本消除外翻现象。在本例中，法向变化的判别标准是与上一层相比，变化超出10°。

9)利用模型边界选择面片操作和面片扩展操作，选择模型边界的具有厚度为2的带状边界，利用一次laplacian smooth操作，对选中的带状边界进行平滑处理，处理后的模型如图8所示，边缘光滑。

10)检测并删除狭长三角形集合，本例中取任意一个角≥160°的三角形为狭长三角形，经过lapliacian smooth操作之后产生的狭长三角面片如图9所示，最终获得边界平滑且形状基本逼近原始模型的牙齿三角网格模型。

由此，本发明利用三角网格模型的曲率分布对牙齿模型进行初步分割，并利用边界增长法，对牙齿模型进行精确的切割。本发明完成所有牙齿分割的处理时间是4.3s。相比传统方法每颗牙齿单独处理，传统方法需要灵活的手动点击操作处理，对于单个牙齿即使操作快的情况也需要几十分钟时间，因此效率大大提高，具有突出显著的技术效果。

Claims (9)

1.一种自动分割全颌牙齿三角网格模型的方法，其特征在于该方法包含下列步骤：

A.由牙颌三角网格模型中各顶点的曲率值求出每个网格顶点P的平均曲率K_m和均方差曲率K_s；

B.对于一个网格顶点P，提取边界特征点，取所有边界特征点及其相邻的所有三角面片的合集作为边界特征区域；

C.根据边界特征区域将牙颌三角网格模型分离为多个独立的网格区域，根据顶点数目将各个独立网格区域区分为牙龈区域、牙齿区域和杂质区域；

D.将分离出来的牙龈区域和牙齿区域作为每个牙齿或牙龈的预分割区域，并进行不同编号标记；

E.根据预分割结果，利用区域增长计算方法处理获取精确分割结果；

F.找出最后一颗磨牙，并通过计算其他牙齿分别到该颗磨牙的测地距离，来依次对其他牙齿进行排序，然后按照排序进行顺序切牙；

G.去除毛刺，删除外翻面片，先采用laplacian平滑方法对牙颌三角模型上的边缘区域进行平滑处理，然后去除狭长三角面片，完成后续边界平滑处理，实现最终牙颌三角网格模型的分割。

2.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述步骤B中，将满足平均曲率K_m<t₁和均方差曲率K_s>t₂中任一条件的网格顶点P作为边界特征点，t₁为平均曲率阈值，t₂为均方差曲率阈值。

3.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述步骤A中的平均曲率K_m和均方差曲率K_s采用以下公式进行计算：

平均曲率均方差曲率

其中，K₁、K₂为网格顶点P的最大曲率和最小曲率。

4.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述的步骤C具体如下：牙颌三角网格模型中除边界特征区域以外的其他区域被边界特征区域分离为多个独立的网格区域，并统计独立网格区域的顶点数目，根据独立网格区域的顶点数目确定该独立网格区域的性质，分为牙龈区域、牙齿区域和杂质区域，牙龈区域的顶点数目>牙齿区域的顶点数目>杂质区域的顶点数目。

5.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述的步骤D具体如下：对于牙龈区域和每个牙齿区域分别标记不同编号，其所包含的顶点全部标记为该顶点所在区域相同的编号，杂质区域与边界特征区域此处不进行标记，一起归为未标记区域。

6.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述的步骤E中区域增长方法具体如下：

E.1)将所有杂质区域与边界特征区域的顶点组成集合U，为杂质区域与边界特征区域的每个顶点t建立一个三元组<t，m，d>，其中m为其标记，m的初始值为u，d为该顶点到相邻顶点的特征距离；

E.2)对于牙龈区域和牙齿区域中的每一个顶点，将其1-邻域中属于杂质区域与边界特征区域的顶点的三元组归属到集合Q′中，并按照特征距离d从小到大排列；

E.3)将集合Q′中特征距离d最小的顶点的标记改为与该顶点相邻的牙龈区域或牙齿区域的标记，并将特征距离d最小的顶点1-邻域中属于杂质区域与边界特征区域的顶点的三元组再加入到集合Q′中，按照特征距离d从小到大重新排列；

E.4)重复上述步骤E.3)，直至牙颌三角网格模型中的所有顶点均被标记，停止增长，得到精确分割结果。

7.根据权利要求6所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述步骤E.2)和E.3)中，对于任一网格顶点p到相邻另一网格顶点t的特征距离具体采用以下方式计算：

d(p,t)＝-K_m

其中，K_m是网格顶点t的平均曲率。

8.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述步骤F中最后一颗磨牙采用以下方式进行寻找：根据步骤E得到的分割结果，分别计算整个牙颌三角网格模型的重心坐标r_0和每一颗牙齿的中心坐标{r_1,…,r_n}，中心坐标离重心坐标r_0欧式距离最远所对应的牙齿为最后一颗磨牙。

9.根据权利要求1所述的一种自动分割全颌牙齿三角网格模型的方法，其特征在于：所述步骤F中的测地距离是指牙颌三角网格模型上的两顶点沿着模型表面路径上距离。