基于形态学方法的牙颌模型自动分割算法

陈 虎 马珂楠 胡毕富

1(北京大学口腔医院口腔医学数字化研究中心 北京 100081)

2(北京航空航天大学机械工程及自动化学院飞行器制造工程系 北京 100191)

0 引 言

在口腔医学领域中,CAD/CAM技术的引入正逐渐改变传统口腔治疗方式,而牙科CAD系统作为口腔数字化医疗中的重要环节,对其设计策略和算法的研究都具有重要的意义[1]。口腔医学包括口腔修复、正畸、种植以及颌面修复等主要方面,每种不同治疗方式需要构建的CAD系统差异性很大,但很多情况下都需要从牙颌的数字化模型分离开牙齿和牙龈部分,实现对牙颌模型的分割,例如口腔正畸中的虚拟排牙等。受不同牙颌的形状差异性,以及牙齿畸形、牙列缺损等因素影响,目前的牙科CAD软件大多是以繁琐、低效的手动交互方式完成分割的,对牙颌模型分割的自动化程度较低。

牙颌模型的分割属于数字几何处理领域中的三维网格分割问题,近年来三维网格分割已成为一个研究热点[2],但已有的网格模型分割方法不能实现对任意形状的模型的自动分割。目前针对具有复杂形状的牙颌模型,已提出了较多的分割方法。文献[3-4]从全景图像中构建了两幅距离图像实现了牙齿的自动分割,避免了直接处理复杂的三维网格数据问题,但也使得这种分割方式不适合于处理三维牙颌模型的牙科CAD系统。目前比较常见的方法是通过分析牙颌模型的离散曲率实现对牙齿和牙龈的分离。Kronfeld等[5]利用牙齿和牙龈边界附近顶点的曲率值较大的特点,构建了活动轮廓线,之后通过能量法局部优化边界线实现牙颌分割。但该方法不易处理活动轮廓线出现不连续的牙颌模型。Kumar等[6]通过设置曲率阈值选取出高曲率值的顶点获得龈缘线轮廓,在分离开牙龈和牙齿后,将龈缘线上的尖点作为相邻牙的分割边界点分离开相邻牙齿。该方法也难以处理牙齿畸形等情形,比较依赖于标准理想的牙颌模型。为了提高分割的质量,在曲率获取龈缘线的基础上,利用形态学操作获取更精细的牙齿边界线。Zhao等[7]在加入较多交互的方式下,利用形态学操作实现了牙颌模型的分割。Yuan等[8]通过形态学操作分离开牙齿和牙龈曲面,并实现了对相邻牙齿的分离以及对牙齿曲面的修补,该方法同样需要较多的手动交互操作。Wu等[9]在形态学操作分割牙颌模型的基础上,实现了少量交互下对牙齿牙龈以及相邻牙齿的分割。除了利用牙颌网格模型的曲率特性实现分割外,还有其他反映网格曲面特征的信息可用于牙颌模型的分割。Zou等[10]通过计算网格上的调和场来反映网格的几何细节,再手动选取具有约束作用的特征点,实现牙颌模型的分割,在之后的改进方法中[11],通过对牙颌模型的多层切分实现了自动选择特征点,从而完成对牙颌模型的自动分割,但该方法在处理相邻牙齿贴靠紧密的情况时难以将牙齿分割开。Hwang等[12]也基于调和场进行牙颌模型的分割,并采用狄利克雷边界条件,凸分割等后处理方法,在具有噪声和网格碎片的临床实际的牙颌模型上取得了一定的效果。Fan等[13]采用离散Morse理论分析曲面结构,把分割问题转化成从曲面上冗余特征线族中提取分割边界线,再通过优化方法从底层分割中聚合产生子部分,获得了较好的分割效果。Kim等[14]通过对牙颌模型进行截面投影,随后应用区域增长算法实现牙齿的自动分割。He等[15]则基于曲率信息和局部分布密度信息初步划分牙齿-牙龈边界和牙齿间边界,然后通过区域增长算法实现牙齿分割。Pupykina等[16]利用可视化工具包进行交互,基于用户选择的少量关键特征点,通过模糊综合评判的方法确定牙齿的边界。近年来也出现了许多基于机器学习的方法。Xu等[17]建立了牙齿分割的深度卷积神经网络,将牙颌网格面片的特征作为输入,经过两层神经网络,第一层用于区分牙齿和牙龈区域,第二层用于区分不同牙齿区域,最终通过模糊聚类对结果进行优化,取得了较好的分割效果。Lian等[18]建立了MeshSegNet,实现了牙颌模型多尺度的特征提取和学习,采用图约束(Graph-Constrained)学习模块进行局部特征提取,并通过密集融合(dense fusion)策略进行更高层次的特征学习,最后采用图割算法对分类结果进行优化调整。Sun等[19]基于图卷积神经网络实现了牙齿和牙龈的同时分割。Zhang等[20]将牙颌三维模型映射到二维的调和参数空间,转换为图像信息,从而利用卷积神经网络进行牙齿分割,随后将二维图像上的分割结果映射回三维实体模型中,并通过改进的模糊聚类和分割方法完成边界的调整和优化。这些基于机器学习的方法虽然在一些数据集上取得了较高的准确率,但是受训练数据规模、质量等影响较大,在处理数据质量较差的临床实际模型时可能无法达到理想的效果。另一方面,由于训练数据通常需要由专业人员进行大规模和高质量的标注,成本高昂,训练和测试时对机器运算能力也有非常高的要求。

现有的牙科CAD系统在处理牙颌模型的分割时,都是以手动交互的方式分出牙齿再进行其他处理,频繁地交互操作难以提高CAD系统的使用效果。但受到牙颌的形状差异性的影响,以及牙齿排列畸形、牙列缺损等口腔病症的出现,使得牙颌模型的自动分割一直没有好的解决方案,还没有通用的方法,已有方法都或多或少存在交互和准确性问题。

本文在现有形态学方法处理牙颌模型分割的基础上[21],提出了一种自动闭合牙颌模型轮廓特征线的方法。牙颌三角网格模型在经过形态学处理得到单个顶点宽度的特征线后,有时存在不连续的情形。首先通过计算特征顶点的复杂度[8,21]判断出特征线的端点,之后通过距离约束需要配对连接的端点以及无效的端点,再逐层分别搜索端点周围的非特征顶点并标记,直到两配对的端点都搜索到同一非特征顶点,以该顶点作为两个配对端点的中点,根据最短路径原则分别与两端点相连,从而使得特征线封闭。该方法解决了原来需要手动处理特征线不封闭的问题,实现了对牙颌模型的自动分割。

1 网格曲面上的形态学操作

1.1 三角网格的描述

对于三角网格模型M通常可表示为M={V,T},其中V是网格顶点的集合,T是网格上三角面片面的集合。

(1)

三角网格拓扑信息的建立主要考虑网格M的顶点V、边E和三角面片T的信息以及相互之间的遍历和访问。

牙颌模型曲面的特征提取主要可归结为三角网格分割问题,特征线一般是分割的边界线。实现对牙颌模型的分割关键在于龈缘线的提取,在牙颌模型中,龈缘线及相邻牙齿之间的特征线附近的顶点平均曲率值比其他区域偏大,可以通过平均曲率值初步区分出龈缘线附近的特征区域。为了计算网格顶点上的离散曲率,对于每个顶点v∈V的主曲率可分别表示为κ1(v)和κ2(v),从而平均曲率表示为：

(2)

在牙颌模型的特征线附近都是凹下去的“谷线”区域,这些区域顶点的平均曲率值都是偏小的负值,通常直接设置曲率阈值ε来获得初始的特征点集F：

F={v∈V|κH(v)≤ε}

(3)

根据曲率阈值初步筛选出的特征点集有较多的噪声点,而且能够表示龈缘线和相邻牙齿分割线的特征点也存在缺陷,需要做后续处理,才能获得适用于牙颌模型分割的只有一个顶点宽度的特征线。

1.2 形态学操作

形态学方法可以处理初始特征点集并获得只有一个顶点宽度可用于模型分割的特征线。Rössl等[21]最早将形态学操作应用于三角网格曲面的特征线提取,并将用于图像处理的形态学操作转化到处理三角网格曲面特征区域上,其中使用的几种算子主要包括：膨胀、腐蚀、开闭操作、骨架化以及裁剪操作等。该方法主要在于利用顶点的邻接关系来表示形态学中的结构元,为此定义了任意顶点vi的邻接关系nhd,对于顶点vi的1-ring邻接关系定义为：

nhd{i}={i}∪{j|∃edge(vi,vj)}

(4)

邻接关系可以逐次扩展为n-ring邻接关系nhdn：

nhdn+1=nhd(nhdn{i})n>1

(5)

1.2.1膨胀和腐蚀操作

膨胀操作可以将比较接近但不连续的特征区域合并在一起,并且特征区域向周围扩展,比原来的特征区域连续性更好。对特征区域F以nhdn为结构元的膨胀可表达为：

dilaten(F)={j|∃i∈F：j∈nhdn{i}}

(6)

膨胀操作是把特征区域F中顶点v的n-ring邻接点全部都加入到特征集F中。

腐蚀操作的作用和膨胀操作相反,可以将特征区域细化,达到简化区域,使得特征结构更清晰的目的。对特征区域F以nhdn为结构元的腐蚀可表达为：

eroden(F)={j|nhdn{j}⊆F}

(7)

膨胀和腐蚀操作通常是组合来应用的,网格上的特征区域需要经过一系列的膨胀或腐蚀操作,两者不同的组合会出现不同的操作效果,其中包括常用的开操作、闭操作等,这些操作会使得特征区域在保持原有基本形状,获得较清晰的骨架特征。

1.2.2骨架化和裁剪操作

骨架化操作可以把网格上的特征区域细化为单个顶点宽度的特征线,骨架化操作的实质是在判断应去除特征区域中的哪些顶点。特征区域中的顶点根据其复杂度的不同可以分为三类[13]。对于顶点vi的复杂度CP(i)定义为：

(8)

式中：Cnei(i)表示顶点vi依次按顺时针的1-ring邻接顶点{v1,v2,…,vm},如果vk∈Cnei(i)且vk∈F,就令Cnei(i)k=1,否则Cnei(i)k=0。当CP(i)≥4时,顶点vi就是复杂点(Complex Vertex),构成的点集记作C;当CP(i)=0时是中心点(Center Vertex),构成的点集记作⊙;其他的特征点为外侧点(Disk Vertex),其点集记作○,通常CP(i)=2。

骨架化操作就是去除外侧点并保留复杂点的过程,可用数学式表示为：

(9)

骨架化操作产生短的无关的分支可以通过裁剪操作消除。令S⊆F表示获得的骨架的特征点集合,裁剪操作为：

(10)

2 牙颌模型自动分割算法

2.1 算法概述

本文的牙颌模型分割算法中使用的牙颌模型是通过三维扫描仪根据牙颌的石膏模型获取的点云数据生成的,牙颌数字化模型是一张三角网格曲面。本文的算法首先使用形态学方法生成龈缘线以及相邻牙交线这两种能用于牙齿分割的特征线,再通过端点的距离约束用最短路径原则连接特征线断开的地方,得到封闭特征线用于牙颌模型的分割[22]。如图1所示,算法的主要步骤有：

图1 算法概述

Step1初始化特征点集。计算每个顶点v的平均曲率值κH(v),对于给定的曲阈值ε,若满足κH(v)≤ε,则把该点v记为特征点加入到集合F。

Step2删除噪点。将集合F中特征点按是否相连分类得到m个特征区域{F1,F2,…,Fm},记num(Fi)为特征区域Fi中顶点个数,删除num(Fi)/F<εd的特征区域,得到新的特征点集F,εd是设置的删除比例。

Step3形态学处理。对特征点集依次做膨胀F=dilate(F),腐蚀F=erode(F),骨架化F=skeletonize(F),修剪F=prune(F)四种操作,得到单个顶点宽度的特征线F。

Step4识别特征线上的端点。对于顶点vi∈F,若vi的复杂度CP(vi)≤2,则把vi放入连接端点的集合ConnectP。

Step5对集合ConnectP中的点配对,计算连接中点MidP,根据最短路径连接配对的端点,得到新的特征点集F。

Step6根据封闭的特征线F分割牙颌模型。

下面具体介绍Step5中的方法。

2.2 自动闭合特征线

特征线的端点集合ConnectP中的顶点既有用于连接间断处的端点,也有一些未被完全删除的分支上的末端点,会造成干扰,如图1(h)所示。需要解决如何从ConnectP中找出需要相互连接的顶点,可以通过设置一个距离范围δmax来约束端点的连接,避免错误连接的出现。

令δ(v)表示顶点u到v的最短距离：

(11)

式中：Pmin(u,v)表示顶点u到v的最短路径。查找相互连接的端点的过程是不断计算每个顶点的vi∈ConnectP与其n-ring邻域内顶点u在δmax范围内的最短距离δ(v),对n依次取n=1,2,…,n,每计算一次都要判断顶点u是否在其他端点vj也被访问,如果访问过,则把顶点u存入连接中点的集合MidP中,并将顶点u与端点vi、vj根据最短路径相连,这个过程如图1(i)、(j)所示。对ConnectP中的端点同时计算,查找中间点,如果是特征线上的分支端,通常在δmax范围内不会找到与之配对的端点,在特殊情况下如果存在,可以在形成封闭轮廓后,对网格模型做分割时剔除这些小的区域,不影响分割效果。

3 实验结果

本研究中用于调参的牙颌模型如图2所示,这些模型是由激光扫描仪对牙颌的石膏模型扫描获取的曲面信息生成的三角网格曲面,图2(a)、(b)、(c)中模型的顶点数量分别为1.53×105、1.32×105、1.56×105。程序中的平均曲率阈值取ε=0.66,在删除噪声点过程中的删除比例取εd=0.01,在闭合特征线过程中的距离范围取δmax=3.5 mm。牙齿与牙龈的分割以及相邻牙齿之间的分割都比较准确完整,表明自动的牙颌模型分割算法具有良好的分割效果。

图2 本文算法自动分割效果

本研究收集了另外5种常见模型(图3)共49颗牙齿进行测试,本算法成功实现了46颗牙齿的自动分割,成功率达到93.88%。

图3 本文算法在5种常见模型上的测试

4 结 语

本文改进了基于形态学方法的牙颌模型分割算法,解决了依靠交互方式选取牙颌模型特征线端点的问题,实现了牙颌模型的自动分割。实验表明,该算法的分割效果较好,适用于不同牙颌模型的自动分割。

相关的后续工作包括：研究其他特殊的牙病患者模型对算法的影响,进一步提高算法的适用性;建立对牙颌模型分割效果的评估方法,判断分割后的模型能否用于患者的治疗。