CN111476887A - 一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 - Google Patents
一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 Download PDFInfo
- Publication number
- CN111476887A CN111476887A CN202010265986.0A CN202010265986A CN111476887A CN 111476887 A CN111476887 A CN 111476887A CN 202010265986 A CN202010265986 A CN 202010265986A CN 111476887 A CN111476887 A CN 111476887A
- Authority
- CN
- China
- Prior art keywords
- inclined plane
- curve
- tooth preparation
- point
- interpolation
- 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.)
- Granted
Links
- 238000002360 preparation method Methods 0.000 title claims abstract description 93
- 238000000034 method Methods 0.000 title claims abstract description 40
- 238000007781 pre-processing Methods 0.000 claims abstract description 12
- 239000013598 vector Substances 0.000 claims description 52
- 238000004364 calculation method Methods 0.000 claims description 25
- 230000008569 process Effects 0.000 claims description 13
- 238000012545 processing Methods 0.000 claims description 9
- 238000013507 mapping Methods 0.000 claims description 7
- 238000005457 optimization Methods 0.000 claims description 5
- 239000004576 sand Substances 0.000 claims description 4
- 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 claims description 3
- 238000013459 approach Methods 0.000 claims description 3
- 238000012217 deletion Methods 0.000 claims description 3
- 230000037430 deletion Effects 0.000 claims description 3
- 238000013178 mathematical model Methods 0.000 claims description 3
- 239000011159 matrix material Substances 0.000 claims description 3
- 230000000452 restraining effect Effects 0.000 claims description 3
- 238000000638 solvent extraction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000002474 experimental method Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 2
- 208000002925 dental caries Diseases 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 210000000746 body region Anatomy 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 208000037921 secondary disease Diseases 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61C—DENTISTRY; APPARATUS OR METHODS FOR ORAL OR DENTAL HYGIENE
- A61C7/00—Orthodontics, i.e. obtaining or maintaining the desired position of teeth, e.g. by straightening, evening, regulating, separating, or by correcting malocclusions
- A61C7/002—Orthodontic computer assisted systems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4023—Scaling of whole images or parts thereof, e.g. expanding or contracting based on decimating pixels or lines of pixels; based on inserting pixels or lines of pixels
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61C—DENTISTRY; APPARATUS OR METHODS FOR ORAL OR DENTAL HYGIENE
- A61C7/00—Orthodontics, i.e. obtaining or maintaining the desired position of teeth, e.g. by straightening, evening, regulating, separating, or by correcting malocclusions
- A61C7/002—Orthodontic computer assisted systems
- A61C2007/004—Automatic construction of a set of axes for a tooth or a plurality of teeth
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/10—Internal combustion engine [ICE] based vehicles
- Y02T10/40—Engine management systems
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Mathematical Physics (AREA)
- Geometry (AREA)
- Computational Mathematics (AREA)
- Software Systems (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computing Systems (AREA)
- Oral & Maxillofacial Surgery (AREA)
- Dentistry (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Epidemiology (AREA)
- Animal Behavior & Ethology (AREA)
- Algebra (AREA)
- Veterinary Medicine (AREA)
- Databases & Information Systems (AREA)
- Public Health (AREA)
- Computer Graphics (AREA)
- Numerical Control (AREA)
Abstract
一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,它涉及机器人辅助牙体预备技术领域,本发明将功能尖斜面备牙阶段规划为机器人的三维曲面轨迹规划,基于医生的标准化后牙手工预备体模型,预处理后得到了功能尖斜面的曲面数据,根据牙齿之间的相邻关系和功能将功能尖斜面划分为三个曲面并分别得到了各个曲面的插值结果,再根据曲面的型值点和曲面片特征完成了三个曲面的备牙轨迹规划,建立了一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法。本发明使得处理后的预备体模型曲率更加光滑,优化了备牙轨迹的运算速率,有效提高了备牙曲线的预备精度。
Description
技术领域
本发明专利涉及一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,属于机器人辅助牙体预备技术领域。
背景技术
龋齿是导致牙体缺损的重要原因,严重的影响了人们的口腔健康。口腔修复是治疗牙体缺损的重要手段,牙体预备是修复过程中的必要治疗环节,是指医生对患龋齿处的牙齿进行定量去除并形成预期三维形状的操作过程。在传统的牙体预备过程中,需要依赖于医生手工的操作并结合丰富的临床经验对患龋牙齿进行大量的重复的微细调整。然而目前我国医患比失衡的现状无法满足当前极大的牙体预备需求量,因此需要引入机器人和辅助软件来代替或辅助医生高效的完成牙体预备工作,不仅能够缓解医患失衡的难题,而且还能有效的提高备牙质量和口腔治疗的效果。
在机器辅助完成牙体预备的过程中,机器人末端的运动轨迹决定了牙体预备的精度,进而影响治疗的效果和成功率。因此,获取合理的功能尖斜面备牙轨迹,着重考虑预备体与后牙修复体的贴合程度,从而减少因牙体预备不当所形成口腔继发疾病的情况,是目前机器人辅助牙体预备领域的研究热点。
发明内容
针对上述问题,本发明提出一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,解决目前机器人辅助牙体预备技术领域缺少获得功能尖斜面备牙轨迹的方法,以提高牙体预备的精度和成功率,进而实现机器人辅助牙体预备。
本发明为解决上述问题所采取的方案为:
一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,所述方法的具体实现过程为:
步骤一、预备体模型预处理
根据医生临床的备牙模型进行扫描,得到obj格式的标准化预备体模型,选用Geomagic Wrap逆向工程软件对obj格式的标准化备牙三维模型进行处理,对生成三角形曲面片和四边形曲面片进行曲面片形状参数处理和曲面特征处理,以获取更加光滑的预备体模型和优化计算速度;
a)建立牙体预备功能尖斜面曲面映射关系:
计算参考二维平面Θ的单位法向量n,其方向垂直于参考二维平面Θ的z轴方向:
式中:Ns为区域Ω内包含的四边形曲面片数量;ni为曲面片i的法向量;
根据法向量u确定任意一个不与u共线的v单位向量,构建投影的参考二维平面Θ,其中,参考二维平面Θ的y方向为Θy=u×v,参考二维平面Θ的x方向为Θx=Θy×u,构造三维曲面区域Ω与参考二维平面Θ之间的映射矩阵H(u,v):
式中:u∈[0,1];v∈[0,1];h1(u)、h2(v)、h3(u)、h4(v)为四边形曲面片的四个边界;a1、a2、a3、a4为四边形的四个顶点;
b)对牙体预备功能尖斜面四边形曲面片形状进行约束:
根据公式(2)、(3)获取理想化曲面片和非理想化曲面片,对于非理想化曲面片需要进行边界夹角θi和边界间距L约束,便于后续整体曲面的生成,其中θi∈[60,120°],由a4做延长线与h2(v)得到b4,b4的 x轴坐标为x4+Lmin,y轴坐标与x4相同,新构成的a4b4为夹角θ1所对应的最小间距推导出的计算公式:
c)对牙体预备功能尖斜面四边形曲面片进行优化处理:
通过移动曲面片的节点位置完成对曲面片的优化,采取优化曲面片网格的拓扑结构的方法,当内部节点a的相邻边单元数Nk(a)=2,则删除该节点;边删除需要检查所有2个内部节点a和b连接成边ef,此时分为两种情况:若Nk(a)=Nk(b)=2且ef的2个相邻单元中不存在边界,则删除ef;若 [Nk(a)+Nk(d)]≥[Nk(d)+Nk(b)],则用bc做新边,反之用ad作为新边;对角线变换是以a和b相连为对角线为例,若[Nk(a)+Nk(b)]≥9,需要考虑三种情况下的最优对角线:当[Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]>[Nk(a)+Nk(b)]- [Nk(e)+Nk(f)]≥3时,用cd代替ab;当[Nk(a)+Nk(b)]-[Nk(e)+Nk(f)]>[Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]≥3时,用ef代替ab;
步骤二、功能尖斜面的区域划分
依据预处理阶段中产生的四边形曲面片S数量Ns、三角形曲面片T数量Nt和特殊交点M数量Nm作为划分参数,控制四边形曲面片S数量Ns和三角形曲面片T数量Nt的比例Rs在0.05以内,其中特殊交点 M是指在拓扑优化阶段相邻边单元数Nk>4的节点;将功能尖斜面以唇侧为参考方向被分为左邻接面、中间尖斜面和右邻接面,其中,左邻接面为Ns=35,Nt=3,Nm=2;中间尖斜面为Ns=64,Nt=0,Nm=0;右邻接面为Ns=31,Nt=3,Nm=2;
步骤三、功能尖斜面的节点矢量U和V的确定和控制点Pi,j的反算
采用积累弦长参数化法确定u方向的第i行型值点Qi,j,j=0、1、...、m,i=0、1、2、...、n,并设u 方向的第i行的参数为ui,j:
则u方向和v方向参数化可以取所有行的相同列参数的算术平均值得到公式(7)、(8):
由于基于NURBS曲面的功能尖斜面有两个方向,因此需要先计算p次,u方向的第i行的型值点Qi,j进行曲线的反算,进而可以构造出第i行的截面线,得到中间控制点Ri,j。然后固定v方向的j坐标,再计算q次,v方向的第j列进行反算,其中Ri,j作为此次反算的数据点,反算结果可以求出第j列曲线的控制点Pi,j,则得到的Pi,j即为曲面的控制点;
步骤四、功能尖斜面的曲面插值
功能尖斜面的插值过程通过MATLAB软件编程实现,通过定义一个插值函数“Functional_Bevel_interpolation”,其中的变量为u方向与v方向的次数p和q、两个方向的Nu和Nv以及曲面型值点Qi,j的xyz坐标,输出得到节点矢量U和V以及控制点Px,Py和Pz,将该函数输出的值输入至“plot_ Functional_Bevel”绘制图形函数,再输入网格生成总量SUM的值就可以完成功能尖斜面的曲面插值;
左邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=13,Nv=20,型值点NQ=260,SUM=100,输出节点矢量U数量为16个,节点矢量V数量为23个,控制点数为260个;
中间尖斜面的输入的参数为u方向和v方向的次数为p=q=3,Nu=9,Nv=25,型值点NQ=225,SUM=100,输出节点矢量U数量为12个,节点矢量V数量为28个,控制点数为225个;
右邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=14,Nv=20,型值点NQ=280,SUM=100,输出节点矢量U数量为17个,节点矢量V数量为23个,控制点数为280个;
步骤五、生成功能尖斜面的备牙轨迹
根据u方向和v方向将曲面划分为两种参数线族,再对两个方向上的曲线长度进行比较,选取其中较长的参考线作为机器人末端轨迹参考线,而另一个参考线方向作为进给方向;根据步骤三得到的型值点 Qi,j,为了使机器人精确的完成预期的邻面备牙曲线,需要将曲线进行划分,理论上划分的段数越多机器人的末端轨迹就越逼近预期曲线,引入NURBS曲线数学模型;
ωi(i=0,1,...,n)为权重因子;Pi(i=0,1,...,n)为控制顶点,其数量为n+1;p为NURBS曲线次数;Ni,p(u)为基函数;公式(9)中的Ni,p(u)为:
式(10)中:ui是非均匀节点向量U的元素,如公式(11)所示;
m+1为节点向量U的长度;m、p、n三者的关系为m=n+p+1;a和b通常为0和1;
插补原理是利用时间序列{t1,t2,...,tk,...,tn-1,tn}分割参数序列{u1,u2,...,uk,uk+1,...,un-1,un},进而得到插补点序列{C(u1),C(u2),...,C(uk),...,C(un-1),C(un)},插补的核心计算就是利用插补周期T求得uk和 uk+1之间的关系,进一步通过C(uk)得到C(uk+1);
首先,对NURBS曲线的导数进行推导,得到曲线对u的一阶导数C(1)(u);
式中:C(1)(u)为曲线对u的一阶导数;
继续对曲线求二阶导数得到C(2)(u):
式中:C(2)(u)为曲线对u的二阶导数,基函数Ni,p (m)(u)的表达式为:
将NURBS曲线上的进给速率表示为:
整理变换公式(15)得到公式(16):
式中:Cx (1)(uk)为曲线x方向的一阶导;Cy (1)(uk)为曲线y方向的一阶导;Cz (1)(uk)为曲线z方向的一阶导。
继续对公式(16)求二阶导数,得到:
式中:Cx (2)(uk)为曲线x方向的二阶导;Cy (2)(uk)为曲线y方向的二阶导;Cz (2)(uk)为曲线z方向的二阶导;
整理得到了基于NURBS曲线的功能尖斜面备牙轨迹的几何特性和运动特性之间的关系:
本发明的有益效果为:
1.本发明提出了一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,在预备体模型的预处理阶段,对三角形曲面片进行筛选和优化,控制四边形曲面片S数量Ns和三角形曲面片T数量Nt的比例Rs在0.05以内,并对非理想化曲面片边界和夹角设定约束条件,优化曲面片网格的拓扑结构,使得处理后的预备体模型曲率更加光滑,便于后续备牙曲线的反算和精度保障。
2.本发明在功能尖斜面的划分的区域划分阶段,以牙齿相邻特点为参考,依据预处理阶段中产生的四边形曲面片S数量Ns、三角形曲面片T数量Nt和特殊交点M数量Nm作为划分参数,发现三角形曲面片T 和特殊交点M均分布在发生弧形变化的两侧,将功能尖斜面以唇侧为参考方向被分为左邻接面、中间尖斜面和右邻接面,得到u方向和v方向的两条具有延伸特性的边界曲线,解决了功能尖斜面的个性化需求,优化了备牙轨迹的运算速率。
3.本发明在功能尖斜面备牙轨迹的反算过程中,考虑到预备过程中对目标牙齿邻牙的无效接触,利用累弦长参数化法所得曲线的特性,反应了离散点Q按照弦长的分布情况,有效提高了备牙轨迹的预备精度,功能尖斜面备牙曲线保持一致且光滑连续。
附图说明
为了易于说明,本发明由下述的具体实施及附图作以详细描述。
图1为功能尖斜面数据的预处理流程图;
图2为功能尖斜面三维曲面区域Ω与参考二维平面Θ的映射过程;
图3为功能尖斜面非理想化曲面片边界间最短间距;
图4为功能尖斜面曲面控制点的反算流程;
图5为功能尖斜面曲面插值的流程图;
图6为标准后牙预备体模型预处理流程
图7为预处理后标准化预备体模型;
图8为功能尖斜面曲面划分结果;
图9为左邻接面插值结果;
图10为中间尖斜面插值结果;
图11为右邻接面插值结果;
图12为左邻接面备牙轨迹;
图13为中间尖斜面备牙轨迹;
图14为右邻接面备牙轨迹;
图15为机器人辅助牙体预备实验系统;
图16为机器人完成功能尖斜面预备阶段的实验过程;
图17为功能尖斜面预备阶段特征点的选取。
具体实施方式
为使本发明专利的目的、技术方案和优点更加清楚明了,下面通过附图中示出的具体实施例来描述本发明专利,但是应该理解,这些描述只是示例性的,而并非要限制本发明专利的范围,此外,在以下说明中,省略了对公知结构和技术的描述,以避免不必要地混淆本发明专利的概念。
实施实例一:
如图1、图2、图3、图4、图5所示,本具体实施方式采用以下技术方案:一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,其特征在于:所述方法的具体实现过程为:
步骤一、预备体模型预处理
根据医生临床的备牙模型进行扫描,得到obj格式的标准化预备体模型,选用Geomagic Wrap逆向工程软件对obj格式的标准化备牙三维模型进行处理,对生成三角形曲面片和四边形曲面片进行曲面片形状参数处理和曲面特征处理,以获取更加光滑的预备体模型和优化计算速度;
a)建立曲面映射关系:
计算参考二维平面Θ的单位法向量n,其方向垂直于参考二维平面Θ的z轴方向:
式中:Ns为区域Ω内包含的四边形曲面片数量;ni为曲面片i的法向量;
根据法向量u确定任意一个不与u共线的v单位向量,构建投影的参考二维平面Θ,其中,参考二维平面Θ的y方向为Θy=u×v,参考二维平面Θ的x方向为Θx=Θy×u,构造三维曲面区域Ω与参考二维平面Θ之间的映射矩阵H(u,v):
式中:u∈[0,1];v∈[0,1];h1(u)、h2(v)、h3(u)、h4(v)为四边形曲面片的四个边界;a1、a2、a3、a4为四边形的四个顶点;
b)对四边形曲面片形状进行约束:
根据公式(2)、(3)获取理想化曲面片和非理想化曲面片,对于非理想化曲面片需要进行边界夹角θi和边界间距L约束,便于后续整体曲面的生成,其中θi∈[60,120°],由a4做延长线与h2(v)得到b4,b4的 x轴坐标为x4+Lmin,y轴坐标与x4相同,新构成的a4b4为夹角θ1所对应的最小间距推导出的计算公式:
c)对四边形曲面片进行优化处理:
通过移动曲面片的节点位置完成对曲面片的优化,采取优化曲面片网格的拓扑结构的方法,当内部节点a的相邻边单元数Nk(a)=2,则删除该节点;边删除需要检查所有2个内部节点a和b连接成边ef,此时分为两种情况:若Nk(a)=Nk(b)=2且ef的2个相邻单元中不存在边界,则删除ef;若 [Nk(a)+Nk(d)]≥[Nk(d)+Nk(b)],则用bc做新边,反之用ad作为新边;对角线变换是以a和b相连为对角线为例,若[Nk(a)+Nk(b)]≥9,需要考虑三种情况下的最优对角线:当 [Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]>[Nk(a)+Nk(b)]-[Nk(e)+Nk(f)]≥3时,用cd代替ab;当 [Nk(a)+Nk(b)]-[Nk(e)+Nk(f)]>[Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]≥3时,用ef代替ab;
步骤二、功能尖斜面的区域划分
依据预处理阶段中产生的四边形曲面片S数量Ns、三角形曲面片T数量Nt和特殊交点M数量Nm作为划分参数,控制四边形曲面片S数量Ns和三角形曲面片T数量Nt的比例Rs在0.05以内,其中特殊交点 M是指在拓扑优化阶段相邻边单元数Nk>4的节点;将功能尖斜面以唇侧为参考方向被分为左邻接面、中间尖斜面和右邻接面;
步骤三、功能尖斜面的节点矢量U和V的确定和控制点Pi,j的反算
采用积累弦长参数化法确定u方向的第i行型值点Qi,j,j=0、1、...、m,i=0、1、2、...、n,并设u 方向的第i行的参数为ui,j:
则u方向和v方向参数化可以取所有行的相同列参数的算术平均值得到公式(7)、(8):
由于基于NURBS曲面的功能尖斜面有两个方向,因此需要先计算p次,u方向的第i行的型值点Qi,j进行曲线的反算,进而可以构造出第i行的截面线,得到中间控制点Ri,j。然后固定v方向的j坐标,再计算q次,v方向的第j列进行反算,其中Ri,j作为此次反算的数据点,反算结果可以求出第j列曲线的控制点Pi,j。则得到的Pi,j即为曲面的控制点;
步骤四、功能尖斜面的曲面插值
功能尖斜面的插值过程通过MATLAB软件编程实现,通过定义一个插值函数“Functional_Bevel_interpolation”,其中的变量为u方向与v方向的次数p和q、两个方向的Nu和Nv以及曲面型值点Qi,j的xyz坐标,输出得到节点矢量U和V以及控制点Px,Py和Pz,将该函数输出的值输入至“plot_ Functional_Bevel”绘制图形函数,再输入网格生成总量SUM的值就可以完成功能尖斜面的曲面插值;
步骤五、生成功能尖斜面的备牙轨迹
根据u方向和v方向将曲面划分为两种参数线族,再对两个方向上的曲线长度进行比较,选取其中较长的参考线作为机器人末端轨迹参考线,而另一个参考线方向作为进给方向;根据步骤三得到的型值点 Qi,j,为了使机器人精确的完成预期的邻面备牙曲线,需要将曲线进行划分,理论上划分的段数越多机器人的末端轨迹就越逼近预期曲线,引入NURBS曲线数学模型;
ωi(i=0,1,...,n)为权重因子;Pi(i=0,1,...,n)为控制顶点,其数量为n+1;p为NURBS曲线次数;Ni,p(u)为基函数;公式(9)中的Ni,p(u)为:
式(10)中:ui是非均匀节点向量U的元素,如公式(11)所示;
m+1为节点向量U的长度;m、p、n三者的关系为m=n+p+1;a和b通常为0和1;
插补原理是利用时间序列{t1,t2,...,tk,...,tn-1,tn}分割参数序列{u1,u2,...,uk,uk+1,...,un-1,un},进而得到插补点序列{C(u1),C(u2),...,C(uk),...,C(un-1),C(un)},插补的核心计算就是利用插补周期T求得uk和 uk+1之间的关系,进一步通过C(uk)得到C(uk+1);
首先,对NURBS曲线的导数进行推导,得到曲线对u的一阶导数C(1)(u);
式中:C(1)(u)为曲线对u的一阶导数;
继续对曲线求二阶导数得到C(2)(u):
式中:C(2)(u)为曲线对u的二阶导数,基函数Ni,p (m)(u)的表达式为:
将NURBS曲线上的进给速率表示为:
整理变换公式(15)得到公式(16):
式中:Cx (1)(uk)为曲线x方向的一阶导;Cy (1)(uk)为曲线y方向的一阶导;Cz (1)(uk)为曲线z方向的一阶导。
继续对公式(16)求二阶导数,得到:
式中:Cx (2)(uk)为曲线x方向的二阶导;Cy (2)(uk)为曲线y方向的二阶导;Cz (2)(uk)为曲线z方向的二阶导;
整理得到了基于NURBS曲线的功能尖斜面备牙轨迹的几何特性和运动特性之间的关系:
实施实例二
如图6、图7、图8、图9、图10、图11、图12、图13、图14、图15、图16、图17所示,下面以具体实验数据为例,来对本发明方法进行举例。
依据预处理阶段中产生的四边形曲面片S数量Ns、三角形曲面片T数量Nt和特殊交点M数量Nm作为划分参数,发现三角形曲面片T和特殊交点M均分布在发生弧形变化的两侧,将功能尖斜面以唇侧为参考方向被分为左邻接面、中间尖斜面和右邻接面,其中三个曲面的参数分别为:左邻接面为Ns=35,Nt=3, Nm=2;中间尖斜面为Ns=64,Nt=0,Nm=0;右邻接面为Ns=31,Nt=3,Nm=2;
插值过程通过MATLAB软件编程实现,通过定义一个插值函数“Functional_Bevel_interpolation”,其中的变量为u方向与v方向的次数p和q、两个方向的Nu和Nv以及曲面型值点Qi,j的xyz坐标,输出得到节点矢量U和V以及控制点Px,Py和Pz,将该函数输出的值输入至“plot_Functional_Bevel”绘制图形函数,再输入网格生成总量SUM的值就可以完成曲面的插值;
左邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=13,Nv=20,型值点NQ=260,SUM=100,输出节点矢量U数量为16个,节点矢量V数量为23个,控制点数为260个;
中间尖斜面的输入的参数为u方向和v方向的次数为p=q=3,Nu=9,Nv=25,型值点NQ=225,SUM=100,输出节点矢量U数量为12个,节点矢量V数量为28个,控制点数为225个;
右邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=14,Nv=20,型值点NQ=280,SUM=100,输出节点矢量U数量为17个,节点矢量V数量为23个,控制点数为280个;
生成的功能尖斜面备牙轨迹是由曲面的两个方向的参考线作为机器人末端的接触点路径所确定,根据 u方向和v方向将曲面划分为两种参数线族,再对两个方向上的曲线长度进行比较,选取其中较长的参考线作为机器人末端轨迹参考线,而另一个参考线方向作为进给方向;
左邻接面的输入的参数轨迹次数为k=3,型值点NQ=260,插补周期T=0.002s,输出节点矢量U数量为266个,控制点数为259个,插补点数为385个;
中间尖斜面的输入的参数轨迹次数为k=3,型值点NQ=225,插补周期T=0.002s,输出节点矢量U数量为231个,控制点数为225个,插补点数为438个;
右邻接面的输入的参数轨迹次数为k=3,型值点NQ=280,插补周期T=0.002s,输出节点矢量U数量为286个,控制点数为280个,插补点数为367个;
将所规划的机器人功能尖斜面预备轨迹的插补点通过空间坐标转换得到机器人在关节坐标系下的末端坐标点并输入到上位机的软件中用于机器人完成牙体预备。实验过程中末端车针直径为1.6mm,实验电压为12V,电流0.25A,末端转速约为46000r/min,机器人移动速度为11.1mm/s为设置存点回放速度的15%。
在咬合面区域选取横向沟壑作为特征曲线,选取特征曲线上的两端节点①和④以及曲线最低端的两个特征点②和③。在牙尖区域选取牙尖区域的一周曲线中的四个尖端作为⑤、⑥、⑦和⑧,上述四点同时也是三个曲面划分部分中的衔接区域。最后,在颊面区域选取两个特征点作为⑨和⑩,将功能尖斜面预备特征点理论数值的按照各区域汇总至表1,为方便后续计算将①点的x轴数值置零处理。在每个点测量5次并进行记录,共计得到50个实验数据点,如表2所示。
表1 功能尖斜面预备特征点理论数值
表2 各个特征点的实验测量数值
计算得到测量均值μ、相对定点误差ε、相对标准偏差RSD和置信区间,系统误差各项参数统计表汇总至表3。各个特征点在X方向上的相对定点误差ε范围为0.06~0.23mm,在Y方向上该误差范围为0.03~0.15mm,在Z方向上该误差范围为0.02~0.97mm,特征点误差控制在1mm以内;在相对标准误差 RSD方面,由于功能尖斜面阶段的特征点是按照图16中预备体区域划分,而三个区域的特点导致了各个特征点在Z方向上的相对标准误差较小(小于3.5%),因此重点讨论其余两个方向的相对标准误差。在咬合面区域中的①②③④特征点中X方向的相对标准误差范围为4.90~25.34%,在Y方向上该误差为 10.38~29.40%;在牙尖区域中的⑤⑥⑦⑧特征点中X方向的相对标准误差范围为4.21~25.96%,在Y方向上该误差为5.58~31.32%;在颊面区域中的⑨⑩特征点中X方向的相对标准误差为17.12%和6.12%,在Y 方向上该误差为3.19%和0.44%。通过数据可以发现相较于邻面预备阶段,该阶段的误差较高,主要原因是机器人在完成曲面规划过程中的精度不足,另外特征点的测量难度相对较大导致多次测量数值之间变化较大,使得相对标准误差升高,但误差维持在32%以内以及后续的置信区间使得测量结果在可接受范围内。最后,各个特征点在X方向上的置信区间宽度范围为0.07~0.48mm,Y方向上的置信区间宽度范围为 0.20~0.46mm,Z方向上的置信区间宽度范围为0.25~0.95mm。在同一个特征点下的不同方向的置信区间宽度平均稳定在0.30mm左右,说明了功能尖斜面阶段各个特征点的系统误差可以稳定在较小范围,保证了机器人在该阶段的预备精度。
表3 系统误差各项参数统计表
以上显示和描述了本发明专利的基本原理和主要特征和本发明专利的优点,本行业的技术人员应该了解,本发明专利不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明专利的原理,在不脱离本发明专利精神和范围的前提下,本发明专利还会有各种变化和改进,这些变化和改进都落入要求保护的本发明专利范围内。本发明专利要求保护范围由所附的权利要求书及其等效物界定。
Claims (1)
1.一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法,其特征在于:所述方法的具体实现过程为:
步骤一、预备体模型预处理
根据医生临床的备牙模型进行扫描,得到obj格式的标准化预备体模型,选用GeomagicWrap逆向工程软件对obj格式的标准化备牙三维模型进行处理,对生成三角形曲面片和四边形曲面片进行曲面片形状参数处理和曲面特征处理,以获取光滑的预备体模型和优化计算速度;
a)建立牙体预备功能尖斜面曲面映射关系:
计算参考二维平面Θ的单位法向量n,其方向垂直于参考二维平面Θ的z轴方向:
式中:Ns为区域Ω内包含的四边形曲面片数量;ni为曲面片i的法向量;
根据法向量u确定任意一个不与u共线的v单位向量,构建投影的参考二维平面Θ,其中,参考二维平面Θ的y方向为Θy=u×v,参考二维平面Θ的x方向为Θx=Θy×u,构造三维曲面区域Ω与参考二维平面Θ之间的映射矩阵H(u,v):
式中:u∈[0,1];v∈[0,1];h1(u)、h2(v)、h3(u)、h4(v)为四边形曲面片的四个边界;a1、a2、a3、a4为四边形四个顶点;
b)对牙体预备功能尖斜面四边形曲面片形状进行约束:
根据公式(2)、(3)获取理想化曲面片和非理想化曲面片,对于非理想化曲面片需要进行边界夹角θi和边界间距L约束,便于后续整体曲面的生成,其中θi∈[60°,120°],由a4做延长线与h2(v)得到b4,b4的x轴坐标为x4+Lmin,y轴坐标与x4相同,新构成的a4b4为夹角θ1所对应的最小间距推导出的计算公式:
c)对牙体预备功能尖斜面四边形曲面片进行优化处理:
通过移动曲面片的节点位置完成对曲面片的优化,采取优化曲面片网格的拓扑结构的方法,当内部节点a的相邻边单元数Nk(a)=2,则删除该节点;边删除需要检查所有2个内部节点a和b连接成边ef,此时分为两种情况:若Nk(a)=Nk(b)=2且ef的2个相邻单元中不存在边界,则删除ef;若[Nk(a)+Nk(d)]≥[Nk(d)+Nk(b)],则用bc做新边,反之用ad作为新边;对角线变换是以a和b相连为对角线为例,若[Nk(a)+Nk(b)]≥9,需要考虑三种情况下的最优对角线:当[Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]>[Nk(a)+Nk(b)]-[Nk(e)+Nk(f)]≥3时,用cd代替ab;当[Nk(a)+Nk(b)]-[Nk(e)+Nk(f)]>[Nk(a)+Nk(b)]-[Nk(c)+Nk(d)]≥3时,用ef代替ab;
步骤二、功能尖斜面的区域划分
依据预处理阶段中产生的四边形曲面片S数量Ns、三角形曲面片T数量Nt和特殊交点M数量Nm作为划分参数,控制四边形曲面片S数量Ns和三角形曲面片T数量Nt的比例Rs在0.05以内,其中特殊交点M是指在拓扑优化阶段相邻边单元数Nk>4的节点;将功能尖斜面以唇侧为参考方向被分为左邻接面、中间尖斜面和右邻接面,其中,左邻接面为Ns=35,Nt=3,Nm=2;中间尖斜面为Ns=64,Nt=0,Nm=0;右邻接面为Ns=31,Nt=3,Nm=2;
步骤三、功能尖斜面的节点矢量U和V的确定和控制点Pi,j的反算
采用积累弦长参数化法确定u方向的第i行型值点Qi,j,j=0、1、...、m,i=0、1、2、...、n,并设u方向的第i行的参数为ui,j:
则u方向和v方向参数化可以取所有行的相同列参数的算术平均值得到公式(7)、(8):
由于基于NURBS曲面的功能尖斜面有两个方向,因此需要先计算p次,u方向的第i行的型值点Qi,j进行曲线的反算,进而可以构造出第i行的截面线,得到中间控制点Ri,j;然后固定v方向的j坐标,再计算q次,v方向的第j列进行反算,其中Ri,j作为此次反算的数据点,反算结果可以求出第j列曲线的控制点Pi,j;则得到的Pi,j即为曲面的控制点;
步骤四、功能尖斜面的曲面插值
功能尖斜面的插值过程通过MATLAB软件编程实现,通过定义一个插值函数“Functional_Bevel_interpolation”,其中的变量为u方向与v方向的次数p和q、两个方向的Nu和Nv以及曲面型值点Qi,j的xyz坐标,输出得到节点矢量U和V以及控制点Px,Py和Pz,将该函数输出的值输入至“plot_Functional_Bevel”绘制图形函数,再输入网格生成总量SUM的值就可以完成功能尖斜面的曲面插值;
左邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=13,Nv=20,型值点NQ=260,SUM=100,输出节点矢量U数量为16个,节点矢量V数量为23个,控制点数为260个;
中间尖斜面的输入的参数为u方向和v方向的次数为p=q=3,Nu=9,Nv=25,型值点NQ=225,SUM=100,输出节点矢量U数量为12个,节点矢量V数量为28个,控制点数为225个;
右邻接面的输入的参数为u方向和v方向的次数为p=q=3,Nu=14,Nv=20,型值点NQ=280,SUM=100,输出节点矢量U数量为17个,节点矢量V数量为23个,控制点数为280个;
步骤五、生成功能尖斜面的备牙轨迹
根据u方向和v方向将曲面划分为两种参数线族,再对两个方向上的曲线长度进行比较,选取其中较长的参考线作为机器人末端轨迹参考线,而另一个参考线方向作为进给方向;根据步骤三得到的型值点Qi,j,为了使机器人精确的完成预期的邻面备牙曲线,需要将曲线进行划分,理论上划分的段数越多机器人的末端轨迹就越逼近预期曲线;引入NURBS曲线数学模型,
ωi(i=0,1,...,n)为权重因子;Pi(i=0,1,...,n)为控制顶点,其数量为n+1;p为NURBS曲线次数;Ni,p(u)为基函数;公式(9)中的Ni,p(u)为:
式(10)中:ui是非均匀节点向量U的元素,如公式(11)所示;
m+1为节点向量U的长度;m、p、n三者的关系为m=n+p+1;a和b通常为0和1;
插补原理是利用时间序列{t1,t2,...,tk,...,tn-1,tn}分割参数序列{u1,u2,...,uk,uk+1,...,un-1,un},进而得到插补点序列{C(u1),C(u2),...,C(uk),...,C(un-1),C(un)},插补的核心计算就是利用插补周期T求得uk和uk+1之间的关系,进一步通过C(uk)得到C(uk+1);
首先,对NURBS曲线的导数进行推导,得到曲线对u的一阶导数C(1)(u);
式中:C(1)(u)为曲线对u的一阶导数;
继续对曲线求二阶导数得到C(2)(u):
式中:C(2)(u)为曲线对u的二阶导数,基函数Ni,p (m)(u)的表达式为:
将NURBS曲线上的进给速率表示为:
整理变换公式(15)得到公式(16):
式中:Cx (1)(uk)为曲线x方向的一阶导;Cy (1)(uk)为曲线y方向的一阶导;Cz (1)(uk)为曲线z方向的一阶导;
继续对公式(16)求二阶导数,得到:
式中:Cx (2)(uk)为曲线x方向的二阶导;Cy (2)(uk)为曲线y方向的二阶导;Cz (2)(uk)为曲线z方向的二阶导;
整理得到了基于NURBS曲线的功能尖斜面备牙轨迹的几何特性和运动特性之间的关系:
左邻接面的输入的参数轨迹次数为k=3,型值点NQ=260,插补周期T=0.002s,输出节点矢量U数量为266个,控制点数为259个,插补点数为385个;
中间尖斜面的输入的参数轨迹次数为k=3,型值点NQ=225,插补周期T=0.002s,输出节点矢量U数量为231个,控制点数为225个,插补点数为438个;
右邻接面的输入的参数轨迹次数为k=3,型值点NQ=280,插补周期T=0.002s,输出节点矢量U数量为286个,控制点数为280个,插补点数为367个。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010265986.0A CN111476887B (zh) | 2020-04-04 | 2020-04-04 | 一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202010265986.0A CN111476887B (zh) | 2020-04-04 | 2020-04-04 | 一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111476887A true CN111476887A (zh) | 2020-07-31 |
CN111476887B CN111476887B (zh) | 2024-06-14 |
Family
ID=71750668
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202010265986.0A Active CN111476887B (zh) | 2020-04-04 | 2020-04-04 | 一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111476887B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117217992A (zh) * | 2023-08-22 | 2023-12-12 | 哈尔滨理工大学 | 一种基于曲率波动率及步长误差判断的牙体预备轨迹插补精度评价方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060098008A1 (en) * | 2002-06-04 | 2006-05-11 | Christof Holberg | Method, device and computer program product for generating a three-dimensional model |
CN104504759A (zh) * | 2014-12-29 | 2015-04-08 | 佛山市诺威科技有限公司 | 一种基于义齿基底冠三角网格的快速过渡缝补的方法 |
CN104699865A (zh) * | 2013-12-09 | 2015-06-10 | 南京智周信息科技有限公司 | 一种数字化口腔固定修复的方法及装置 |
CN108549325A (zh) * | 2018-05-23 | 2018-09-18 | 合肥工业大学 | 一种自由曲面弧长参数曲线加工轨迹生成方法 |
CN109035410A (zh) * | 2018-07-19 | 2018-12-18 | 浙江大学 | 一种基于离散化的多重曲面建筑网格划分方法 |
-
2020
- 2020-04-04 CN CN202010265986.0A patent/CN111476887B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20060098008A1 (en) * | 2002-06-04 | 2006-05-11 | Christof Holberg | Method, device and computer program product for generating a three-dimensional model |
CN104699865A (zh) * | 2013-12-09 | 2015-06-10 | 南京智周信息科技有限公司 | 一种数字化口腔固定修复的方法及装置 |
CN104504759A (zh) * | 2014-12-29 | 2015-04-08 | 佛山市诺威科技有限公司 | 一种基于义齿基底冠三角网格的快速过渡缝补的方法 |
CN108549325A (zh) * | 2018-05-23 | 2018-09-18 | 合肥工业大学 | 一种自由曲面弧长参数曲线加工轨迹生成方法 |
CN109035410A (zh) * | 2018-07-19 | 2018-12-18 | 浙江大学 | 一种基于离散化的多重曲面建筑网格划分方法 |
Non-Patent Citations (2)
Title |
---|
JINGANG JIANG 等: "Registration Technology of Augmented Reality in Oral Medicine: A Review", IEEE ACCESS, vol. 7, 23 April 2019 (2019-04-23) * |
周立英;韩栋伟;: "纤维桩核修复下颌第二前磨牙三维有限元模型", 同济大学学报(医学版), no. 04, 15 August 2008 (2008-08-15) * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117217992A (zh) * | 2023-08-22 | 2023-12-12 | 哈尔滨理工大学 | 一种基于曲率波动率及步长误差判断的牙体预备轨迹插补精度评价方法 |
CN117217992B (zh) * | 2023-08-22 | 2024-06-07 | 哈尔滨理工大学 | 一种基于曲率波动率及步长误差判断的牙体预备轨迹插补精度评价方法 |
Also Published As
Publication number | Publication date |
---|---|
CN111476887B (zh) | 2024-06-14 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113385486B (zh) | 一种基于线结构光的激光清洗路径自动生成系统及方法 | |
CN108227630B (zh) | 一种采用时间参数多项式插补的自由曲面数控加工方法 | |
CN102785166B (zh) | 一种基于运动学变换的数控砂轮磨削加工方法 | |
CN113741426B (zh) | 一种基于局部点云曲线拟合的机器人加工路径规划方法 | |
Ji et al. | An adaptive real-time NURBS curve interpolation for 4-axis polishing machine tool | |
CN112022382B (zh) | 一种牙套的自动切割方法及装置 | |
CN110188423A (zh) | 一种基于有限元网格划分的线性工程结构快速bim建模方法 | |
CN113276130B (zh) | 一种基于点云切片的自由曲面喷涂路径规划方法及系统 | |
CN108682043A (zh) | 一种基于参数映射的复杂曲面测量规划方法 | |
CN112033338B (zh) | 一种叶片类曲面接触式扫描测量测头半径面补偿方法 | |
CN114237161B (zh) | 一种基于数字滤波的工业机器人nurbs曲线插补方法 | |
CN111476887A (zh) | 一种用于机器人辅助牙体预备功能尖斜面备牙轨迹生成方法 | |
CN110412941A (zh) | 螺旋曲面数控包络铣削方法及其集成控制系统 | |
CN111489437B (zh) | 一种用于机器人辅助牙体预备的邻面备牙曲线生成方法 | |
Zhou et al. | A novel method to generate the tooth surface model of face-milled generated spiral bevel gears | |
CN112387995B (zh) | 一种自由曲面超精密车削后表面形貌预测方法 | |
CN111487929A (zh) | 一种基于双向比例调整的多约束数控加工进给率定制方法 | |
CN109991921A (zh) | 一种平顺b样条轨迹直接生成方法 | |
CN113345112B (zh) | 一种长骨骨折断面点云预处理及配准方法 | |
CN111880480B (zh) | 一种基于cnc铣床的铣刀切割路径生成方法及系统 | |
CN110458822B (zh) | 一种复杂曲面零件非接触式三维匹配检测方法 | |
CN113835397A (zh) | 基于b样条曲线和路径积分的线性数控加工路径平滑方法 | |
Zheng et al. | Robot path optimization for laser cladding forming | |
JP2527230B2 (ja) | 工具径路生成方法 | |
Lu et al. | A novel mathematical model for the accurate measurement of face gears by considering the geometric deviations of multiple teeth |
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 |