CN107689049B - 一种牙齿预备体修复模型特征线提取方法 - Google Patents
一种牙齿预备体修复模型特征线提取方法 Download PDFInfo
- Publication number
- CN107689049B CN107689049B CN201610634704.3A CN201610634704A CN107689049B CN 107689049 B CN107689049 B CN 107689049B CN 201610634704 A CN201610634704 A CN 201610634704A CN 107689049 B CN107689049 B CN 107689049B
- Authority
- CN
- China
- Prior art keywords
- point
- node
- path
- characteristic
- points
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000002360 preparation method Methods 0.000 title claims abstract description 66
- 238000000605 extraction Methods 0.000 title claims abstract description 43
- 238000005457 optimization Methods 0.000 claims abstract description 24
- 238000010845 search algorithm Methods 0.000 claims abstract description 10
- 238000000034 method Methods 0.000 claims description 56
- 239000013598 vector Substances 0.000 claims description 55
- 238000004364 calculation method Methods 0.000 claims description 20
- 239000004576 sand Substances 0.000 claims description 20
- 230000008569 process Effects 0.000 claims description 14
- 230000000694 effects Effects 0.000 claims description 10
- 230000000877 morphologic effect Effects 0.000 claims description 8
- 230000008859 change Effects 0.000 claims description 7
- 239000011159 matrix material Substances 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 4
- 238000009499 grossing Methods 0.000 claims description 4
- 238000004458 analytical method Methods 0.000 claims description 3
- 238000013461 design Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 4
- 238000010276 construction Methods 0.000 description 3
- 238000013507 mapping Methods 0.000 description 3
- 235000010575 Pueraria lobata Nutrition 0.000 description 2
- 241000219781 Pueraria montana var. lobata Species 0.000 description 2
- 210000003464 cuspid Anatomy 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 210000004283 incisor Anatomy 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- LTXREWYXXSTFRX-QGZVFWFLSA-N Linagliptin Chemical group N=1C=2N(C)C(=O)N(CC=3N=C4C=CC=CC4=C(C)N=3)C(=O)C=2N(CC#CC)C=1N1CCC[C@@H](N)C1 LTXREWYXXSTFRX-QGZVFWFLSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007797 corrosion Effects 0.000 description 1
- 238000005260 corrosion Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000010191 image analysis Methods 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000010355 oscillation Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0012—Biomedical image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20072—Graph-based image processing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30036—Dental; Teeth
Landscapes
- Engineering & Computer Science (AREA)
- Quality & Reliability (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Health & Medical Sciences (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Dental Tools And Instruments Or Auxiliary Dental Instruments (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种牙齿预备体修复模型特征线提取方法，包括以下步骤：在预备体的特征区域拾取关键特征点，然后根据综合约束控制函数采用启发式搜索算法搜索出关键特征点间的所有最优候选节点，最后顺序连接所有最优候选节点，以构成关键特征点间的初始特征路径，综合约束控制函数以距离、方向和特征作为综合约束条件；对关键特征点间的初始特征路径进行编辑；采用改进的特征路径光顺算法对编辑后的初始特征路径进行形态优化，得到牙齿预备体修复模型的特征线，改进的特征路径光顺算法直接根据预备体网格的空间拓扑形态在三维空间内完成初始特征路径的形态优化。本发明具有准确、稳定、效率高和提取质量有保证的优点，可广泛应用于义齿修复领域。
Description
技术领域
本发明涉及义齿修复领域，尤其是一种牙齿预备体修复模型特征线提取方法。
背景技术
名词解释：
三角网格：多边形网格的一种，是计算机图形学中用于为各种不规则物体建立模型的一种数据结构。现实世界中的物体表面直观上看都是由曲面构成的；而在计算机世界中，由于只能用离散的结构去模拟现实中连续的事物。所以现实世界中的曲面实际上在计算机里是由无数个小的多边形面片组成的。
牙齿预备体：在冠、嵌体的数字化制作过程中，口腔医生需要首先对患牙进行预先制备(简称预备)，而牙齿的预备过程是由牙科医生在患牙上，利用高速电钻根据医学知识，预备一个符合形态要求的牙齿预备体。
特征信息：是指能够表现模型显著形状的信息，主要包括脊线、谷线、拐点和边界线等要素。在微分几何中，一般通过主曲率、主方向以及主曲率极值系数等微分量来计算曲面的特征信息。
特征线：由包含特征信息的特征点连接所形成的一条线。
近年来，CAD/CAM技术在口腔修复领域中的应用取得了飞速的发展。相比于传统的手工修复方式，口腔CAD/CAM技术显著地降低了牙科技师的劳动强度，极大地提高了修复效率，具备更加优良的修复质量，同时满足了个性化修复体的设计需求。牙齿预备体的特征线包括冠预备体颈缘线和嵌体预备体洞缘线，是全冠修复以及嵌体修复的基准线，其提取质量直接影响修复体的建模精度，关系到患者佩戴修复体的舒适度，对牙龈健康具有重要影响。因此，牙齿预备体特征线的提取是数字化口腔修复体造型设计的关键环节。预备体一般具有复杂的几何形貌特征，且其特征质量取决于医生预备水平和测量所获得的网格质量，表面缺陷、噪音数据等都会对特征线的提取造成影响。目前网格曲面上的特征线提取方法主要分为自动提取方法和半自动提取方法这两大类。
1、自动提取方法。
邱彦杰等人在名为“基于Morse-Smale复形的三角网格特征线提取”的论文中提出了基于Morse-Smale(MS)复形的特征线提取算法，该方法通过计算网格顶点曲率来构造指标函数，并以此为依据建立MS复形，定义显著度作为判断特征线重要程度的控制参数，最后通过复形简化过程依次删除次要特征，从而实现相互连接的特征线的自动提取。该方法中特征线强度与显著度的设置也需通过试验来确定，对于具有特征干涉的模型无法提取出准确的特征线。
2、半自动提取方法。
刘胜兰等人在名为“用主动轮廓模型优化网格曲面上的特征线”的论文中提出了一种基于主动轮廓模型的特征线优化提取方法：首先交互地选取几个初始特征点，利用追踪投影法快速确定初始特征线；然后计算初始特征线上特征点的主动轮廓模型能量；最后，将初始特征线上的特征点经多次迭代后移动到能量极小处，实现优化。
葛闪等人在名为“一种新的网格曲面上的特征线提取方法”的论文中提出了一种与刘胜兰等人类似的特征线提取方法，但其并未用主动轮廓模型能量进行优化，而是找出特征采样点n环邻域内平均曲率最大的点作为新的特征点，再利用这些新的特征点拟合B样条曲线并将B样条曲线投影在三角网格曲面上。
刘胜兰等人和葛闪等人所提出的这两种方法都需要映射到二维平面，且仅使用了曲率信息，只适合小范围的形态优化，导致人工选取的两个特征点间的形态不能发生较大变化，难以保证特征线的提取质量；需要选取较多的特征点才能提取出一条准确的特征线，效率较低。
综上所述，目前网格曲面上的特征线提取方法主要存在以下缺点：
(1)特征线自动提取方法在局部特征干涉区域无法实现牙齿预备体特征线的准确提取，且易受网格模型质量的影响，不够稳定。
(2)特征线半自动提取方法需人工选取多个特征点才能提取出一条准确的特征线，且两个人工选取的特征点间网格模型的形态不能发生剧烈变化。
发明内容
为解决上述技术问题，本发明的目的在于：提供一种准确、稳定、效率高和提取质量有保证的牙齿预备体修复模型特征线提取方法。
本发明所采取的技术方案是：
一种牙齿预备体修复模型特征线提取方法，包括以下步骤：
S1、在预备体的特征区域拾取关键特征点，然后根据综合约束控制函数采用启发式搜索算法搜索出关键特征点间的所有最优候选节点，最后顺序连接所有最优候选节点，以构成关键特征点间的初始特征路径，所述综合约束控制函数以距离、方向和特征作为综合约束条件；
S2、对关键特征点间的初始特征路径进行编辑，以修正初始特征路径的偏差；
S3、采用改进的特征路径光顺算法对编辑后的初始特征路径进行形态优化，得到牙齿预备体修复模型的特征线，所述改进的特征路径光顺算法直接根据预备体网格的空间拓扑形态在三维空间内完成初始特征路径的形态优化。
进一步，所述步骤S1包括：
S11、在牙齿预备体的特征区域拾取特征路径的起始点Vs和末端点Ve，将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中，其中，Vs的代价信息根据综合约束控制函数计算得到；
S12、取出NodeMap的第一个元素Vtop，并判断Vtop＝Ve是否成立，若是，则结束搜索过程，转到步骤S14；反之，则执行步骤S13；
S13、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空，然后遍历Neibourhood中的每一个点，将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中，最后返回步骤S12，其中，o为Vtop的一环邻域点总数；
S14、顺次连接起始点Vs和末端点Ve之间的所有路径节点，构造出初始特征路径。
进一步，所述步骤S11包括：
在牙齿预备体的特征区域标记特征路径可能的起始点和末端点；
判断标记的起始点和末端点是否为网格顶点，若是，则以标记的起始点和末端点作为拾取的起始点Vs和末端点Ve，反之，则从标记的起始点或末端点所在网格边的两端点和所在三角片的三个顶点中选择一个主曲率值最高的点作为拾取的起始点Vs或末端点Ve；
根据综合约束控制函数计算起始点Vs的代价信息Fs，所述综合约束控制函数f(g)的表达式为：f(g)＝α*fdis+β*fdir+γ*(fc+ft)，Fs＝f(s)，其中，f(g)为当前候选点Vg的综合约束控制函数，fdis为当前候选点Vg与末端点Ve之间的欧式几何距离，fdir为当前候选点Vg与当前搜索中心点Vc构建的当前路径方向相对于前一步路径方向的角度变化值；fc为当前候选点Vg最大主曲率的倒数，ft为当前候选点Vg相对于当前搜索中心点Vc最小主方向Tc的变化值，α、β和γ分别为距离控制因子、方向控制因子和特征控制因子；
将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中。
进一步，所述综合约束控制函数f(g)中fc和ft的计算过程为：
估算当前点p的法矢Ni，所述当前点p的法矢Ni的估算公式为：
其中，当前点p为当前候选点Vg或当前搜索中心点Vc，lj和nj分别表示当前点p在一环邻域内第j个三角片的对边和法矢，m表示当前点p在一环邻域内网格顶点的个数；
根据当前点p的法矢Ni在网格顶点处建立局部坐标系p-uvw，所述局部坐标系p-uvw由全局坐标系O-XYZ执行以下操作得到：将O点平移到p点，将Z轴旋转至与法矢Ni重合，将X和Y轴分别取为u和v；
在局部坐标系p-uvw中根据网格顶点的2环邻域点确定相应的局部三次曲面方程f(x,y)，然后通过最小二乘拟合方法求解f(x,y)的系数值，进而确定相应的Weingarten曲率矩阵W1，所述局部三次曲面方程f(x,y)的表达式为：
其中，A、B、C、D、E、F和G均为f(x,y)的系数值，
对Weingarten曲率矩阵W1进行奇异值分解，得到与当前网格顶点的最大主曲率和最小主曲率相对应的两个特征值，然后根据这两个特征值对应的两个特征向量t1和t2确定当前网格顶点的最大主方向tmax和最小主方向tmin，其中，t1、t2、tmax和tmin的表达式分别为：t1＝(t11,t12)T，t2＝(t21,t22)T，tmax＝t11x+t12y，tmin＝t21x+t22y；
根据t1、t2、tmax和tmin确定综合约束控制函数f(g)中fc和ft的值。
进一步，所述步骤S13包括：
S131、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空；
S132、以Neibourhood中的任意一点Vn作为当前候选点进行遍历分析，最后将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中；
S133、判断Neibourhood中所有点的遍历分析是否都已完成，若是，则返回步骤S12，反之，则返回步骤S132。
进一步，所述步骤S132包括：
根据综合约束控制函数f(g)计算当前候选点Vn的代价信息Cost；
判断当前候选点Vn是否已在候选队列NodeQueue中，若是，则表明当前候选点Vn之前已作为候选点被搜索过，此时对该点不做任何处理，继续选取Neibourhood中的下一个候选点来计算其代价信息；反之，则将Vn加入到候选队列NodeQueue中，并令Fn＝Cost，最后将Vn以及Fn存入关联容器NodeMap中。
进一步，所述步骤S2具体为：
判断关键特征点间的初始特征路径是否与期望的特征线存在偏差，若是，则采用点编辑操作对关键特征点间的初始特征路径进行编辑，反之，则执行步骤S3，所述点编辑操作包括但不限于删除点操作、***点操作和移动点操作。
进一步，所述步骤S3包括：
S31、遍历编辑后的初始特征路径path{V1,V2,…,Vh}中的所有节点，计算path中每一个节点Vi对应的关联点集Connect(Vi-1,Vi,Vi+1)，并构造Vi到前一节点Vi-1的方向向量Dir1以及Vi到后一节点Vi+1的方向向量Dir2，其中，h为path的节点总数；
S32、计算Dir1与Dir2间的夹角Angle(Vi)大小；
S33、根据夹角Angle(Vi)大小判断是否需要在节点Vi-1和Vi+1之间构建新的局部路径进行局部路径的光顺，若是，则在节点Vi-1和Vi+1之间构建新的局部路径；反之，则不在节点Vi-1和Vi+1之间构建新的局部路径。
进一步，所述步骤S33中：
S331、根据夹角Angle(Vi)大小计算节点Vi的光顺权值P(Vi)，并将节点Vi的光顺权值P(Vi)存入光顺权值集合W中，所述节点Vi的光顺权值P(Vi)的计算公式为：
S332、取出权值集合W中的最大权值wmax，并得到wmax对应的节点Vi及关联点集Connect(Vi-1,Vi,Vi+1)；
S333、根据wmax对应的节点Vi在节点Vi-1和Vi+1之间构建新的局部路径；
S334、继续从余下的权值集合W中取出最大权值wmax，然后返回步骤S333，直到权值集合W为空，从而完成初始特征路径的形态优化。
进一步，所述步骤S333包括：
S3331、判断节点Vi属于网格顶点还是属于网格边上的点，若属于网格顶点，则执行步骤S3332，反之，则执行步骤S3333；
S3332、遍历节点Vi的一环邻域点集合Ei{X1,X2,…,Xq}中的每一点，计算向量ViXi与向量Vi-1Xi间的夹角和向量ViXi与向量Vi+1Xi间的夹角这两个夹角之和θi，然后判断θi是否小于π，若是，则先得到节点Vi相交边ViXi的集合Ci，然后构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，接着，计算平面Pane与集合Ci中相交边的交点，最后，根据计算的交点构建新的局部路径：若平面Pane与集合Ci中的部分相交边有交点，则计算集合Ci中与Pane无交点的相交边非节点端点以及节点端点Vi到平面Pane的距离大小，将距离平面Pane最近的端点作为这些相交边与Pane的交点，并将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；若平面Pane与集合Ci中的全部相交边均有交点，则将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；反之，则直接将节点Vi-1和Vi+1连接起来作为新的局部路径，其中，q为Ei的点的总数，Xi为Vi对应的的一环邻域点，集合Ci内的元素按照θi从小到大进行排序，给定的平面法矢
的值为相交边所在三角片法矢的加权平均；
S3333、以节点Vi所在的边作为相交边，构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，然后计算平面Pane与相交边的交点，若平面Pane与相交边的交点不存在，则计算节点端点Vi以及节点Vi所在边的两个端点到平面P的距离大小，将距离平面最近的端点作为交点，所述给定的平面法矢
的值为相交边所在三角片法矢的加权平均；接着，将节点Vi从特征路径Path中删除，并在当前位置将交点顺序加入到特征路径Path中，最后，对与节点Vi相关的节点及其权值信息进行相应的更新：首先从权值集合W中删除节点Vi的关联点集Connect(Vi-1,Vi,Vi+1)中各节点的光顺权值，然后更新节点Vi的关联点集中两相邻节点Vi-1和Vi+1对应的关联点集及光顺权值，最后将更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值存入权值集合W中，其中，更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值计算公式为：
其中，delete为删除操作，updadte为更新操作，Vi+2为Vi+1的后一节点。
本发明的有益效果是：基于启发式特征线提取技术，只需在预备体的特征区域拾取关键特征点，然后以距离、方向以及特征作为综合约束设计综合约束控制函数，结合启发式搜索算法提取出关键特征点间初始特征线，就能够实现存在局部特征干涉区域的预备体模型特征线的提取，可以更好的应对各种质量状况的预备体模型，提高了特征线提取的准确性和稳定性；采用了启发式搜索算法，无需大量拾取关键特征点便可实现预备体特征线的完整提取，效率更高；采用改进的特征路径光顺算法对编辑后的初始特征路径进行形态优化，根据网格空间拓扑形态在三维空间内对初始特征线进行形态优化，无需映射到二维平面，保证了特征线的提取质量，使得最终提取的特征线形态准确、光顺自然。
附图说明
图1为本发明一种牙齿预备体修复模型特征线提取方法的整体流程图；
图2为本发明综合约束控制函数的示意图；
图3为本发明局部曲面拟合计算离散微分量时建立的局部坐标系示意图；
图4为本发明特征路径形态优化的四种情况示意图；
图5为本发明冠预备体模型提取出的特征线优化前后的形态对比图；
图6为切牙预备体采用本发明的方法进行颈缘线提取的效果图；
图7为尖牙预备体采用本发明的方法进行颈缘线提取的效果图；
图8为磨牙预备体采用本发明的方法进行颈缘线提取的效果图；
图9为三面嵌体预备体采用本发明的方法进行洞缘线提取的效果图。
具体实施方式
参照图1，一种牙齿预备体修复模型特征线提取方法，包括以下步骤：
S1、在预备体的特征区域拾取关键特征点，然后根据综合约束控制函数采用启发式搜索算法搜索出关键特征点间的所有最优候选节点，最后顺序连接所有最优候选节点，以构成关键特征点间的初始特征路径，所述综合约束控制函数以距离、方向和特征作为综合约束条件；
S2、对关键特征点间的初始特征路径进行编辑，以修正初始特征路径的偏差；
S3、采用改进的特征路径光顺算法对编辑后的初始特征路径进行形态优化，得到牙齿预备体修复模型的特征线，所述改进的特征路径光顺算法直接根据预备体网格的空间拓扑形态在三维空间内完成初始特征路径的形态优化。
进一步作为优选的实施方式，所述步骤S1包括：
S11、在牙齿预备体的特征区域拾取特征路径的起始点Vs和末端点Ve，将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中，其中，Vs的代价信息根据综合约束控制函数计算得到；
S12、取出NodeMap的第一个元素Vtop，并判断Vtop＝Ve是否成立，若是，则结束搜索过程，转到步骤S14；反之，则执行步骤S13；
S13、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空，然后遍历Neibourhood中的每一个点，将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中，最后返回步骤S12，其中，o为Vtop的一环邻域点总数；
S14、顺次连接起始点Vs和末端点Ve之间的所有路径节点，构造出初始特征路径。
进一步作为优选的实施方式，所述步骤S11包括：
在牙齿预备体的特征区域标记特征路径可能的起始点和末端点；
判断标记的起始点和末端点是否为网格顶点，若是，则以标记的起始点和末端点作为拾取的起始点Vs和末端点Ve，反之，则从标记的起始点或末端点所在网格边的两端点和所在三角片的三个顶点中选择一个主曲率值最高的点作为拾取的起始点Vs或末端点Ve；
根据综合约束控制函数计算起始点Vs的代价信息Fs，所述综合约束控制函数f(g)的表达式为：f(g)＝α*fdis+β*fdir+γ*(fc+ft)，Fs＝f(s)，其中，f(g)为当前候选点Vg的综合约束控制函数，fdis为当前候选点Vg与末端点Ve之间的欧式几何距离，fdir为当前候选点Vg与当前搜索中心点Vc构建的当前路径方向相对于前一步路径方向的角度变化值；fc为当前候选点Vg最大主曲率的倒数，ft为当前候选点Vg相对于当前搜索中心点Vc最小主方向Tc的变化值，α、β和γ分别为距离控制因子、方向控制因子和特征控制因子；
将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中。
进一步作为优选的实施方式，所述综合约束控制函数f(g)中fc和ft的计算过程为：
估算当前点p的法矢Ni，所述当前点p的法矢Ni的估算公式为：
其中，当前点p为当前候选点Vg或当前搜索中心点Vc，lj和nj分别表示当前点p在一环邻域内第j个三角片的对边和法矢，m表示当前点p在一环邻域内网格顶点的个数；
根据当前点p的法矢Ni在网格顶点处建立局部坐标系p-uvw，所述局部坐标系p-uvw由全局坐标系O-XYZ执行以下操作得到：将O点平移到p点，将Z轴旋转至与法矢Ni重合，将X和Y轴分别取为u和v；
在局部坐标系p-uvw中根据网格顶点的2环邻域点确定相应的局部三次曲面方程f(x,y)，然后通过最小二乘拟合方法求解f(x,y)的系数值，进而确定相应的Weingarten曲率矩阵W1，所述局部三次曲面方程f(x,y)的表达式为：
其中，A、B、C、D、E、F和G均为f(x,y)的系数值，
对Weingarten曲率矩阵W1进行奇异值分解，得到与当前网格顶点的最大主曲率和最小主曲率相对应的两个特征值，然后根据这两个特征值对应的两个特征向量t1和t2确定当前网格顶点的最大主方向tmax和最小主方向tmin，其中，t1、t2、tmax和tmin的表达式分别为：t1＝(t11,t12)T，t2＝(t21,t22)T，tmax＝t11x+t12y，tmin＝t21x+t22y；
根据t1、t2、tmax和tmin确定综合约束控制函数f(g)中fc和ft的值。
进一步作为优选的实施方式，所述步骤S13包括：
S131、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空；
S132、以Neibourhood中的任意一点Vn作为当前候选点进行遍历分析，最后将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中；
S133、判断Neibourhood中所有点的遍历分析是否都已完成，若是，则返回步骤S12，反之，则返回步骤S132。
进一步作为优选的实施方式，所述步骤S132包括：
根据综合约束控制函数f(g)计算当前候选点Vn的代价信息Cost；
判断当前候选点Vn是否已在候选队列NodeQueue中，若是，则表明当前候选点Vn之前已作为候选点被搜索过，此时对该点不做任何处理，继续选取Neibourhood中的下一个候选点来计算其代价信息；反之，则将Vn加入到候选队列NodeQueue中，并令Fn＝Cost，最后将Vn以及Fn存入关联容器NodeMap中。
进一步作为优选的实施方式，所述步骤S2具体为：
判断关键特征点间的初始特征路径是否与期望的特征线存在偏差，若是，则采用点编辑操作对关键特征点间的初始特征路径进行编辑，反之，则执行步骤S3，所述点编辑操作包括但不限于删除点操作、***点操作和移动点操作。
进一步作为优选的实施方式，所述步骤S3包括：
S31、遍历编辑后的初始特征路径path{V1,V2,…,Vh}中的所有节点，计算path中每一个节点Vi对应的关联点集Connect(Vi-1,Vi,Vi+1)，并构造Vi到前一节点Vi-1的方向向量Dir1以及Vi到后一节点Vi+1的方向向量Dir2，其中，h为path的节点总数；
S32、计算Dir1与Dir2间的夹角Angle(Vi)大小；
S33、根据夹角Angle(Vi)大小判断是否需要在节点Vi-1和Vi+1之间构建新的局部路径进行局部路径的光顺，若是，则在节点Vi-1和Vi+1之间构建新的局部路径；反之，则不在节点Vi-1和Vi+1之间构建新的局部路径。
进一步作为优选的实施方式，所述步骤S33中：
S331、根据夹角Angle(Vi)大小计算节点Vi的光顺权值P(Vi)，并将节点Vi的光顺权值P(Vi)存入光顺权值集合W中，所述节点Vi的光顺权值P(Vi)的计算公式为：
S332、取出权值集合W中的最大权值wmax，并得到wmax对应的节点Vi及关联点集Connect(Vi-1,Vi,Vi+1)；
S333、根据wmax对应的节点Vi在节点Vi-1和Vi+1之间构建新的局部路径；
S334、继续从余下的权值集合W中取出最大权值wmax，然后返回步骤S333，直到权值集合W为空，从而完成初始特征路径的形态优化。
进一步作为优选的实施方式，所述步骤S333包括：
S3331、判断节点Vi属于网格顶点还是属于网格边上的点，若属于网格顶点，则执行步骤S3332，反之，则执行步骤S3333；
S3332、遍历节点Vi的一环邻域点集合Ei{X1,X2,…,Xq}中的每一点，计算向量ViXi与向量Vi-1Xi间的夹角和向量ViXi与向量Vi+1Xi间的夹角这两个夹角之和θi，然后判断θi是否小于π，若是，则先得到节点Vi相交边ViXi的集合Ci，然后构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，接着，计算平面Pane与集合Ci中相交边的交点，最后，根据计算的交点构建新的局部路径：若平面Pane与集合Ci中的部分相交边有交点，则计算集合Ci中与Pane无交点的相交边非节点端点以及节点端点Vi到平面Pane的距离大小，将距离平面Pane最近的端点作为这些相交边与Pane的交点，并将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；若平面Pane与集合Ci中的全部相交边均有交点，则将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；反之，则直接将节点Vi-1和Vi+1连接起来作为新的局部路径，其中，q为Ei的点的总数，Xi为Vi对应的的一环邻域点，集合Ci内的元素按照θi从小到大进行排序，给定的平面法矢
的值为相交边所在三角片法矢的加权平均；
S3333、以节点Vi所在的边作为相交边，构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，然后计算平面Pane与相交边的交点，若平面Pane与相交边的交点不存在，则计算节点端点Vi以及节点Vi所在边的两个端点到平面P的距离大小，将距离平面最近的端点作为交点，所述给定的平面法矢
的值为相交边所在三角片法矢的加权平均；接着，将节点Vi从特征路径Path中删除，并在当前位置将交点顺序加入到特征路径Path中，最后，对与节点Vi相关的节点及其权值信息进行相应的更新：首先从权值集合W中删除节点Vi的关联点集Connect(Vi-1,Vi,Vi+1)中各节点的光顺权值，然后更新节点Vi的关联点集中两相邻节点Vi-1和Vi+1对应的关联点集及光顺权值，最后将更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值存入权值集合W中，其中，更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值计算公式为：
其中，delete为删除操作，updadte为更新操作，Vi+2为Vi+1的后一节点。
下面结合说明书附图和具体实施例对本发明作进一步解释和说明。
实施例一
针对现有特征线提取技术的缺点，本发明提出了一种面向口腔医生的启发式牙齿预备体模型特征线提取方法，只需人工选取3～5个特征点便能够实现各类常见预备体模型特征线的提取，并在局部特征干涉区域仍然可以实现特征线的提取，还采用了新的形态优化算法使提取的特征线形态更加准确和光顺自然。
如图1所示，本发明的牙齿预备体修复模型特征线提取方法包括以下步骤：
步骤一、在预备体的颈缘或洞缘特征区域交互拾取关键特征点，然后进行相邻两个关键特征点间的初始特征路径提取。
本发明任意两相邻关键特征点间的路径节点都是通过启发式搜索算法得到的。启发式搜索算法的原理为：搜索过程从起始点出发，利用控制函数评估每一个候选点的代价，通过比较代价的大小决定下一个搜索点的位置并跳转到其上，再对当前搜索位置做进一步的搜索，渐进地到达最终的目标点。而起始点与目标点之间的路径由一系列中间节点顺次连接而成，所有中间节点都是利用控制函数比较代价大小后得到的最优候选点。
两相邻关键特征点间的特征路径提取问题可以转化为：构造三角网格曲面上两个顶点之间的特征路径。用Vs表示搜索起始点，Ve表示搜索目标点，Vc表示当前搜索中心点。首次搜索时，以Vs作为搜索中心点。而在搜索过程中，评估搜索中心点的一环邻域点集中每个候选点通过控制函数计算得到的代价大小。以当前搜索中心点Vc搜索代价最小的候选点Vn作为下一步新的搜索中心点，搜索中心点Vc随搜索过程的进行不断更新。当搜索中心点Vc为搜索目标点Ve时，结束整个搜索过程。最优候选点决定了路径方向，而如何准确的筛选最优候选点，取决于控制函数f(g)的设计，所以该问题的核心为控制函数的设计。
本发明结合当前候选点与探索中心点、探索目标点的位置信息、特征信息设计了距离控制函数fdis、方向控制函数fdir、特征控制函数fc和ft，三个控制函数的具体定义如下：
fdis＝||vn,ve|| (1)
fdir＝θdir (2)
fdis：表示当前候选点与目标点之间的欧式几何距离。该函数确保了搜索最终能达到目标点，保证了搜索规则满足收敛性，避免出现振荡发散的死循环。fdis的值越小，表示特征路径的搜索越来越接近于最终目标点，该参数越优。
fdir：表示当前候选点与搜索中心点构建的当前路径方向相对于前一步路径方向的角度变化，如图2中的角度θdir所示。该函数可以保持搜索路径的光顺程度，减少锯齿状路径的产生。θdir越小，表示当前路径方向越接近于前一步路径方向，特征路径节点间的变化比较平缓；此外，当θdir的值大于90°时，表明当前候选路径方向与路径整体方向相反，此时需直接将此类候选点从候选点集中删除，以加快搜索效率，并防止回溯探索。
fc：表示当前候选点最大主曲率K的倒数。ft：表示当前候选点的最小主方向Tn相对于当前搜索中心点最小主方向Tc的变化，如图2中的角度θt所示。特征控制函数确保搜索到的路径能够最大程度地描述特征区域的显著特征特。(fc+ft)越小，该参数越优。
三角网格曲面是一种分片线性曲面，无法用参数方程或隐式方程精确表示，一般可通过曲面拟合的方法逼近网格形态，用拟合曲面的微分几何特征近似表示三角网格的特征信息。本实施例采用了基于局部曲面拟合的方法计算预备体网格模型的离散微分量主曲率、主方向，具体过程如下：
(1)在网格顶点处建立局部坐标系。
局部拟合方法需在每个网格顶点建立局部坐标系p-uvw，然后在局部坐标系中用多项式曲面方程f(x,y)拟合当前点p及其邻域点，利用多项式曲面计算离散网格顶点的微分量，如图3所示。本实施例的局部坐标系可以通过如下方法获得：
1)估算p点的法矢Ni，这里Ni的计算方法采用下式的“对边平衡”的简化方法：
其中，lj、ni分别表示当前点一环邻域内第j三角片的对边和法矢；m表示当前点一环邻域内网格顶点的个数。
2)令全局坐标系O-XYZ执行如下操作来建立局部坐标系p-uvw：将O点平移到p点，将Z轴旋转至与法矢Ni重合，将X和Y轴分别取为u和v。
(2)在网格顶点处定义局部三次曲面方程
为了更为精确的计算网格顶点的特征信息，本发明采用了当前顶点的2环邻域点来拟合三次多项式曲面，其采用的三次多项式曲面方程f(x,y)可表示为：
f(x,y)对应的Weingarten曲率矩阵W1可表示为：
为了求解W1中的未知项，可将局部范围内的数据点{V1,V2,…,Vl}(当前数据点周围的二环邻域数据点)变换到当前点的局部坐标系下，然后利用最小二乘拟合方法求解f(x,y)的系数值，有：
X＝(A,B,C,D,E)T
其中，{nx,ny,nz}为该曲面上任意一点{xi,yi,zi}的法矢。
(3)计算网格顶点的主曲率与主方向。
将X中的系数A、B、C代入W1中，对W1进行奇异值分解可以得到两个特征值，以及与两个特征值对应的两个特征向量t1、t2。其中，最大的特征值对应于当前网格顶点的最大主曲率以及特征向量t1，最小特征值对应于当前网格顶点的最小主曲率以及特征向量t2。故当前网格顶点的最大、最小主方向tmax、tmin采用如下方法计算：
则综合约束控制函数f(g)可表示为：
f(g)＝α*fdis+β*fdir+γ*(fc+ft) (10)
α、β、γ表示三个控制函数的权因子，由于fdis、fdir、(fc+ft)三个控制函数分别代表长度、角度、曲率三个物理量，表征了不同的数量意义，所以确定这三个控制函数之间的权因子是一个很重要的问题。为保证在同一个数量级上，可分别对三个控制函数做归一化处理，归一化的对象范围为当前探索中心点的一环邻域点集，此时这三个控制函数的权因子均为1。归一化后的三个控制函数可以改写为以下形式：
式(11)、(12)和(13)中的最大值、最小值表示当前搜索中心点的一环邻域点集中相关数据(包括欧式几何距离d、角度变化dir、最大主曲率的倒数c和最小主方向变化t)的最大值、最小值。
基于上述综合约束控制函数的介绍，结合启发式搜索算法，本发明初始特征路径的提取步骤如下：
(1)在牙齿预备体的特征区域拾取特征路径起始点Vs、末端点Ve，将Vs以及Vs的代价信息Fs存入关联容器NodeMap(其内部已按照节点的代价大小进行排序)中，代价信息Fs由公式(10)计算得到，并将Vs加入候选队列NodeQueue中。
Vs、Ve均位于三角网格曲面上，实际提取过程中，拾取得到的标记点可能是网格顶点，也有可能是网格边上的点或者三角片内部的点。根据本发明启发式搜索算法的需求，标记点应该属于网格顶点，当标记点位于网格边上或者三角片内部时，可以从该点所在网格边的两端点和所在三角片的三个顶点中选择一个主曲率值最高的点作为间接搜索始、末点。
(2)取出NodeMap的第一个元素Vtop。若Vtop＝Ve，则结束搜索过程，转到步骤(4)；否则以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空，令Vn表示Neibourhood中的任意一点，由式(10)计算Vn的代价信息Cost，接着对Neibourhood中的每一个点Vn，均作如下的分析：
(2.1)若Vn已在候选队列NodeQueue中，则表明Vn之前已作为候选点被检索过，此时对此点不做任何处理，继续评估Neibourhood中下一个候选点的代价情况；
(2.2)若Vn不在候选队列NodeQueue中，则将Vn加入到队列NodeQueue中，并Fn＝Cost，最后将Vn以及Vn的代价信息Fn存入关联容器NodeMap中；
(3)判断对Neibourhood中的所有点遍历是否评估完成，若是，则执行步骤(4)，反之，则返回执行步骤(2)；
(4)完成特征路径的搜索，顺次连接目标点到末端点之间的所有路径节点，以构造出初始特征路径。
步骤二、特征路径编辑：为初步提取的特征路径中局部路径的偏差部分加入点编辑操作进行修正。
在所有关键特征点的两相邻关键特征点间分别生成特征路径可得到一条完整的特征路径。考虑到部分牙齿预备体的特线较长、形态变化较大的情况，尤其在局部特征干涉区域中，初步拾取的关键特征点可能无法提取出一条准确的特征线，这时就需要采用点编辑的交互方式来修正特征线的局部偏差。当提取出的特征线与操作者所期望的特征线存在偏差时，可以加入相对应的点编辑操作来修正特征线。其中，点编辑操作共有三种类型：删除点操作、***点操作和移动点操作。
步骤三、特征线路径形态优化。
初步提取到的特征路径由有序相连的网格边构成，故其形态上缺乏良好的光滑效果并极有可能存在锯齿状路径。为满足牙齿预备体特征线的光顺要求，本发明设计了新的特征路径光顺算法来优化特征路径的形态，使得光顺后的特征路径完全贴合在网格曲面上。
特征路径用Path{V1,V2,…,Vh}表示，将Path中任意一节点Vi和其顺次相连的前一个节点Vi-1、后一个节点Vi+1，定义为节点Vi的关联点集Connect(Vi-1,Vi,Vi+1)，则有每个节点与其关联点集的一一对应关系：Vi→Connect(Vi-1,Vi,Vi+1)。构造节点Vi到节点Vi-1、Vi+1的方向向量Dir1、Dir2，计算Dir1与Dir2的夹角Angle(Vi)，可以通过Angle(Vi)的大小判断是否需要在节点Vi-1与Vi+1之间构造新的路径进行局部路径光顺。所以新的特征路径光顺算法在构建光顺权值函数时，只需要考虑当前节点与前后节点的夹角因素Angle(Vi)即可。
定义光顺权值函数P的计算公式如下：
ε为设定的角度参数，控制整体的光顺强度。ε越小，光顺效果越好，但相对应的效率越低，故ε一般取值为175/180，可以得到很好的光顺处理效果。
基于上述理论基础，本发明的特征路径形态优化的步骤如下：
(1)遍历特征路径Path中的所有节点，计算每一个节点Vi对应的关联点集Connect(Vi-1,Vi,Vi+1)，并将节点Vi代入式(14)，计算权值P(Vi)，然后将权值P(Vi)存入权值集合W中；
(2)取出权值集合W中的最大值wmax，得到其对应的节点Vi及关联点集Connect(Vi-1,Vi,Vi+1)，接下来需要将wmax对应的节点Vi移动到合适的位置以完成局部路径的光顺，Vi新位置的寻找方法如下：
(2.1)当节点Vi属于网格顶点时，如图4(a)所示，将节点Vi的一环邻域点集合记为Ei{X1,X2,…,Xq}，遍历Ei中的每一点，计算向量ViXi与Vi-1Xi、Vi+1Xi,的夹角之和θi。如果θi小于π，就将边ViXi称为相交边，得到节点Vi的相交边集合Ci，相交边集合内的元素是按照向量ViXi与Vi-1Xi的夹角从小到大排序的。然后构建经过节点Vi-1和Vi+1并平行于给定的平面法矢Ni的平面Pane；计算平面Pane与相交边的交点，如果Ci中存在与平面Pane无交点的部分相交边，如图4(b)所示，则计算这些相交边非节点端点以及节点端点Vi到平面Pane的距离大小，将距离平面最近的端点作为交点。这里会存在一种特殊的情况，相交边的集合Ci为空，如图4(c)所示，此时不存在交点。若平面Pane与集合Ci中的部分相交边有交点，则将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径，新的局部路径如图4(b)的灰色粗黑线所示；若平面Pane与集合Ci中的全部相交边均有交点，则将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径，新的局部路径如图4(a)的灰色粗黑线所示；若平面Pane与集合Ci中的全部相交边均无交点，则直接将节点Vi-1和Vi+1连接起来作为新的局部路径，新的局部路径如图4(c)的灰色粗黑线所示。
对于平面Pane的构造，
的值为相交边所在三角片法矢的加权平均。Ni能反映局部特征路径的拓扑形态，保持局部特征路径的形态变化趋势；Vi-1和Vi+1可以保证局部特征路径光顺要求，使节点Vi与其前后节点的夹角大于所设定的阈值。
(2.2)当节点Vi位于网格边上时，如图4(d)所示，以节点Vi所在的边作为相交边，平面Pane的构建方法与步骤(2.1)的相同，计算平面Pane与相交边的交点，如果交点不存在，则计算节点端点Vi以及节点Vi所在边的两个端点到平面Pane的距离大小，将距离平面最近的端点作为交点。接着，将节点Vi从特征路径Path中删除，并在当前位置将交点顺序加入到特征路径Path中，这样与节点Vi相关的节点及其权值信息就需要做相应的更新：首先从权值集合W中删除节点Vi的关联点集中各节点的权值，其次更新节点Vi的关联点集中两相邻节点Vi-1和Vi+1对应的关联点集及权值，同时将更新后的关联点集的权值存入权值集合W中，更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值计算公式为：
(3)迭代执行步骤(2)，直到权值集合W为空，从而完成特征路径的光顺处理。
为了防止极端扭曲特征路径可能产生的无穷迭代趋势，可对特征路径上的每个节点附加迭代标记，以记录节点移动的次数，避免出现重复访问的情况。
按照上述算法，冠预备体模型提取出的特征线在优化前后的形态图分别如图5(a)和如图5(b)所示，可以看出形态优化后的特征线比优化前更加光顺和自然。
而采用本发明的方法分别对切牙、尖牙、磨牙和三面嵌体预备体模型特征线提取最终效果如图6-9所示，可以看出本发明提取的特征线十分光顺和自然，且适合于不同质量的牙齿模型的提取。
本发明提出了一种面向口腔医生的启发式牙齿预备体模型特征线提取方法，具有以下优点：
(1)采用基于局部三次多项式拟合曲面的方法计算预备体网格模型顶点的主曲率和主方向，并以最大主曲率和最小主方向设计特征控制函数，使得预备体模型的特征信息能够得到准确的表达。
(2)根据当前候选节点和关键特征点的相关信息设计了距离控制函数、方向控制函数和特征控制函数，通过控制函数计算出最优候选节点，并结合启发式搜索策略计算两个关键特征点间的所有最优候选节点，能够实现存在局部特征干涉区域的预备体模型特征线的提取，无需大量拾取关键特征点，只需3～5关键特征点即可实现特征线的完整提取，提取的特征线形态更加准确。
(3)直接根据预备体网格空间拓扑形态在三维空间内实现初始特征路径的形态优化，无需映射到二维平面，相比现有“构造平面—投影到平面—平面内光顺—映射回网格”的特征线优化方式，计算过程直接、误差小，保证了特征线的提取质量，使得最终提取的特征线形态准确、光顺自然。
以上是对本发明的较佳实施进行了具体说明，但本发明并不限于所述实施例，熟悉本领域的技术人员在不违背本发明精神的前提下还可做作出种种的等同变形或替换，这些等同的变形或替换均包含在本申请权利要求所限定的范围内。
Claims (7)
1.一种牙齿预备体修复模型特征线提取方法，其特征在于：包括以下步骤：
S1、在预备体的特征区域拾取关键特征点，然后根据综合约束控制函数采用启发式搜索算法搜索出关键特征点间的所有最优候选节点，最后顺序连接所有最优候选节点，以构成关键特征点间的初始特征路径，所述综合约束控制函数以距离、方向和特征作为综合约束条件；
S2、对关键特征点间的初始特征路径进行编辑，以修正初始特征路径的偏差；
S3、采用改进的特征路径光顺算法对编辑后的初始特征路径进行形态优化，得到牙齿预备体修复模型的特征线，所述改进的特征路径光顺算法直接根据预备体网格的空间拓扑形态在三维空间内完成初始特征路径的形态优化；
所述步骤S3包括：
S31、遍历编辑后的初始特征路径path{V1,V2,…,Vh}中的所有节点，计算path中每一个节点Vi对应的关联点集Connect(Vi-1,Vi,Vi+1)，并构造Vi到前一节点Vi-1的方向向量Dir1以及Vi到后一节点Vi+1的方向向量Dir2，其中，h为path的节点总数；
S32、计算Dir1与Dir2间的夹角Angle(Vi)大小；
S33、根据夹角Angle(Vi)大小判断是否需要在节点Vi-1和Vi+1之间构建新的局部路径进行局部路径的光顺，若是，则在节点Vi-1和Vi+1之间构建新的局部路径；反之，则不在节点Vi-1和Vi+1之间构建新的局部路径；
所述“根据夹角Angle(Vi)大小判断是否需要在节点Vi-1和Vi+1之间构建新的局部路径进行局部路径的光顺，若是，则在节点Vi-1和Vi+1之间构建新的局部路径”这一步骤具体包括：
根据夹角Angle(Vi)大小计算节点Vi的光顺权值P(Vi)，并将节点Vi的光顺权值P(Vi)存入光顺权值集合W中；
取出权值集合W中的最大值wmax，得到其对应的节点Vi及关联点集Connect(Vi-1,Vi,Vi+1)，接下来需要将wmax对应的节点Vi移动到合适的位置以完成局部路径的光顺，Vi新位置的寻找方法如下：
若节点Vi属于网格顶点，则遍历节点Vi的一环邻域点集合Ei{X1,X2,…,Xq}中的每一点，计算向量ViXi与向量Vi-1Xi间的夹角和向量ViXi与向量Vi+1Xi间的夹角这两个夹角之和θi，然后判断θi是否小于π，若是，则得到节点Vi相交边ViXi的集合Ci，然后构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，接着，计算平面Pane与集合Ci中相交边的交点，最后，根据计算的交点构建新的局部路径：若平面Pane与集合Ci中的部分相交边有交点，则计算集合Ci中与Pane无交点的相交边非节点端点以及节点端点Vi到平面Pane的距离大小，将距离平面Pane最近的端点作为这些相交边与Pane的交点，并将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；若平面Pane与集合Ci中的全部相交边均有交点，则将节点Vi-1和Vi+1与集合Ci中所有相交边与Pane的交点依次连接起来作为新的局部路径；反之，则直接将节点Vi-1和Vi+1连接起来作为新的局部路径，其中，q为Ei的点的总数，Xi为Vi对应的的一环邻域点，集合Ci内的元素按照θi从小到大进行排序，给定的平面法矢
的值为相交边所在三角片法矢的加权平均；
若节点Vi属于网格边上的点，以节点Vi所在的边作为相交边，构建经过节点Vi-1和Vi+1并平行于给定的平面法矢
的平面Pane，然后计算平面Pane与相交边的交点，若平面Pane与相交边的交点不存在，则计算节点端点Vi以及节点Vi所在边的两个端点到平面P的距离大小，将距离平面最近的端点作为交点，所述给定的平面法矢
的值为相交边所在三角片法矢的加权平均；接着，将节点Vi从特征路径Path中删除，并在当前位置将交点顺序加入到特征路径Path中，最后，对与节点Vi相关的节点及其权值信息进行相应的更新：首先从权值集合W中删除节点Vi的关联点集Connect(Vi-1,Vi,Vi+1)中各节点的光顺权值，然后更新节点Vi的关联点集中两相邻节点Vi-1和Vi+1对应的关联点集及光顺权值，最后将更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值存入权值集合W中，其中，更新后的关联点集Connect(Vi-1,Vi,Vi+1)对应的光顺权值计算公式为：
其中，delete为删除操作，updadte为更新操作，Vi+2为Vi+1的后一节点。
2.根据权利要求1所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述步骤S1包括：
S11、在牙齿预备体的特征区域拾取特征路径的起始点Vs和末端点Ve，将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中，其中，Vs的代价信息根据综合约束控制函数计算得到；
S12、取出NodeMap的第一个元素Vtop，并判断Vtop＝Ve是否成立，若是，则结束搜索过程，转到步骤S14；反之，则执行步骤S13；
S13、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空，然后遍历Neibourhood中的每一个点，将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中，最后返回步骤S12，其中，o为Vtop的一环邻域点总数；
S14、顺次连接起始点Vs和末端点Ve之间的所有路径节点，构造出初始特征路径。
3.根据权利要求2所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述步骤S11包括：
在牙齿预备体的特征区域标记特征路径可能的起始点和末端点；
判断标记的起始点和末端点是否为网格顶点，若是，则以标记的起始点和末端点作为拾取的起始点Vs和末端点Ve，反之，则从标记的起始点或末端点所在网格边的两端点和所在三角片的三个顶点中选择一个主曲率值最高的点作为拾取的起始点Vs或末端点Ve；
根据综合约束控制函数计算起始点Vs的代价信息Fs，所述综合约束控制函数f(g)的表达式为：f(g)＝α*fdis+β*fdir+γ*(fc+ft)，Fs＝f(s)，其中，f(g)为当前候选点Vg的综合约束控制函数，fdis为当前候选点Vg与末端点Ve之间的欧式几何距离，fdir为当前候选点Vg与当前搜索中心点Vc构建的当前路径方向相对于前一步路径方向的角度变化值；fc为当前候选点Vg最大主曲率的倒数，ft为当前候选点Vg相对于当前搜索中心点Vc最小主方向Tc的变化值，α、β和γ分别为距离控制因子、方向控制因子和特征控制因子；
将Vs以及Vs的代价信息Fs存入关联容器NodeMap中，并将Vs加入候选队列NodeQueue中。
4.根据权利要求3所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述综合约束控制函数f(g)中fc和ft的计算过程为：
估算当前点p的法矢Ni，所述当前点p的法矢Ni的估算公式为：
其中，当前点p为当前候选点Vg或当前搜索中心点Vc，lj和nj分别表示当前点p在一环邻域内第j个三角片的对边和法矢，m表示当前点p在一环邻域内网格顶点的个数；
根据当前点p的法矢Ni在网格顶点处建立局部坐标系p-uvw，所述局部坐标系p-uvw由全局坐标系O-XYZ执行以下操作得到：将O点平移到p点，将Z轴旋转至与法矢Ni重合，将X和Y轴分别取为u和v；
在局部坐标系p-uvw中根据网格顶点的2环邻域点确定相应的局部三次曲面方程f(x,y)，然后通过最小二乘拟合方法求解f(x,y)的系数值，进而确定相应的Weingarten曲率矩阵W1，所述局部三次曲面方程f(x,y)的表达式为：
其中，A、B、C、D、E、F和G均为f(x,y)的系数值，
对Weingarten曲率矩阵W1进行奇异值分解，得到与当前网格顶点的最大主曲率和最小主曲率相对应的两个特征值，然后根据这两个特征值对应的两个特征向量t1和t2确定当前网格顶点的最大主方向tmax和最小主方向tmin，其中，t1、t2、tmax和tmin的表达式分别为：t1＝(t11,t12)T，t2＝(t21,t22)T，tmax＝t11x+t12y，tmin＝t21x+t22y；
根据t1、t2、tmax和tmin确定综合约束控制函数f(g)中fc和ft的值。
5.根据权利要求3所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述步骤S13包括：
S131、以Vtop作为当前搜索中心点，查找出Vtop的一环邻域点集Neibourhood{V1,V2,…,Vo}，并将关联容器NodeMap清空；
S132、以Neibourhood中的任意一点Vn作为当前候选点进行遍历分析，最后将Neibourhood中未在候选队列NodeQueue中的点及相应的代价信息存入关联容器NodeMap中；
S133、判断Neibourhood中所有点的遍历分析是否都已完成，若是，则返回步骤S12，反之，则返回步骤S132。
6.根据权利要求5所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述步骤S132包括：
根据综合约束控制函数f(g)计算当前候选点Vn的代价信息Cost；
判断当前候选点Vn是否已在候选队列NodeQueue中，若是，则表明当前候选点Vn之前已作为候选点被搜索过，此时对该点不做任何处理，继续选取Neibourhood中的下一个候选点来计算其代价信息；反之，则将Vn加入到候选队列NodeQueue中，并令Fn＝Cost，最后将Vn以及Fn存入关联容器NodeMap中。
7.根据权利要求1所述的一种牙齿预备体修复模型特征线提取方法，其特征在于：所述步骤S2具体为：
判断关键特征点间的初始特征路径是否与期望的特征线存在偏差，若是，则采用点编辑操作对关键特征点间的初始特征路径进行编辑，反之则执行步骤S3，所述点编辑操作包括但不限于删除点操作、***点操作和移动点操作。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610634704.3A CN107689049B (zh) | 2016-08-03 | 2016-08-03 | 一种牙齿预备体修复模型特征线提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610634704.3A CN107689049B (zh) | 2016-08-03 | 2016-08-03 | 一种牙齿预备体修复模型特征线提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107689049A CN107689049A (zh) | 2018-02-13 |
| CN107689049B true CN107689049B (zh) | 2020-06-16 |
Family
ID=61151652
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610634704.3A Expired - Fee Related CN107689049B (zh) | 2016-08-03 | 2016-08-03 | 一种牙齿预备体修复模型特征线提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107689049B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108491850B (zh) * | 2018-03-27 | 2020-04-10 | 北京正齐口腔医疗技术有限公司 | 三维牙齿网格模型的特征点自动提取方法及装置 |
| CN110335297B (zh) * | 2019-06-21 | 2021-10-08 | 华中科技大学 | 一种基于特征提取的点云配准方法 |
| CN111063031B (zh) * | 2019-12-06 | 2021-09-28 | 南京航空航天大学 | 一种陶瓷基复合材料铺层纱线铺敷区域的限定方法 |
| CN113592763A (zh) * | 2020-04-30 | 2021-11-02 | 深圳云甲科技有限公司 | 基于曲率方向的桩核检测方法和装置 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101477707A (zh) * | 2008-12-18 | 2009-07-08 | 浙江大学 | 一种插值给定若干闭合曲线的曲面造型方法 |
| CN104699865A (zh) * | 2013-12-09 | 2015-06-10 | 南京智周信息科技有限公司 | 一种数字化口腔固定修复的方法及装置 |
| EP3018461A1 (en) * | 2014-11-07 | 2016-05-11 | 3M Innovative Properties Company | A method of making a dental restoration |
-
2016
- 2016-08-03 CN CN201610634704.3A patent/CN107689049B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101477707A (zh) * | 2008-12-18 | 2009-07-08 | 浙江大学 | 一种插值给定若干闭合曲线的曲面造型方法 |
| CN104699865A (zh) * | 2013-12-09 | 2015-06-10 | 南京智周信息科技有限公司 | 一种数字化口腔固定修复的方法及装置 |
| EP3018461A1 (en) * | 2014-11-07 | 2016-05-11 | 3M Innovative Properties Company | A method of making a dental restoration |
Non-Patent Citations (1)
| Title |
|---|
| "口腔修复嵌体造型关键技术研究";张长东;《中国博士学位论文全文数据库医药卫生科技辑》;20150115;论文第3.2-3.4节 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107689049A (zh) | 2018-02-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107689049B (zh) | 一种牙齿预备体修复模型特征线提取方法 | |
| Pauly et al. | Example-based 3d scan completion | |
| CN108510577B (zh) | 基于已有动作数据的真实感动作迁移和生成方法及*** | |
| WO2018059155A1 (zh) | 带有几何误差的三维实体模型的构建方法及计算机可读存储介质 | |
| Alliez et al. | Isotropic surface remeshing | |
| JP2753085B2 (ja) | 形状モデリング方法及びその装置 | |
| US8004517B1 (en) | Methods, apparatus and computer program products that model three-dimensional surface structures | |
| CN108470365B (zh) | 一种基于上下牙颌的牙弓线绘制方法 | |
| JP4832991B2 (ja) | 所要の幾何的連続を有するパラメトリック曲面を生成するプロセス | |
| JP4832990B2 (ja) | パラメータ化された曲面のアイソトポロジックな集合をメッシュから生成する方法 | |
| CN108711194B (zh) | 一种基于三次b样条插值的三维网格模型拼接方法 | |
| WO2019214339A1 (zh) | 牙齿三维数字模型的局部坐标系设定方法 | |
| Yeh et al. | Template-based 3d model fitting using dual-domain relaxation | |
| Veltkamp | Closed object boundaries from scattered points | |
| WO2015188445A1 (zh) | 点云三维模型重建方法及*** | |
| US20120203513A1 (en) | Reconstruction of non-visible part of tooth | |
| CN110176073B (zh) | 三维缺陷模型自动建模和自适应分层方法 | |
| CN109993751B (zh) | 基于凹陷感知调和标量场的牙颌半自动精确分割算法 | |
| CN110689620B (zh) | 一种多层次优化的网格曲面离散样条曲线设计方法 | |
| Qiu et al. | An efficient and collision-free hole-filling algorithm for orthodontics | |
| CN108242056A (zh) | 一种基于调和场算法的三维牙齿网格数据的分割方法 | |
| CN110097642A (zh) | 一种基于半边结构的模型网格补全的方法 | |
| Wuttke et al. | Quality preserving fusion of 3d triangle meshes | |
| CN117172399B (zh) | 一种基于启发式算法的自动铺丝轨迹规划方法 | |
| CN115830287B (zh) | 基于激光口扫与cbct重建的牙齿点云融合方法、设备及介质 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20200616 |