CN105069777A - 一种预备体网格模型的颈缘线自动提取方法 - Google Patents

Abstract

本发明公开了一种预备体网格模型的颈缘线自动提取方法，包括：S1、计算预备体网格模型顶点的平均曲率；S2、采用基于K-means聚类的方法处理所得曲率值，得到特征点集F；S3、采用一种改进的折线生长技术，得到折线点集合F_z，进而生成特征线；S4、优化特征线；S5、以优化后的特征线上的点作为控制点，拟合成颈缘线。本发明具有提取精度高、鲁棒性好且自动化程度高的优点，可广泛应用于口腔修复领域。

Description

一种预备体网格模型的颈缘线自动提取方法

技术领域

本发明属于口腔修复领域，涉及计算机应用技术和生物医学工程。

背景技术

随着光学测量、CAD/CAM技术的发展，从患者牙颌模型数字化测量到快速计算机辅助设计制作的应用，口腔修复体制作的效率和质量不断提高，口腔CAD/CAM修复技术及相应系统的研制成为研究的热点。在口腔CAD/CAM修复过程中，预备体颈缘线的提取质量直接影响修复体的建模精度及产品质量，而且义齿在口腔内就位时，其边缘的位置、形态等严重影响患者佩戴修复体的舒适度和牙周组织健康。因而，准确、灵活地提取预备体的颈缘线对于修复体的成功制作至关重要，预备体颈缘线的提取技术研究意义重大。

目前，国外已有一批先进的口腔修复CAD/CAM系统，但颈缘线提取算法大多未知。国内也出现一些颈缘线提取算法，例如，戴宁等在《牙齿预备体颈缘线自适应提取算法的研究与实现》论文中借助数理统计学计算概率引导线、搜索球半径和质心修正，动态搜索关键特征点，最终拟合成颈缘特征线。该方法实现了颈缘线的自动提取，但提取精度有待提高。而且构造局部的引导线会被附近的非颈缘线但是符合特征的点所引导，会导致个别牙齿无法正确地完成颈缘线提取，不够可靠；另外，聂欣在《口腔修复金属基底冠CAD关键技术研究与实现》论文中利用“方向追踪”法产生初始轮廓线，再用Snake能量模型将轮廓线进化成特征线，张长东等在《基于启发式搜索策略的牙齿生物特征线提取技术》中提出了一种基于启发式搜索策略的牙齿生物特征线提取技术，也都实现了颈缘线的提取，但在提取过程中都需要人工判断颈缘线的位置并交互拾取一个甚至多个特征点，这种提取方式对操作人员的有一定要求，提取效率也相对较低。

发明内容

为克服上述现有技术问题，本发明提出一种精度高、鲁棒性好且自动化程度高的，预备体网格模型的颈缘线提取方法。

本发明解决其技术问题所采用的技术方案是：

一种预备体网格模型的颈缘线自动提取方法，其特征在于，包括：

S1、计算预备体网格模型顶点的平均曲率；

S2、采用基于K-means聚类的方法处理所得曲率值，得到特征点集F；

S3、采用一种改进的折线生长技术，得到折线点集合F_z，进而生成特征线；

S4、优化特征线；

S5、以优化后的特征线上的点作为控制点，拟合成颈缘线。

进一步，所述步骤S1采用一种改进的Voronoi法计算预备体网格模型顶点的平均曲率，包括：

S11、根据Laplace-Beltrami算子性质和高斯定理求顶点p_i的平均曲率向量，所述的平均曲率向量计算公式为：

$K\left(p_i\right)=\frac{1}{2A_M}\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}(\cot {\mathrm{\α}}_{\mathrm{ij}}+\cot {\mathrm{\β}}_{\mathrm{ij}})(p_i-p_j)$

其中，A_M为顶点p_i的1-邻域三角面片的Voronoi区域面积之和；N(i)为与顶点p_i直接相连的顶点索引的集合；p_j为与顶点p_i直接相连的顶点；α_ij＝＜p_ip_j-1p_j，β_ij＝＜p_ip_j+1p_j，α_ij、β_ij为边p_ip_j的两对角。

S12、求顶点p_i处的单位法矢，所述的单位法矢计算公式为：

$n_{\mathrm{pi}}=\frac{{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}}{\left|\right|{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}\left|\right|}$

其中，N_fk为三角面片f_k的法矢；A_fk表示三角面片f_k的Voronoi区域面积。

S13、根据Laplace-Beltrami算子K(p_i)＝2K_H(p_i)n_pi计算出顶点p_i处的平均曲率，所述的平均曲率的计算公式为：

$K_H\left(p_i\right)=\frac{1}{2}\left|\right|K\left(p_i\right)\left|\right|\mathrm{sgn}(K\left(p_i\right)\·n_{\mathrm{pi}})$

其中，sgn(·)为符号函数。

进一步，所述步骤S11中所述的顶点p_i的1-邻域任意三角面片PQR(P为顶点p_i)的Voronoi区域面积计算公式为：

其中，S_Δ为三角面片PQR的面积。

进一步，所述步骤S2，包括：

S21、初始化：选定平均曲率作为数据空间进行聚类，从中选择K个正曲率值作为初始聚类中心，并设定聚类中心变化阈值；

S22、对于数据空间中的平均曲率，根据它们与这些聚类中心的距离，按距离最近的准则分别将它们分配给与其最相似的聚类中心所代表的类；

S23、计算每个类别中所有曲率的均值作为该类别的新聚类中心，计算数据空间中所有曲率到其所在类别聚类中心的距离平方和，即J(C)值；

S24、得到新的聚类中心后，将依据新的聚类中心对待测数据进行重新聚类，并不断重复此过程，直到聚类中心的变化达到某一设定阈值为止；

S25、从聚类后的结果中选择其中曲率较高的几类作为特征点。

进一步，所述步骤S3，包括：

S31、将所得特征点排序存储至队列Q中；

S32、计算生长方向v_g；

S33、在v_g方向上进行折线生长，直到找到新的生长点；

S34、找新的生长点，将其添加到集合F_z中；

S35、改变生长方向继续生长，直到算法结束；

S36、将集合F_z中的折线点连接成特征线。

进一步，所述步骤S4，包括：

S41、连接端点接近、走向相符的特征线；

S42、删除不闭合的和闭合但不属于颈缘线的特征线分支；

S43、对特征线进行光顺处理。

进一步，所述步骤S41，包括：

对特征线的端点进行检测，若与之相邻特征线的端点为其邻域点，且两特征线夹角满足阈值，则连接两端点。

进一步，所述步骤S42，包括：

S421、遍历所有特征线，得到分支点；

S422、分别以各分支点为起点，对其每个分支进行搜索，若在某一分支搜索到分支点则搜索另一分支，若在某一分支搜索到端点则删除这一分支；

S423、重复S422直至删除所有不闭合分支；

S424、对于闭合分支的分支点，通过判断分支点处各分支的走向来删除不属于颈缘线的分支。

进一步，所述步骤S43，包括：

对局部过度弯折的特征线通过替换和翻转操作进行处理。

进一步，所述步骤S5，包括：

以优化后的特征线上点作为控制点，使用k次Bezier曲线b_j∈R^k+1拟合这些点得到颈缘线。

本发明的有益效果是：采用改进的Voronoi法计算离散网格平均曲率，相比局部曲面拟合的方法降低了算法的复杂度，同时又避免了狭长及不规则三角面片对顶点曲率的影响，提高了精度和鲁棒性；采用基于K-means聚类的方法得到特征点，再采用改进的折线生长技术生成特征线，然后经过特征线优化，可以得到一条闭合且比较光顺的特征线，最终经过曲线拟合得到光顺的颈缘线。采用本发明的方法能够稳定且自动地提取出高精度的预备体颈缘线。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1为本发明一种预备体网格模型的颈缘线自动提取方法的整体流程图；

图2为本发明步骤S1的流程图；

图3为本发明步骤S2的流程图；

图4为本发明步骤S3的流程图；

图5为本发明步骤S4的流程图；

图6为本发明步骤S42的流程图；

图7为本发明的一个实施例整体流程图；

图8为本发明Voronoi法中α_ij，β_ij及A_M的定义示意图；

图9为本发明顶点p_i的1-邻域任意三角面片PQR(P为顶点p_i)的Voronoi区域面积示意图；

图10为本发明基于K-means聚类的特征点获取算法流程图；

图11为本发明改进的折线生长技术算法流程图；

图12为本发明折线生长示意图；

图13为本发明连接端点接近、走向相符的特征线的一个实施例示意图；

图14为本发明对特征线进行光顺处理的一个实施例示意图；

图15为本发明最终提取的颈缘线效果图。

具体实施方式

参照图1，一种预备体网格模型的颈缘线自动提取方法，包括：

S1、计算预备体网格模型顶点的平均曲率；

S2、采用基于K-means聚类的方法处理所得曲率值，得到特征点集F；

S3、采用一种改进的折线生长技术，得到折线点集合F_z，进而生成特征线；

S4、优化特征线；

S5、以优化后的特征线上的点作为控制点，拟合成颈缘线。

参照图2，进一步作为优选的实施方式，所述步骤S1采用一种改进的Voronoi法计算预备体网格模型顶点的平均曲率，包括：

S11、根据Laplace-Beltrami算子性质和高斯定理计算顶点p_i的平均曲率向量，所述的平均曲率向量计算公式为：

$K\left(p_i\right)=\frac{1}{2A_M}\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}(\cot {\mathrm{\α}}_{\mathrm{ij}}+\cot {\mathrm{\β}}_{\mathrm{ij}})(p_i-p_j)$

其中，A_M为顶点p_i的1-邻域三角面片的Voronoi区域面积之和；N(i)为与顶点p_i直接相连的顶点索引的集合；p_j为与顶点p_i直接相连的顶点；α_ij＝∠p_ip_j-1p_j，β_ij＝∠p_ip_j+1p_j，α_ij、β_ij为边p_ip_j的两对角。

S12、求顶点p_i处的单位法矢，所述的单位法矢计算公式为：

$n_{\mathrm{pi}}=\frac{{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}}{\left|\right|{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}\left|\right|}$

其中，N_fk为三角面片f_k的法矢；A_fk表示三角面片f_k的Voronoi区域面积。

S13、根据Laplace-Beltrami算子K(p_i)＝2K_H(p_i)n_pi计算出顶点p_i处的平均曲率，所述的平均曲率的计算公式为：

$K_H\left(p_i\right)=\frac{1}{2}\left|\right|K\left(p_i\right)\left|\right|\mathrm{sgn}(K\left(p_i\right)\·n_{\mathrm{pi}})$

其中，sgn(·)为符号函数。

进一步作为优选的实施方式，所述步骤S11中所述的顶点p_i的1-邻域任意三角面片PQR(P为顶点p_i)的Voronoi区域面积计算公式为：

其中，S_Δ为三角面片PQR的面积。

参照图3，进一步作为优选的实施方式，所述步骤S2，包括：

S21、初始化：选定平均曲率作为数据空间进行聚类，从中选择k个正曲率值作为初始聚类中心，并设定聚类中心变化阈值；

S22、对于数据空间中的平均曲率，根据它们与这些聚类中心的距离，按距离最近的准则分别将它们分配给与其最相似的聚类中心所代表的类；

S23、计算每个类别中所有曲率的均值作为该类别的新聚类中心，计算数据空间中所有曲率到其所在类别聚类中心的距离平方和，即J(C)值；

S24、得到新的聚类中心后，将依据新的聚类中心对待测数据进行重新聚类，并不断重复此过程，直到聚类中心的变化达到某一设定阈值为止；

S25、从聚类后的结果中选择其中曲率较高的几类作为特征点。

参照图4，进一步作为优选的实施方式，所述步骤S3，包括：

S31、将所得特征点排序存储至队列Q中；

S32、计算生长方向v_g；

S33、在v_g方向上进行折线生长，直到找到新的生长点；

S34、找新的生长点，将其添加到集合F_z中；

S35、改变生长方向继续生长，直到算法结束；

S36、将集合F_z中的折线点连接成特征线。

参照图5，进一步作为优选的实施方式，所述步骤S4，包括：

S41、连接端点接近、走向相符的特征线；

S42、删除不闭合的和闭合但不属于颈缘线的特征线分支；

S43、对特征线进行光顺处理。

进一步作为优选的实施方式，所述步骤S41，包括：

对特征线的端点进行检测，若与之相邻特征线的端点为其邻域点，且两特征线夹角满足阈值，则连接两端点。

参照图6，进一步作为优选的实施方式，所述步骤S42，包括：

S421、遍历所有特征线，得到分支点；

S422、分别以各分支点为起点，对其每个分支进行搜索，若在某一分支搜索到分支点则搜索另一分支，若在某一分支搜索到端点则删除这一分支；

S423、重复S422直至删除所有不闭合分支；

S424、对于闭合分支的分支点，通过判断分支点处各分支的走向来删除不属于颈缘线的分支。

进一步作为优选的实施方式，所述步骤S43，包括：

对局部过度弯折的特征线通过替换和翻转操作进行处理。

进一步作为优选的实施方式，所述步骤S5，包括：

以优化后的特征线上点作为控制点，使用k次Bezier曲线b_j∈R^k+1拟合这些点得到颈缘线。

下面结合说明书附图和具体实施例对本发明做进一步详细说明。

参照图7，本发明一种预备体网格模型的颈缘线自动提取方法的一个实施例包括以下步骤：

步骤一：读取模型

本发明中颈缘线提取的处理对象是预备体三角网格模型，其文件格式可以为.stl、.obj、.off和.ply等。

步骤二：计算顶点曲率

本发明采用一种改进的Voronoi法计算预备体网格模型顶点的平均曲率，计算时遍历网格的所有顶点，对于任意当前顶点p_i，其具体计算过程为：

(1)根据Laplace-Beltrami算子性质和高斯定理求顶点p_i的平均曲率向量，其计算公式为：

$K\left(p_i\right)=\frac{1}{2A_M}\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}(\cot {\mathrm{\α}}_{\mathrm{ij}}+\cot {\mathrm{\β}}_{\mathrm{ij}})(p_i-p_j)$

其中，A_M为顶点p_i的1-邻域三角面片的Voronoi区域面积之和；N(i)为与顶点p_i直接相连的顶点索引的集合；p_j为与顶点p_i直接相连的顶点；α_ij＝∠p_ip_j-1p_j，β_ij＝∠p_ip_j+1p_j，α_ij、β_ij为边p_ip_j的两对角，如图8所示。

如图9所示，顶点p_i的1-邻域任意三角面片PQR(P为顶点p_i)的Voronoi区域面积根据该三角形形状的不同而有所区别，具体地，根据ΔPQR是否为钝角三角形及∠P、∠Q还是∠R为钝角，其Voronoi区域面积计算公式为：

其中，S_Δ为三角面片PQR的面积。

(2)求顶点p_i处的单位法矢，所述的单位法矢计算公式为：

$n_{\mathrm{pi}}=\frac{{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}}{\left|\right|{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}\left|\right|}$

其中，N_fk为三角面片f_k的法矢；A_fk表示三角面片f_k的Voronoi区域面积，其计算公式与A_voronoi相同。

(3)根据Laplace-Beltrami算子K(p_i)＝2K_H(p_i)n_pi计算出顶点p_i处的平均曲率，其计算公式为：

$K_H\left(p_i\right)=\frac{1}{2}\left|\right|K\left(p_i\right)\left|\right|\mathrm{sgn}(K\left(p_i\right)\·n_{\mathrm{pi}})$

其中，sgn(·)为符号函数。

步骤三：获取特征点

本发明采用基于K-means聚类的方法处理所得曲率值，得到特征点集F，具体过程为：

(1)初始化：选定平均曲率作为数据空间进行聚类，从中选择K个正曲率值作为初始聚类中心，并设定聚类中心变化阈值。

具体地，由平均曲率构成的数据空间为X＝{x₁，x₂，…，x_i，…，x_n}，其中，n为网格模型顶点数目，x_i为点p_i的平均曲率；设定K₀(K₀＞0)作为阈值，从步骤一所得平均曲率选定曲率大于K₀的顶点曲率(K₀～K_max)作为K-means聚类的待测数据，将其均分为K个梯度，则相应的聚类中心然后设定聚类中心的变化阈值。作为优选的实施参数，K₀＝0.3，K＝7，将变化阈值设为0.000001。

(2)对于数据空间中的平均曲率，根据它们与这些聚类中心的距离，按距离最近的准则分别将它们分配给与其最相似的聚类中心所代表的类。

具体地，K-Means聚类算法将数据对象组织为K个划分C＝{c_k，i＝1，2，…，K}。每个划分代表一个类c_k，每个类c_k有一个类别中心μ_i。

(3)计算每个类别中所有曲率的均值作为该类别的新聚类中心，计算数据空间中所有曲率到其所在类别聚类中心的距离平方和，即J(C)值。

具体地，选取欧氏距离作为相似性和距离判断准则，计算该类内各点到聚类中心μ_i的距离平方和

$J\left(C\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}J\left(c_k\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{d}_{\mathrm{ki}}{\left|\right|x_i-{\mathrm{\μ}}_k\left|\right|}^2$

其中，

(4)得到新的聚类中心后，将依据新的聚类中心对待测数据进行重新聚类，并不断重复此过程，直到聚类中心的变化达到某一设定阈值为止。

(5)从聚类后的结果中选择其中曲率较高的几类作为特征点。

基于K-means聚类的特征点获取的算法流程图如图10所示。

步骤四：提取特征线

本发明采用一种改进的折线生长技术，得到折线点集合F_z，进而生成特征线。改进的折线生长技术算法流程图如图11所示。提取特征线的具体过程为：

(1)将所得特征点排序存储至队列Q中。

具体地，先通过PCA分析方法计算出每一个特征点在其r-邻域N_E(p)内的特征值，然后计算该点的相关系数。对于三角网格曲面上的任意点p，在其邻域N_E(p)内通过PCA分析可以得到三个特征值相关系数n依据下式计算得出

$\mathrm{\η}=\frac{{\mathrm{\λ}}_0^E}{{\mathrm{\λ}}_0^E\≥{\mathrm{\λ}}_1^E\≥{\mathrm{\λ}}_2^E}$

用同样的方法计算出F中所有特征点的相关系数，并将所有的特征点按其相关系数从大到小的顺序放入队列Q中。由此可避免选择一个具有较小相关系数的角点作为种子点进行折线生长的可能。

(2)计算生长方向v_g。

依次从Q中取出远离角点的顶点p，作为特征折线的初始生长点，即种子点。然后，以p为中心，以r_max为半径构造生长圆，取该圆内NBHD(p)＝{p_j|p_j∈F，||pj-p||≤r_max}范围内的点进行PCA分析，即计算出协方差矩阵最大特征值对应的特征向量v₀。接着，将NBHD(p)内的每个点p_j投影到向量v₀上得到p′_j，并计算距离d_j＝||p_j-p′_j||，按下式计算生长方向v_g，如图12(a)所示。

$v_g=\underset{p_j\∈\mathrm{NBHD}\left(p\right)}{\mathrm{\Σ}}{e}^{-\left(\frac{\left|\right|p_j^{\′}-p\left|\right|d_j}{r_{\max }^2}\right)}\·\frac{p_j-p}{\left|\right|p_j-p\left|\right|}$

(3)在v_g方向上进行折线生长，直到找到新的生长点。

(4)找新的生长点，将其添加到集合F_z中。

新的生长点应选在种子点p的r半径内距离种子点最远的点p_j，可使用p_j＝p+r_max·v_g计算得出，此式中r_max≤r将该点添加到特征折线点集合F_z中，如图12(b)所示。下一次生长将从新点p_j，开始采用以上方法重复执行，直到在生长方向上查找不到邻域点为止，如图12(c)所示。

(5)改变生长方向继续生长，直到算法结束。

选择初始种子点p生长方向的反方向-v_g进行折线生长操作，如图12(d)所示。

在折线生长的过程中，可能会出现在同一条特征线上多次设置种子点进行折线生长、重复收集这条边上的点的问题。为解决这一问题，每当从队列Q中取出一个种子点p时，就要相应的删除队列中点p的所有邻域点。执行折线生长算法会从队列Q中不断的取出种子点进行生长，直到队列为空，算法才结束。

(6)将集合F_z中的折线点连接成特征线。

步骤五：优化特征线

本发明对特征线进行优化，从而获得一条闭合且较为光顺的的特征线。优化特征线的具体过程为：

(1)连接端点接近、走向相符的特征线。

理想的预备体颈缘线是一条闭合的曲线，但限于预备体的形态和精度，在实际提取特征线过程中，部分预备体的属于颈缘线的特征线可能在曲率平缓的位置不连续，为此本发明通过连接端点接近、走向相符的特征线来使其连续且闭合。

本发明连接端点接近、走向相符的特征线的一个实施例为：对特征线的端点进行检测，若与之相邻特征线的端点为其邻域点，且两特征线夹角满足阈值，则连接两端点。

具体地，检测特征折线端点和其1-邻域顶点。如图13所示，对每个包含特征线的1-邻域顶点，如果满足α≤π/3，β≤π/3，γ≤π/2，则连接两个端点。其中α、β和γ为两结束段的夹角或结束段与连接两结束点直线的夹角。

(2)删除不闭合的和闭合但不属于颈缘线的特征线分支。

前面通过连接端点接近、走向相符的特征线使属于颈缘线的特征线闭合，但此时的特征线还包含一些不闭合的和闭合但不属于颈缘线的特征线分支，这些都都是无用的特征线，为此本发明对其逐一删除，从而得到一条闭合的属于颈缘线的特征线。

本发明删除不闭合的和闭合但不属于颈缘线的特征线分支的一个实施例为：

首先，遍历所有特征线，得到分支点；然后，分别以各分支点为起点，对其每个分支进行搜索，若在某一分支搜索到分支点则搜索另一分支，若在某一分支搜索到端点则删除这一分支；接着，重复S422直至删除所有不闭合分支；最后，对于闭合分支的分支点，通过判断分支点处各分支的走向来删除不属于颈缘线的分支。

(3)对特征线进行光顺处理。

为了得到较为理想的颈缘线，本发明还对特征线进行光顺处理，减少离散点的波动，以有利于最终拟合出光顺的颈缘线。

本发明对特征线进行光顺处理的一个实施例为：对局部过度弯折的特征线通过替换和翻转操作进行处理。其中，所述的替换和翻转操作如图14所示。处理时，优先采用替换操作，若没有适当的边可供替换时，则采用翻转操作。

步骤六：拟合颈缘线

本发明以优化后的特征线上的点作为控制点，拟合成颈缘线。

本发明拟合成颈缘线的一个实施例为：以优化后的特征线上点作为控制点，使用k次Bezier曲线b_j∈R^k+1拟合这些点得到颈缘线。

本发明提取的颈缘线效果图如图15所示。

与现有技术相比，本发明采用改进的Voronoi法降低了离散网格曲率计算算法的复杂度，同时又避免了狭长及不规则三角面片对顶点曲率的影响，提高了精度和鲁棒性；结合K-means聚类的方法和改进的折线生长技术，并经过特征线优化和曲线拟合，最终得到一条光顺的特征线。采用本发明的方法能够稳定且自动地提取出高精度的预备体颈缘线。

以上是对本发明的较佳实施进行了具体说明，但本发明创造并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做出种种的等同变形或替换，这些等同的变型或替换均包含在本申请权利要求所限定的范围内。

Claims (10)

1.一种预备体网格模型的颈缘线自动提取方法，其特征在于，包括：

S1、计算预备体网格模型顶点的平均曲率；

S2、采用基于K-means聚类的方法处理所得曲率值，得到特征点集F；

S3、采用一种改进的折线生长技术，得到折线点集合F_z，进而生成特征线；

S4、优化特征线；

S5、以优化后的特征线上的点作为控制点，拟合成颈缘线。

2.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S1采用一种改进的Voronoi法计算预备体网格模型顶点的平均曲率，包括：

S11、根据Laplace-Beltrami算子性质和高斯定理求顶点p_i的平均曲率向量，所述的平均曲率向量计算公式为：

$K\left(p_i\right)=\frac{1}{2A_M}\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}(\cot {\mathrm{\α}}_{\mathrm{ij}}+\cot {\mathrm{\β}}_{\mathrm{ij}})(p_i-p_j)$

其中，A_M为顶点p_i的1-邻域三角面片的Voronoi区域面积之和；N(i)为与顶点p_i直接相连的顶点索引的集合；p_j为与顶点p_i直接相连的顶点；α_ij＝∠p_ip_j-1p_j，β_ij＝∠p_ip_j+1p_j，α_ij、β_ij为边p_ip_j的两对角；

S12、求顶点p_i处的单位法矢，所述的单位法矢计算公式为：

$n_{\mathrm{pi}}=\frac{{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}}{\left|\right|{\mathrm{\Σ}}_{f_k\∈F_{\mathrm{pi}}}{A}_{\mathrm{fk}}{N}_{\mathrm{fk}}\left|\right|}$

其中，N_fk为三角面片f_k的法矢；A_fk表示三角面片f_k的Voronoi区域面积；

S13、根据Laplace-Beltrami算子K(p_i)＝2K_H(p_i)n_pi计算出顶点p_i处的平均曲率，所述的平均曲率的计算公式为：

$K_H\left(p_i\right)=\frac{1}{2}\left|\right|K\left(p_i\right)\left|\right|\mathrm{sgn}(K\left(p_i\right)\·n_{\mathrm{pi}})$

其中，sgn(·)为符号函数。

3.根据权利要求2所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S11中所述的顶点p_i的1-邻域任意三角面片PQR(P为顶点p_i)的Voronoi区域面积计算公式为：

其中，S_Δ为三角面片PQR的面积。

4.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S2，包括：

S21、初始化：选定平均曲率作为数据空间进行聚类，从中选择K个正曲率值作为初始聚类中心，并设定聚类中心变化阈值；

S22、对于数据空间中的平均曲率，根据它们与这些聚类中心的距离，按距离最近的准则分别将它们分配给与其最相似的聚类中心所代表的类；

S23、计算每个类别中所有曲率的均值作为该类别的新聚类中心，计算数据空间中所有曲率到其所在类别聚类中心的距离平方和，即J(C)值；

S24、得到新的聚类中心后，将依据新的聚类中心对待测数据进行重新聚类，并不断重复此过程，直到聚类中心的变化达到某一设定阈值为止；

S25、从聚类后的结果中选择其中曲率较高的几类作为特征点。

5.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S3，包括：

S31、将所得特征点排序存储至队列Q中；

S32、计算生长方向v_g；

S33、在v_g方向上进行折线生长，直到找到新的生长点；

S34、找新的生长点，将其添加到集合F_z中；

S35、改变生长方向继续生长，直到算法结束；

S36、将集合F_z中的折线点连接成特征线。

6.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S4，包括：

S41、连接端点接近、走向相符的特征线；

S42、删除不闭合的和闭合但不属于颈缘线的特征线分支；

S43、对特征线进行光顺处理。

7.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S41，包括：

对特征线的端点进行检测，若与之相邻特征线的端点为其邻域点，且两特征线夹角满足阈值，则连接两端点。

8.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S42，包括：

S421、遍历所有特征线，得到分支点；

S422、分别以各分支点为起点，对其每个分支进行搜索，若在某一分支搜索到分支点则搜索另一分支，若在某一分支搜索到端点则删除这一分支；

S423、重复S422直至删除所有不闭合分支；

S424、对于闭合分支的分支点，通过判断分支点处各分支的走向来删除不属于颈缘线的分支。

9.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S43，包括：

对局部过度弯折的特征线通过替换和翻转操作进行处理。

10.根据权利要求1所述的一种预备体网格模型的颈缘线自动提取方法，其特征在于，所述步骤S5，包括：

以优化后的特征线上点作为控制点，使用k次Bezier曲线b_j∈R^k+1拟合这些点得到颈缘线。