CN109920058B - 一种基于各向异性度量的牙齿分割方法 - Google Patents

Abstract

本发明提供了一种基于各向异性度量的牙齿分割方法，属于医疗辅助器材设计领域。该方法包括：(1)读取牙齿模型：通过3D扫描获得牙齿模型，所述牙齿模型是三角形网格；(2)计算三角形网格中各个顶点的各向异性度量张量；(3)利用鼠标选择种子点；(4)计算三角形网格中各个顶点与所述种子点之间的基于各向异性度量的最短距离；(5)根据各个顶点与所述种子点之间的基于各向异性度量的最短距离实时获得鼠标当前位置与种子点之间的路径；(6)对步骤(5)得到的路径进行平滑得到平滑的路径；(7)利用鼠标点击确认路径；(8)判断路径是否封闭，如果是，则转入步骤(9)，如果否，则返回步骤(4)；(9)利用所述平滑的路径分割牙齿。

Description

一种基于各向异性度量的牙齿分割方法

技术领域

本发明属于医疗辅助器材设计领域，具体涉及一种基于各向异性度量的牙齿分割方法。

背景技术

随着工业4.0时代的到来，三维扫描技术，3D打印技术和计算机辅助设计技术得到了迅速发展。这使得获取牙齿的三维数字模型和制作个性化口腔种植和正畸医疗辅助器材成为可能。这些年基于数字化技术的牙齿种植和正畸设计系统发展迅速。

典型的牙齿种植和正畸过程是：(1)通过传统印模然后3D扫描或者直接口内3D扫描获取牙齿数字模型，获取患者的头部CT；(2)牙齿分割和CT组织分割和重建；(3)个性化种植牙导板和正畸矫治器设计；(4)进行工程分析和生产制造；(5)实施治疗。而牙齿的分割是后续工作(步骤2,3,4)的基础，对数字化口腔系统十分重要。牙齿模型的质量参差不齐，这是由牙齿模型获取方法(三维扫描技术)决定的，因此设计一个方便实用，并且可以通过交互控制分割质量的工具显得十分重要。

针对牙齿模型的分割算法可以分为自动化的和半自动化的。自动化的方法不需要用户交互，使用非常方便；半自动化的方法仅仅需要一些简单的交互，程序可以根据这些交互自动地进行牙齿分割。然而这些方法由于缺少对分割效果的控制，鲁棒性较差，即在模型质量较高时可以获得较好的分割效果，然而在模型质量较低时分割效果较差，甚至出现无法分割的现象，并且分割精度较差。

为了给予用户足够的分割控制权，Zhuang等人(请参考Zhuang Y,Zou M,Carr N,et al.Anisotropic geodesics for live-wire mesh segmentation[J].ComputerGraphics Forum.2014,33(7):111-120.)通过用户输入的两个端点定义一条线作为分割边界，该分割边界是曲面上基于各向异性度量的最短路径。他的主要贡献是设计一个度量，使路径更接近突出和凹陷的特征，或更好地与曲率方向对齐。该工具为用户提供了对结果的充分控制，同时最大限度地减少了他们的工作量，但是，该方法需要在初始化时进行网格嵌入和局部细分，并且计算分割线时间复杂度较高，由于它是实时交互式工具，因此对时间复杂度高是无法容忍的。此外，该工具有时还无法进行牙齿分割。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于各向异性度量的牙齿分割方法，方便了分割边界的绘制，时间复杂度低，并且可以控制分割质量，能够实时交互，获得平滑的分割边界并针对牙齿分割进行优化。

本发明是通过以下技术方案实现的：

一种基于各向异性度量的牙齿分割方法，所述方法包括：

(1)读取牙齿模型：通过3D扫描获得牙齿模型，所述牙齿模型是三角形网格；

(2)计算三角形网格中各个顶点的各向异性度量张量；

(3)利用鼠标选择种子点；

(4)计算三角形网格中各个顶点与所述种子点之间的基于各向异性度量的最短距离；

(5)根据各个顶点与所述种子点之间的基于各向异性度量的最短距离实时获得鼠标当前位置与种子点之间的路径；

(6)对步骤(5)得到的路径进行平滑得到平滑的路径；

(7)利用鼠标点击确认路径；

(8)判断路径是否封闭，如果是，则转入步骤(9)，如果否，则返回步骤(4)；

(9)利用所述平滑的路径分割牙齿。

所述步骤(2)的操作包括：

(21)计算三角形网格中的每个顶点处的最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂和法向量；

(22)根据最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂获得三角形网格中的每个顶点的各向异性度量张量。

所述步骤(22)的操作包括：

三角形网格中的每个顶点对应一个各向异性度量张量，每个所述各向异性度量张量的特征向量是e₁和e₂，每个所述各向异性度量张量的特征值是λ₁和λ₂；

如果进行齿龈分割，则使用公式(3)计算得到特征值λ₁和λ₂：

如果进行牙齿咬合面分割，则使用公式(4)计算得到特征值λ₁和λ₂：

所述步骤(3)的操作包括：

利用鼠标在分割边缘处选择一个点作为种子点，具体如下：

如果进行齿龈分割，则在齿龈交界处选择一个点作为种子点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点作为种子点。

所述步骤(4)的操作包括：

采用half-edge结构存储所述三角形网格，所述三角形网格中的每个顶点包含一个三元组(dist,final,pred)，其中dist表示该顶点与所述种子点之间的基于各向异性度量的最短距离，final是一个Boolean类型的标记，用于防止重复计算，pred表示在一条路径中该顶点的前一个顶点；

求取各个顶点的三元组(dist,final,pred)的步骤如下：

(41)将种子点的dist字段设置为0，将其他顶点的dist字段设置为无穷大；

(42)创建一个优先队列Q，Q按照各个顶点的dist字段进行排序；

(43)将种子点加到Q中；

(44)当Q不为空时进行如下操作：

(441)从Q中取出一个顶点v，将该顶点v的final字段设置为true；

(442)对于所有与顶点v相邻的顶点w，如果w的final字段为false，则进行如下操作：如果w.dist>v.dist+l_g(v,w)，则w.dist＝v.dist+l_g(v,w),w.pred＝v；并且，如果w不在Q中，则将w加入到Q中；

其中，l_g(v,w)表示两个顶点v、w之间的距离，两个顶点v、w对应的各向异性度量张量分别是g_v,g_w，采用公式(5)计算v、w两点间的距离l_g(v,w)：

其中，

均采用公式(2)计算得到：

其中，x表示从顶点b到顶点a的向量，θ是x和e₁之间的夹角。

所述步骤(5)的操作包括：

找到与鼠标当前位置最近的顶点，从该顶点开始，通过各个顶点的pred字段遍历各个顶点，直到找到的顶点是种子点为止，这样就获得了鼠标当前位置与种子点之间的路径。

所述步骤(6)的操作包括：

(61)将步骤(5)得到的路径记作曲线C，将曲线C以及曲线C经过的边的k邻域分割出来，记作k_ring_mesh；

(62)将k_ring_mesh通过保形映射的方法展开到平面上,记作k_ring_mesh_2D,这样，曲线C变成了2D曲线C^2D；

(63)在平面上使用中点法对2D曲线C^2D进行平滑处理，平滑强度记为s，平滑处理后的2D曲线为S_C^2D；

(64)将S_C^2D映射回原始的三角形网格上，成为曲线S_C，即平滑后的路径。

所述步骤(64)的操作包括：

计算S_C^2D与k_ring_mesh_2D每条边的交点，如果S_C^2D的一条边与k_ring_mesh_2D相交的边的两个顶点为index1,index2，交点为pt_int^2D，则将交点信息记在一个三元组(index1,index2,α)中，其中α表示交点与两个顶点

和

的位置关系，利用公式(6)计算得到α：

pt_int^2D＝(1-α)·pt^2D(index1)+α·pt^2D(index2) (6)

然后根据公式(7)计算出pt_int^3D：

pt_int^3D＝(1-α)·pt^3D(index1)+α·pt^3D(index2) (7)

其中，pt^2D,pt^3D是k_ring_mes_2D和k_ring_mesh的顶点；

所有的pt_int^3D构成了曲线S_C。

所述步骤(7)的操作包括：

利用鼠标在分割边缘处选择一个点，并把该点作为种子点；

利用鼠标在分割边缘处选择一个点的操作包括：

如果进行齿龈分割，则在齿龈交界处选择一个点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点。

所述步骤(8)中的判断路径是否封闭的操作包括：

判断步骤(7)得到的种子点与第一个种子点是否重合，如果重合，则判定路径封闭，如果不重合，则判定路径不封闭。

与现有技术相比，本发明的有益效果是：利用本发明能够进行实时交互，并且需要更少的交互，方便了分割边界的绘制，实现了对分割质量的控制，使获得的分割边界更加精确，而且获得了平滑的分割边界，且大大降低了时间复杂度。

附图说明

图1本发明方法的步骤框图；

图2曲面上的度量张量的示意图；

图3(a)|k_max|-|k_min|的可视化结果；

图3(b)|k_max|+|k_min|.的可视化结果；

图4(a)曲面的最大和最小曲率方向，深色短线的方向为最大曲率的方向，浅色短线的方向为最小曲率的方向。深色的长曲线为需要获取的齿龈分界线；

图4(b)曲面的最大和最小曲率方向，深色短线的方向为最大曲率的方向，浅色短线的方向为最小曲率的方向。深色的长曲线为我们需要获取的牙齿咬合面分界线；

图5(a)使用本发明方法获取的路径；

图5(b)对图5(a)中的路径进行平滑后的效果；

图6本发明方法和zhuang的方法的分割结果对比图；

图7本发明方法中的度量(公式3)与zhuang的Min和Max度量获得的线的比较；

图8本发明方法中的度量(公式4)与zhuang的Min度量获得的线的比较；

图9(a)未平滑过的原始结果；

图9(b)s＝2平滑后的结果；

图9(c)s＝5平滑后的结果；

图9(d)s＝10平滑后的结果；

图9(e)s＝20平滑后的结果；

图9(f)s＝80平滑后的结果；

图10-1本发明方法在第一种非牙齿模型上获得的分割线的效果；

图10-2本发明方法在第二种非牙齿模型上获得的分割线的效果。

具体实施方式

下面结合附图对本发明作进一步详细描述：

牙齿的分割对数字口腔系统十分重要。本发明提出了一个基于live-wire算法的用于牙齿分割的三维智能剪刀工具。原始live-wire算法的缺陷是：无法实时根据交互计算分割边界；需要的交互较多。为了克服原有算法的缺点，本发明提出了两个基于张量的各向异性度量，这个张量可以使获得的分割边界更加精确，并且需要更少的交互。另外，本发明提出了一个基于各向异性度量的Dijkstra算法，该算法可以快速地估计出测地距离。为了解决基于Dijkstra算法获取的边界不平滑问题，本发明提出了三维中点平滑算法。实验表明，本发明提出的分割方法对于牙齿模型的分割任务是有效的，并且本发明方法在时间复杂度和交互便利性上更优。

牙齿模型是三角形网格(即mesh或者triangle mesh)，可以通过3D扫描仪进行3D扫描获得。关于基于三角形网格的分割算法有两种思路：一种是对每个三角形进行分类，同一类别的三角形构成分割的一部分；另一种是找到分割边界并将三角形网格分割成不同的部分。本发明的智能剪刀工具采用了第二种思路。

本发明提出了一种改进的三维三角网格分割的Live-Wire算法，特别针对牙齿分割进行了优化，并设计了一种3D智能剪刀工具。用于三角形网格分割的Live-Wire算法是在三角形网格上定义两个点之间的线作为分割边界，并且关键问题是能够快速计算出这条线。另一个问题是这条线应尽可能光滑。此外，这条线应该尽量靠近突出的脊和凹陷的槽。

基于良好各向异性张量的测地线可用作为这条线。测地线是曲面上一种特殊的最短路径(可参考“Sharir M,Schorr A.On Shortest Paths in Polyhedral Spaces[J].Siam Journal on Computing.1984,15(1):144-153”)。计算精确的测地线(可参考“辛士庆，刘永进，等.网格模型上的离散测地线[J].中国科学:信息科学.2015(03):313-335”)在分割应用中是无意义的并且时间复杂度高。为了实时交互，本发明采用了一种基于图搜索(可参考“Dijkstra E W.A note on two problems in connexion with graphs[M].Springer-Verlag New York,Inc.,1959”)的估计方法计算测地距离。由于基于图搜索的路径要沿着图的边缘，因此它是不平滑的，而本发明需要的是光滑的线。

为了应对这一挑战，本发明提出了一种基于图搜索的各向异性测地线计算方法，以及一种平滑三维曲线方法，以解决计算测地距离复杂度高和基于图搜索路径的不平滑问题。这项工作的主要贡献是三方面：

·开发了一种用于牙齿分割的智能剪刀工具，该工具使用了改进的Live-wire算法，时间复杂度较低。该工具满足实时交互的要求，并且特别针对牙齿分割进行了优化。

·该工具需要较少的交互，可以获得更好的分割效果。

·提出了一种三角网格上三维离散曲线平滑算法，巧妙地将曲面曲线平滑化问题转化到平面上。

基于live-wire算法进行网格模型分割的主要思想是在曲面上定义两个点之间的一条线，这条线作为分割边界，这个方法的关键在于快速地计算出这条线。基于“合适的”各向异性度量的测地线可以作为这条线。测地线是曲面上两点之间的最短曲线。首要问题是计算测地线。在三角网格模型上计算测地线有两类方法：精确的方法和估计的方法。

O’Rourke et al.(可参考“O'Rourke J,Suri S,Booth H.Shortest paths onpolyhedral surfaces[C].Symposium on Theoretical Aspects of Computer Science,1985.1985:243-254.”)首先提出了一个普通多边形网格上计算精确测地距离的算法，其时间复杂度是O(n⁵)，由于这个原因该算法并没有被实现。Chen et al.(可参考“Chen J,HanY.Shortest paths on a polyhedron[C].symposium on computational geometry,1990.1990:360-369”)基于一个角产生一个窗口(one angle one split)的思想提出了一个方法，将时间复杂度降为O(n²)。MMP(可参考“Surazhsky V,Surazhsky T,Kirsanov D,etal.Fast exact and approximate geodesics on meshes[J].ACM Transactions onGraphics.2005,24(3):553”)算法是基于连续Dijkstra思想的计算精确测地距离的方法，该算法的平均时间复杂度是O(n^1.5logn)，但是在最坏情况下时间复杂度是O(n²logn)。Yinet al.(可参考“Surazhsky V,Surazhsky T,Kirsanov D,et al.Fast exact andapproximate geodesics on meshes[J].ACM Transactions on Graphics.2005,24(3):553”)提出了SVG(saddle vertex graph)算法，该算法预构建一个鞍点图，然后可以通过Dijkstra算法计算最短路径，他的时间复杂度是O(Dnlogn),D＝n。

总之，精确的测地线计算方法时间复杂度高，尽管SVG算法在时间复杂度上要低很多，然而这些方法都需要满足三角不等式，并且在较大各向异性的情况下，这个条件通常不能被满足。

估计的方法计算速度较快，Geodesic based Heat Method(可参考“Crane K,Weischedel C,Wardetzky M.Geodesics in heat[J].ACM Transactions onGraphics.2013,32(5):1-11”和“Crane K,Weischedel C,Wardetzky M.The heat methodfor distance computation[J].Communications of the ACM.2017,60(11):90-99”)是当前最先进的估计的方法。Crane等人提出了一种基于Heat Equation的方法求解测地线，该算法是一种近似算法，在使用各项同性度量时计算速度快，并且得到的测地线是光滑的，但测地线不能沿着曲面特征线走。Yang等人(可参考“Yang F,Cohen L D.Geodesic Distanceand Curves Through Isotropic and Anisotropic Heat Equations on Images andSurfaces[J].Journal of Mathematical Imaging and Vision.2016,55(2):210-228”)基于各向同性度量，使用可变的热传导系数，获取了可以沿着“脊”和“凹槽”的测地线，但是计算不稳定，在复杂模型上，回溯测地线时，容易导致测地线中断。Andreux等人(可参考“Andreux M,RodolàE,Aubry M,et al.Anisotropic Laplace-Beltrami Operators forShape Analysis[C].The Workshop on Non-Rigid Shape Analysis&Deformable ImageAlignment,2014.2014:299-312”)提出了各向异性的拉普拉斯算子，结合Heat Method，可以通过求解各向异性热传导方程，得到特征敏感的测地线，然而较大的各向异性(具有明显脊和凹槽的网格模型通常在这些地方都具有较大的各向异性)，将导致计算错误，算法鲁棒性不足。

另外一种估计的方法是Short-Term Vector Dijkstra(STVD)(可参考“Campen M,Heistermann M,Kobbelt L.Practical Anisotropic Geodesy[J].Computer GraphicsForum.2013,32(5):63-71”)，把三角形网格看作一个无向图，使用Dijkstra算法，计算源点到其它所有顶点的最短距离，每条边使用各向异性度量，这样可以估计出最短测地距离(往往这个距离是偏大的)，然后可以回溯出任意一点到源点的测地路径。由于使用各项异性度量，所以可以保证获取到的测地线能够沿着脊和凹槽走。但是这种方法得到的测地路径是沿着三角形网格的边的，不光滑，但是计算效率较高，算法复杂度为O(nlogn)。

从上面的分析可以看出，现有的方法要不就是计算代价高，要不就是需要满足三角不等式，要不就是结果不平滑。为了克服这些缺点，本发明提出了一个各向异性Dijkstra算法和一个平滑三角网格上3D曲线的算法。

类似于图像分割中的“智能剪刀工具”，本发明的分割工具以用户交互地在三角形网格上选择一个种子点开始，然后实时地计算出种子点和当前鼠标位置之间的路径。这个路径应该尽量地贴合三角形网格的特征位置。本发明使用曲面上两点之间的测地线，作为用户交互输入的路径。为了让该路径贴合在曲面的特征位置(脊和凹槽)，本发明使用基于各向异性度量的方法计算三角形网格中每条边的长度。

由于这是一个实时的交互式分割工具，因此需要较高的计算效率，并且计算得到的路径需要尽量平滑。

基于以上约束本发明提出了一种基于各向异性度量的平滑测地线计算方法，其原理如下：

·计算三角形网格的最大、最小曲率(包括曲率的大小和方向)和法向量(可参考“Panozzo,D.,Puppo,E.,&Rocca,L.(2010).Efficient multi-scale curvature andcrease estimation.In 2nd International Workshop on Computer Graphics,ComputerVision and Mathematics,GraVisMa 2010-Workshop Proceedings(pp.9-16)”和“Jin,S.,Lewis,R.R.,&West,D.(2005).A comparison of algorithms for vertex normalcomputation.The Visual Computer,21(1),71-82.”)，采用下面的公式(3)和(4)计算各向异性度量张量；

·将三角形网格看做一个无向图，使用Dijkstra算法，计算种子点与其他所有点的距离(即下面的各向异性Dijkstra算法)，边的权重使用基于各向异性度量的方法计算(采用下面的公式(6))。

·回溯两点之间的最短路径path。

·计算path的k邻域，得到一个新的三角网格，记为k_ring_mesh。

·使用保形映射(可参考“Lévy B,Petitjean S,Ray N,et al.Least squaresconformal maps for automatic texture atlas generation[J].ACM Transactions onGraphics.2002,21(3):362-371”)的方法，将k_ring_mesh映射为二维的三角网格，记为k_ring_mesh_2D。

·在二维平面上平滑得到的路径，然后再将该路径映射到原始的三维网格上。

对于细节较为简单的曲面，可以输入较少的点就能够得到较好的分割结果；对于复杂的细节丰富的曲面，用户可以通过输入更多的点来获取完全满足用户想法的分割结果。

具体如下：

各向异性度量：

本发明使用曲面上两点之间的“最短”路径(即：测地线)作为分割边界的一部分。显然，这里的“最短”不是通常意义(欧氏度量)上的“最短”。为了定义“最短”的概念，引入度量张量g。在黎曼流形M上，点a的切平面T_aM是对该流形的近似表示，如果点b是点a邻域上的另外一点，则在该流形上a和b之间的基于度量张量g_a的距离可以定义为：

其中，x是点a和点b这两点构成的向量，表示从b到a的向量，为了理解度量张量g_a，可以把它看做是在切平面T_aM上定义了一个局部正交坐标系(如图2所示)。e₁和e₂是g_a的特征向量，λ₁和λ₂是它的特征值。θ是x和e₁之间的夹角，那么公式(1)可以被重写为：

从上式可以看出，如果特征向量固定，则两点之间距离的大小将由λ₁和λ₂.决定。因此可以通过控制λ₁和λ₂.的大小，使得两点之间的距离在特征区域较小。

下面将以上想法拓展到三角网格曲面上。特征向量应该与曲面局部的曲率(最大和最小曲率)方向对齐。在本发明中，将e₁的方向与最小曲率的方向对齐，即e₁为最小曲率的方向，将e₂的方向与最大曲率的方向对齐，即e₂为最大曲率的方向。k_min和k_max代表最小曲率和最大曲率的大小。|k_max|-|k_min|在脊处较小，在凹槽处较大，并且|k_max|+|k_min|在凹槽和平坦区域较小(如图3(a)和图3(b)所示)。凹槽处曲线的切线方向与最小曲率方向相近(如图4(a)所示)，因此，为了让凹槽处的曲线更短，将λ₁和λ₂按照公式(3)设置。

除以k_max可以使λ₁和λ₂是尺度无关的。从图3(b)可以看出在凹槽区域λ₁更小，从图4(a)可以看出在凹槽区域sin(θ)更小。通过公式(2)可知，使用公式(3)的特征值可以使得沿着凹槽处的曲线更短。为了让曲线能够沿着脊走，可以将λ₁和λ₂按照公式(4)设置。图4(a)和图4(b)中的长曲线是本发明需要获取的分割边界线。图4(a)和图4(b)中的深色短线是最小曲率的方向，浅色短线是最大曲率的方向。

从图3(a)可以看出，在脊的区域|k_max|和|k_min|的大小是相近的，因此在这个区域λ₂更小，并且在其它区域更大。另外，沿着脊的线的切线方向与最大曲率方向相近(如图4(b)所示)，较小的λ₂使得沿着脊的线更短。在公式(4)中ε(本发明中设置为10^-4)是一个较小的实数，用于防止被零除。从以上定义中可以看出，λ₁和λ₂是尺度无关的。

为了高效地计算所有顶点到源点的最短路径，本发明提出了一种各向异性Dijkstra算法，具体如下：

算法的一个输入为三角网格模型，其采用half-edge(可参考“Botsch M,Steinberg S,Bischoff S,et al.OpenMesh:A generic and efficient polygon meshdata structure[J].Opensg Symposium.2002”)结构存储。三角网格的每个顶点包含一个三元组(dist,final,pred)，其中dist表示该顶点到源点的距离，final是一个Boolean类型的标记，用于防止重复计算，pred表示在一条路径中该顶点的前一个顶点。

·算法首先将源点的距离(dist字段)设置为0，将其他顶点离源点的距离(dist字段)设置为无穷大；

·算法创建一个优先队列Q，Q按照顶点的dist字段进行排序；

·将源点加到Q中；

·当Q不为空时进行如下操作；

■从Q中取出一个顶点v，将该顶点v的final字段设置为true；

■对于所有顶点v相邻的顶点w，如果w的final字段为false，进行如下操作；

◆如果w.dist>v.dist+l_g(v,w)，则w.dist＝v.dist+l_g(v,w),w.pred＝v。并且，如果w不在Q中，则将w加入到Q中。

算法根据Q是否为空进行迭代操作，不断地更新每个顶点的dist和pred字段，直到迭代结束。在迭代过程中需要计算相邻顶点之间的距离，如果两个顶点分别是v,w，其对应的度量张量分别是g_v,g_w，则v,w两点间的距离可以使用如下公式计算：

在路径回溯时，通过pred字段遍历顶点，直到到达源点，这样就可以方便地找到任意顶点到源点的最短路径，即基于各向异性度量的测地线，效果如图5(a)所示。

具体使用本发明方法时，是将种子点作为源点，dist即该顶点与所述种子点之间的基于各向异性度量的最短距离。

上述各向异性Dijkstra算法的伪代码如下：

平滑牙齿网格模型上的三维曲线：

从图5(a)可以看到，得到的测地线是锯齿状的，这是由Dijkstra算法决定的，回溯的路径只能沿着三角形网格的边走，而三角形网格分割(mesh segmentation)希望得到平滑的分割结果，也就是要获取平滑的测地线。为了获取平滑的测地线，Polthier等人(可参考“Polthier K,Schmies M.Straightest geodesics on polyhedral surfaces[C].2006.ACM,2006:30-38”)基于测地距离场的梯度信息获取平滑的测地线，从而提出了Geodesic Euler Method和Geodesic Runge-Kutta Method。这两个方法将测地线回溯问题看做是一个带初值常微分方程的数值解问题，这类算法对测地距离场的要求较高，需要测地距离场较为平滑，而由于较强的各向异性以及Dijkstra算法精确性不足，使得使用各向异性Dijkstra算法计算出的测地距离场无法满足平滑条件，因此导致这类算法无法得到正确的测地线。

(1)从图5(a)可以看到，测地线虽然是锯齿状的，但是其走势是正确的。

(2)在平面上，不断取折线的每条折线段的中点可以使该折线平滑，然而在三角网格上，曲线无法直接使用类似于平面曲线平滑方法的原因在于空间曲线两点的中点不一定在三角网格的面上。

正是基于以上两点，本发明提出了一个用于平滑三角网格模型上曲线的三维中点平滑算法，具体步骤如下：

首先将曲线以及曲线经过的边的k邻域分割出来(采用现有的方法即可实现)，记作k_ring_mesh；

然后将k_ring_mesh通过保形映射(可参考“Lévy B,Petitjean S,Ray N,etal.Least squares conformal maps for automatic texture atlas generation[J].ACMTransactions on Graphics.2002,21(3):362-371”)的方法展开到平面上(记作k_ring_mesh_2D)，这样这条曲线(记作C)也就随之变成了2D曲线，记作C^2D；

在平面上使用中点法平滑这条曲线，平滑强度记为s，平滑后的2D曲线记作S_C^2D；

最后再将S_C^2D映射回原始的三角网格上，这条曲线记作S_C。计算方法如下：

计算S_C^2D与k_ring_mesh_2D每条边的交点，如果S_C^2D的一条边与k_ring_mesh_2D相交的边的两个顶点为index1,index2，交点为pt_int^2D，则将交点信息记在一个三元组(index1,index2,α)中，其中α(指交点与两个顶点

和

的位置关系)可以通过公式(6)计算出来，然后根据公式7计算出pt_int^3D，所有的pt_int^3D构成曲线S_C。

pt_int^2D＝(1-α)·pt^2D(index1)+α·pt^2D(index2) (6)

pt_int^3D＝(1-α)·pt^3D(index1)+α·pt^3D(index2) (7)

公式(6)和(7)中pt^2D,pt^3D是k_ring_mes_2D和k_ring_mesh的顶点，所有的pt_int^3D构成了S_C。平滑后的效果如图5(b)所示。

其伪代码如下：

用户交互接口：

本发明方法使用C++语言开发，依赖于ligigl(可参考“Jacobson A,Panozzo D,Others.libigl-A simple C++geometry processing library[Z].2018.”)，其主要提供了数据读写、渲染、UI组件和普通的离散几何量的计算(Discrete Geometric Quantities)等功能。用户通过Ctrl键和鼠标左键点击在三角网格上选择一个种子点，开始一条分割边界的输入。当Ctrl键按下并且鼠标在模型上移动时，该工具寻找鼠标当前点和种子点之间的最短路径，并显示出来。当用户认为这条路径是最优时，通过鼠标左键点击来确认这条路径。当路径封闭时，结束一条分割边界的输入，程序根据输入的分割边界进行分割。程序保存着用户输入的所有种子点，当鼠标悬停在种子点上时，种子点将高亮显示，用户可以通过拖动这些种子点进行分割边界的修改，拖动的过程中分割边界实时更新显示。在没有长按Ctrl键时，用户可以使用鼠标进行普通的3D交互，包括：平移、旋转和缩放。牙齿与牙龈、牙齿的颌面部的分割都是牙齿分割的任务。在这些分割任务中，有的分割边界是沿着脊的线，有的需要沿着凹槽的线。因此，在程序界面上有一个单选按钮，用于选择不同的模式。通过操作界面，可以改变程序内部使用的各向异性度量。

如图1所示，实际使用本发明方法时的具体步骤如下：

(1)读取牙齿模型：通过3D扫描获得牙齿模型，所述牙齿模型是三角形网格；

(2)计算三角形网格中各个顶点的各向异性度量张量，具体如下：

(21)计算三角形网格中的每个顶点处的最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂和法向量；

(22)根据最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂获得三角形网格中的每个顶点的各向异性度量张量。

所述步骤(22)的操作包括：

三角形网格中的每个顶点对应一个各向异性度量张量，每个所述各向异性度量张量的特征向量是e₁和e₂，每个所述各向异性度量张量的特征值是λ₁和λ₂；

如果进行齿龈分割，则使用公式(3)计算得到特征值λ₁和λ₂：

如果进行牙齿咬合面分割，则使用公式(4)计算得到特征值λ₁和λ₂：

(3)利用鼠标选择种子点：利用鼠标在分割边缘处(不用精确到边界上的某个位置，在边界附近就可以)选择一个点作为种子点，具体如下：如果进行齿龈分割，则在齿龈交界处选择一个点作为种子点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点作为种子点。

(4)计算三角形网格中各个顶点与所述种子点之间的基于各向异性度量的最短距离，具体如下：

采用half-edge结构存储所述三角形网格，所述三角形网格中的每个顶点包含一个三元组(dist,final,pred)，其中dist表示该顶点与所述种子点之间的基于各向异性度量的最短距离，final是一个Boolean类型的标记，用于防止重复计算，pred表示在一条路径中该顶点的前一个顶点；

求取各个顶点的三元组(dist,final,pred)的步骤如下：

(41)将种子点的dist字段设置为0，将其他顶点的dist字段设置为无穷大；

(42)创建一个优先队列Q，Q按照各个顶点的dist字段进行排序；

(43)将种子点加到Q中；

(44)当Q不为空时进行如下操作：

(441)从Q中取出一个顶点v，将该顶点v的final字段设置为true；

(442)对于所有与顶点v相邻的顶点w，如果w的final字段为false，则进行如下操作：

(4421)如果w.dist>v.dist+l_g(v,w)，则w.dist＝v.dist+l_g(v,w),w.pred＝v；并且，如果w不在Q中，则将w加入到Q中；

其中，l_g(v,w)表示两个顶点v、w之间的距离，两个顶点v、w对应的各向异性度量张量分别是g_v,g_w，采用公式(5)计算v、w两点间的距离l_g(v,w)：

其中，

均采用公式(2)计算得到：

其中，x表示从顶点b到顶点a的向量，θ是x和e₁之间的夹角。

(5)根据各个顶点与所述种子点之间的基于各向异性度量的最短距离实时获得鼠标当前位置与种子点之间的路径：找到与鼠标当前位置最近的顶点，从该顶点开始，通过各个顶点的pred字段遍历各个顶点，直到找到的顶点是种子点为止，这样就获得了鼠标当前位置与种子点之间的路径。

(6)对步骤(5)得到的路径进行平滑得到平滑的路径，具体如下：

(61)将步骤(5)得到的路径记作曲线C，将曲线C以及曲线C经过的边的k邻域分割出来，记作k_ring_mesh；

(62)将k_ring_mesh通过保形映射的方法展开到平面上,记作k_ring_mesh_2D,这样，曲线C变成了2D曲线C^2D；

(63)在平面上使用中点法对2D曲线C^2D进行平滑处理，平滑强度记为s，平滑处理后的2D曲线为S_C^2D；

(64)将S_C^2D映射回原始的三角形网格上，成为曲线S_C，即平滑后的路径，具体如下：

计算S_C^2D与k_ring_mesh_2D每条边的交点，如果S_C^2D的一条边与k_ring_mesh_2D相交的边的两个顶点为index1,index2，交点为pt_int^2D，则将交点信息记在一个三元组(index1,index2,α)中，其中α表示交点与两个顶点

和

的位置关系，利用公式(6)计算得到α：

pt_int^2D＝(1-α)·pt^2D(index1)+α·pt^2D(index2) (6)

然后根据公式(7)计算出pt_int^3D：

pt_int^3D＝(1-α)·pt^3D(index1)+α·pt^3D(index2) (7)

其中，pt^2D,pt^3D是k_ring_mes_2D和k_ring_mesh的顶点；

所有的pt_int^3D构成了曲线S_C。

(7)利用鼠标点击确认路径：利用鼠标在分割边缘处选择一个点，并把该点作为种子点；利用鼠标在分割边缘处选择一个点的方法如下：如果进行齿龈分割，则在齿龈交界处选择一个点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点

(8)判断路径是否封闭，如果是，则转入步骤(9)，如果否，则返回步骤(4)；所述判断路径是否封闭的操作包括：

判断步骤(7)选择的种子点与第一个种子点是否重合，如果重合，则判定路径封闭，如果不重合，则判定路径不封闭。

(9)利用所述平滑的路径分割牙齿。

利用本发明方法进行的实验如下：

本发明的实验数据全部来自于临床病例。通过传统的石膏印模技术获取病人的牙齿模型，然后通过三维扫描仪(AutoScan-DS100+)获取数字牙齿模型。程序使用C++编写，通过Microsoft Visual Studio 2017编译生成可执行程序。运行本实验中测试程序的PC配置为：8GB RAM、Intel^R Core^TM i7-4790 3.6GHz CPU、Windows 10 64位操作系统。实时交互性和易用性对这个分割工具非常重要，因此实验也主要从这两个方面对本发明的方法和Zhuang的方法进行评估。

名称	顶点的数量	三角面的数量	边的数量	位置
					DM1	35473	70942	106413	下颌
DM2	113489	226843	340325	上颌
					DM3	84091	168194	252291	上颌
DM4	140251	280498	420747	下颌
					DM5	156257	312510	468765	上颌
DM6	226232	452460	678690	下颌
					DM7	62301	124598	186897	下颌

表1

为了评估本发明的方法，在7个牙齿模型(分别为DM1到DM7)上测试了本发明的方法，这些模型的基本信息如表1所示。本发明的方法获得的结果和zhuang的方法获取的结果如图6所示(图6中的视图1、视图2分别表示两个视图，是为了多视角展示实验结果。由于三维空间中存在遮挡，一个视角无法看到全部的结果，所以使用了两个视角。)。在该测试中，本发明方法的参数分别设置为k＝3和s＝5。Zhuang的方法使用他提出的Max度量，其他参数按照他论文中的建议进行设置。这两个方法都包含三个操作：初始化、交互式选择种子点和回溯路径。每个操作耗费的时间如表2所示。在该实验中，选择种子点的原则是：两个相邻的种子点要尽量地远，并且要保证得到的相邻两个种子点之间的线尽量贴合在牙齿和牙龈的交界处。如图6所示，本发明方法获取的分割线和zhuang的方法获取的分割线基本相同，然而本发明方法需要的种子点要少于zhuang的方法需要的种子点(如表3所示)(由于本方法提出的度量(参见公式(3))，所以需要的种子点更少。)。从表2中可以看出，本发明的方法的seeding步骤(指选择种子点后，“计算所有顶点与种子点的距离”步骤；对应于“各向异性Dijkstra算法”)需要的时间要远小于zhuang的方法(因为“各向异性Dijkstra算法”是估计类方法，所以需要的时间更少，效率更高)，虽然本发明方法在路径回溯步骤中需要的时间略高于zhuang的方法，但是需要的时间也仅仅是几毫秒，不影响实时交互。因此本发明的方法是“所见即所得”的，具有较好的用户体验。另外，本发明方法的初始化时间也比zhuang的方法要少，因为本发明的方法不需要对三角网格模型进行局部细分，并且由于度量时是尺度无关的，也不需要将模型归一化。

表2

	DM1	DM2	DM3	DM4	DM5	DM6	DM7
								本发明的方法	31	33	35	28	30	26	36
Zhuang的方法	51	47	47	48	57	39	48

表3

从图7可以看出，使用本发明提出的第一个度量(即公式(3))获得的分割线离牙齿和牙龈的分界线更近(3#线是使用本发明提出的第一个度量得到的，1#和2#线分别是使用zhuang的Min度量和Max度量得到的)。在图7左上角中zhuang的一条线(2#)和本发明的线(1#)是重合的，表示此时这两种方法效果是相同的。但是从图8其它部分可以看出，当种子点离得较远时，本发明的方法已经明显优于zhuang的方法。并且当两个种子点较远时，zhuang的方法已经无法获得正确的结果(图8下方中zhuang的方法的两条线1#和2#均与牙龈处相差很远)，因此本发明提出的第一个度量更加有效。图8为牙齿咬合面的分割结果，其中1#线是zhuang方法的结果(Zhuang提出的Min度量)，2#线是我们方法的结果(本发明提出的第二个度量，即公式(4))，可见这两个度量的效果是相似的，但是从上边分析可知，本发明的方法在计算效率上更高。

本发明的三维曲线平滑算法包含两个参数。参数k控制着平滑曲线时使用的k_ring_mesh的大小(曲线C的k邻域构成k_ring_mesh，因此k越小k_ring_mesh的规模就越小)，这个参数一般设置为3到5之间。较小的k可能导致平滑后的曲线的点不在k_ring_mesh上，最终导致部分曲线不在曲面上。但是较大的k将导致计算效率降低，因为需要计算曲线与k_ring_mesh每条边的交点。通过实验我们发现k＝3是一个较为合适的参数。参数s代表平滑强度，该参数影响取中点的次数(参见上述的平滑算法“Algorithm:3D Midpointsmoothing algorithm”)。在平面上，不断取一条折线的中点，这条折线将最终变为直线，在曲面向效果也是相似的。因此一个合适的参数s对于获取平滑的分割线时十分重要的。从图9(a)到图9(f)中可以看出，过大的s将导致获得的线远离实际的边界。通过实验，我们发现s最好设置在2到10之间。从图9(a)到图9(f)中可以看出，只要s在一个合理的范围内，本发明的方法对s是不敏感的。

虽然本发明的方法是针对牙齿分割提出的，但是该方法也能用于有明显脊和凹槽的模型的分割，效果如图10-1和图10-2所示，其中的点为种子点，线条为利用本发明方法得到的分割线。

本发明改进了zhuang的live-wire算法，提出了一个用于牙齿模型的分割的智能剪刀工具。为了获取分割边界，本发明定义了两个度量。为了快速地计算测地距离，本发明提出了一个基于各向异性度量的Dijkstra算法。为了解决基于Dijkstra算法获得的路径不平滑的问题，本发明提出了一个三维中点平滑算法，用于平滑三角网格曲面上的曲线。实验表明，本发明方法对于牙齿分割任务是有效的。与zhuang的方法相比，本发明的方法在时间复杂度和可相互性上更优。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (10)

1.一种基于各向异性度量的牙齿分割方法，其特征在于：所述方法包括：

(1)读取牙齿模型：通过3D扫描获得牙齿模型，所述牙齿模型是三角形网格；

(2)计算三角形网格中各个顶点的各向异性度量张量；

(3)利用鼠标选择种子点；

(4)计算三角形网格中各个顶点与所述种子点之间的基于各向异性度量的最短距离；

(5)根据各个顶点与所述种子点之间的基于各向异性度量的最短距离实时获得鼠标当前位置与种子点之间的路径；

(6)对步骤(5)得到的路径进行平滑得到平滑的路径；

(7)利用鼠标点击确认路径；

(8)判断路径是否封闭，如果是，则转入步骤(9)，如果否，则返回步骤(4)；

(9)利用所述平滑的路径分割牙齿。

2.根据权利要求1所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(2)的操作包括：

(21)计算三角形网格中的每个顶点处的最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂和法向量；

(22)根据最大曲率k_max、最小曲率k_min、最小曲率的方向e₁、最大曲率的方向e₂获得三角形网格中的每个顶点的各向异性度量张量。

3.根据权利要求2所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(22)的操作包括：

三角形网格中的每个顶点对应一个各向异性度量张量，每个所述各向异性度量张量的特征向量是e₁和e₂，每个所述各向异性度量张量的特征值是λ₁和λ₂；

如果进行齿龈分割，则使用公式(3)计算得到特征值λ₁和λ₂：

如果进行牙齿咬合面分割，则使用公式(4)计算得到特征值λ₁和λ₂：

4.根据权利要求3所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(3)的操作包括：

利用鼠标在分割边缘处选择一个点作为种子点，具体如下：

如果进行齿龈分割，则在齿龈交界处选择一个点作为种子点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点作为种子点。

5.根据权利要求4所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(4)的操作包括：

采用half-edge结构存储所述三角形网格，所述三角形网格中的每个顶点包含一个三元组(dist,final,pred)，其中dist表示该顶点与所述种子点之间的基于各向异性度量的最短距离，final是一个Boolean类型的标记，用于防止重复计算，pred表示在一条路径中该顶点的前一个顶点；

求取各个顶点的三元组(dist,final,pred)的步骤如下：

(41)将种子点的dist字段设置为0，将其他顶点的dist字段设置为无穷大；

(42)创建一个优先队列Q，Q按照各个顶点的dist字段进行排序；

(43)将种子点加到Q中；

(44)当Q不为空时进行如下操作：

(441)从Q中取出一个顶点v，将该顶点v的final字段设置为true；

(442)对于所有与顶点v相邻的顶点w，如果w的final字段为false，则进行如下操作：如果w.dist>v.dist+l_g(v,w)，则w.dist＝v.dist+l_g(v,w),w.pred＝v；并且，如果w不在Q中，则将w加入到Q中；

其中，l_g(v,w)表示两个顶点v、w之间的距离，两个顶点v、w对应的各向异性度量张量分别是g_v,g_w，采用公式(5)计算v、w两点间的距离l_g(v,w)：

其中，

均采用公式(2)计算得到：

其中，x表示从顶点b到顶点a的向量，θ是x和e₁之间的夹角。

6.根据权利要求5所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(5)的操作包括：

找到与鼠标当前位置最近的顶点，从该顶点开始，通过各个顶点的pred字段遍历各个顶点，直到找到的顶点是种子点为止，这样就获得了鼠标当前位置与种子点之间的路径。

7.根据权利要求6所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(6)的操作包括：

(61)将步骤(5)得到的路径记作曲线C，将曲线C以及曲线C经过的边的k邻域分割出来，记作k_ring_mesh；

(62)将k_ring_mesh通过保形映射的方法展开到平面上,记作k_ring_mesh_2D,这样，曲线C变成了2D曲线C^2D；

(63)在平面上使用中点法对2D曲线C^2D进行平滑处理，平滑强度记为s，平滑处理后的2D曲线为S_C^2D；

(64)将S_C^2D映射回原始的三角形网格上，成为曲线S_C，即平滑后的路径。

8.根据权利要求7所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(64)的操作包括：

计算S_C^2D与k_ring_mesh_2D每条边的交点，如果S_C^2D的一条边与k_ring_mesh_2D相交的边的两个顶点为index1,index2，交点为pt_int^2D，则将交点信息记在一个三元组(index1,index2,α)中，其中α表示交点与两个顶点pt^2D _index1和pt^2D _index2的位置关系，利用公式(6)计算得到α：

pt_int^2D＝(1-α)·pt^2D(index1)+α·pt^2D(index2) (6)

然后根据公式(7)计算出pt_int^3D：

pt_int^3D＝(1-α)·pt^3D(index1)+α·pt^3D(index2) (7)

其中，pt^2D,pt^3D是k_ring_mes_2D和k_ring_mesh的顶点；

所有的pt_int^3D构成了曲线S_C。

9.根据权利要求8所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(7)的操作包括：

利用鼠标在分割边缘处选择一个点，并把该点作为种子点；

利用鼠标在分割边缘处选择一个点的操作包括：

如果进行齿龈分割，则在齿龈交界处选择一个点；如果进行牙齿咬合面分割，则在牙齿咬合面的分割处选择一个点。

10.根据权利要求9所述的基于各向异性度量的牙齿分割方法，其特征在于：所述步骤(8)中的判断路径是否封闭的操作包括：

判断步骤(7)得到的种子点与第一个种子点是否重合，如果重合，则判定路径封闭，如果不重合，则判定路径不封闭。

Images