一种基于数值积分法的分层牙体预备轨迹优化方法

本发明专利涉及一种基于数值积分法的分层牙体预备轨迹优化方法，属于牙体预备轨迹优化。

背景技术：

1、龋齿作为人类三大重要防治疾病之一是导致牙体缺损的重要原因，严重的影响了人们的口腔健康。口腔修复是治疗牙体缺损的重要手段，牙体预备是修复过程中的必要治疗环节，是指医生对患龋齿处的牙齿进行定量去除并形成预期三维形状的操作过程。在传统的牙体预备过程中，需要依赖于医生手工的操作并结合丰富的临床经验对患龋牙齿进行大量的重复的微细调整。然而目前无法满足当前极大的牙体预备需求量，因此需要引入机器人和辅助软件来代替或辅助医生高效的完成牙体预备工作；

2、在使用机器人进行牙体预备情况下，将轨迹分成多层能够有效减少牙体预备过程中牙体的受力，降低预备过程中牙体出现裂纹，以及本来存在隐裂纹的牙体出现裂纹扩展的风险；此外在牙体预备过程中，一方面由于功能尖斜面凹凸不平，曲率变化很大，另一方面由于车针中轴面外轮廓具有二次函数的形状，导致预备体表面会具有一定的残余牙体，残余牙体的存在会影响后续的牙冠安装以及牙冠对牙体的保护；本发明通过控制牙体预备轨迹的密度，从而在最小预备轨迹密度的同时，使残余牙体在评价截面上的轮廓的面积达到要求的阙值，实现对分层预备轨迹的优化，在保证牙体预备的精度同时，能够提高牙体预备的效率，为牙体预备机器人实现临床应用提供理论基础；

3、综上，目前牙体预备轨迹优化技术领域亟待一种功能尖斜面分层牙体预备轨迹优化方法。

技术实现思路

1、针对上述问题，本发明提出了一种基于数值积分法的分层牙体预备轨迹优化方法，解决目前牙体预备轨迹优化技术领域中缺少针对一种功能尖斜面分层牙体预备轨迹优化方法，实现了对功能尖斜面分层牙体预备轨迹的优化。

2、一种基于数值积分法的分层牙体预备轨迹优化方法，其特征在于：所述方法的具体实现过程为：

3、步骤一、输入车针行走路径信息集：

4、在三维软件中设计得到三维理论牙齿预备体，以右手定则建立o-xyz三维预备轨迹标定坐标系op，设定总体预备深度，用符号h表示，总体预备深度表示牙冠的最高点与预备体最高点之间的距离；设定每层预备深度为d；a为分层预备轨迹的层数，a＝[h/d]，并使用平行于yoz平面的c个评价截面x＝i,jx在a层预备轨迹上截取刀触点信息，每个评价截面截取b条路径；其中每个评价截面第k个刀触点连接成一条路径，其中1≤k≤c；计算并输入具有a层b条路径c个评价截面的刀触点信息集：

5、

6、其中，i,jpk为第i层第j条路径的第k个评价截面上的刀触点，1≤i≤a，1≤j≤b，i，j为整数，相对于预备轨迹标定坐标系op的位置信息为(i,jxk,i,jyk,i,jzk)，其中i,jxk是刀触点i,jpk相对于预备轨迹标定坐标系op中的x坐标，i,jyk是刀触点i,jpk相对于预备轨迹标定坐标系op中的y坐标，i,jzk是刀触点i,jpk相对于预备轨迹标定坐标系op中的z坐标，初始化i，j，k，令i＝1，j＝1，k＝1；

7、步骤二、计算车针中轴面外轮廓在刀触点i,jpk处往x＝i,jxk平面投影的投影方程：

8、定义切向向量，用符号表示，切向向量表示刀触点i,jpk指向刀触点i,jpk+1的单位向量，其位置信息为设定i,jdpk,k+1为刀触点i,jpk指向刀触点i,jpk+1的距离，其中：

9、

10、

11、a)判断k＝1是否成立：

12、若k＝1成立，则以刀触点i,jp1为原点，为x轴，和的向量积的向量的单位向量为y轴，以x轴和y轴的向量积的向量的单位向量为z轴，建立刀触点坐标系i,jh1，其中x轴用i,juk表示，坐标形式为(i,juxk,i,juyk,i,juzk)，y轴用i,jvk表示，坐标形式为(i,jvxk,i,jvyk,i,jvzk)，z轴用i,jwk表示，坐标形式为(i,jwxk,i,jwyk,i,jwzk)，跳转至步骤二c)；

13、若k＝1不成立，则跳转至步骤二b)；

14、b)判断k＝c是否成立：

15、若k＝c成立，则以刀触点i,jpc为原点，为x轴，和的向量积的向量的单位向量为y轴，以x轴和y轴的向量积的向量的单位向量为z轴，建立刀触点坐标系i,jhc，其中x轴用i,juk表示，坐标形式为(i,juxk,i,juyk,i,juzk)，y轴用i,jvk表示，坐标形式为(i,jvxk,i,jvyk,i,jvzk)，z轴用i,jwk表示，坐标形式为(i,jwxk,i,jwyk,i,jwzk)，跳转至步骤二c)；

16、若k＝c不成立，则以刀触点i,jpk为原点，为x轴，和的向量积的向量的单位向量为y轴，以x轴和y轴的向量积的向量的单位向量为z轴，建立刀触点坐标系i,jhk，其中x轴用i,juk表示，坐标形式为(i,juxk,i,juyk,i,juzk)，y轴用i,jvk表示，坐标形式为(i,jvxk,i,jvyk,i,jvzk)，z轴用i,jwk表示，坐标形式为(i,jwxk,i,jwyk,i,jwzk)，跳转至步骤二c)；

17、c)设定在三维预备轨迹标定坐标系u下的车针中轴面外轮廓数学表达式为f(x0,y0,z0)，f(x0,y0,z0)是关于x0,y0,z0的函数，(x0,y0,z0)是f(x0,y0,z0)上的任意点，具体表达式为设定三维预备轨迹标定坐标系u下f(x0,y0,z0)的坐标变换到刀触点坐标系i,jhk后，在三维预备轨迹标定坐标系u下的车针中轴面外轮廓数学表达式为是关于的函数，是上的任意点；因此上的点与f(x0,y0,z0)上的点(x0,y0,z0)之间存在其中r为车针坐标变换矩阵，具体为由得设定则可以得到方程将代入方程得的具体表达式为

18、设定在x＝i,jxk平面投影的表达式为是关于i,jxk，的函数，是上的任意点；由方程得

19、将代入方程：

20、得到方程：

21、将方程：

22、在x＝i,jxk平面进行投影得到刀触点i,jpk处的投影表达式

23、

24、d)判断1≤k＜c是否成立：

25、若1≤k＜c成立，则k＝k+1跳转至步骤二a)；

26、若1≤k＜c不成立，则跳转至步骤二e)；

27、e)判断1≤j＜b是否成立：

28、若1≤j＜b成立，则j＝j+1，k＝1，跳转至步骤二a)；

29、若1≤j＜b不成立，则跳转至步骤二f)；

30、f)判断1≤i＜a是否成立：

31、若1≤i＜a成立，则i＝i+1，j＝1，k＝1，跳转至步骤二a)；

32、若1≤i＜a不成立，则i＝1，j＝1，k＝1，跳转至步骤三；

33、步骤三、计算刀触点i,jpk和刀触点i,j+1pk间的z轴坐标信息集i,j,j+1zk：

34、设定集合i,j,j+1zk，i,j,j+1zk是刀触点z轴坐标信息集，具体信息为：其中i,jzk和i,j+1zk是刀触点i,jpk和刀触点i,j+1pk相对于预备轨迹标定坐标系op中的z坐标，设定i,j,j+1zk(t)，i,j,j+1zk(t)表示i,j,j+1zk中的第t个值，其中1≤t≤n+1；设定集合i,j,j+1yk，i,j,j+1yk是将i,jyk与i,j+1yk之间进行n等分的刀触点y轴坐标信息集，具体信息为：其中i,jyk和i,j+1yk是刀触点i,jpk和刀触点i,j+1pk相对于预备轨迹标定坐标系op中的y坐标；是i,jyk与i,j+1yk之间进行n等分的第m个等分刀触点y轴坐标，设定i,j,j+1yk(t)，i,j,j+1yk(t)表示i,j,j+1yk中的第t个值；

35、已知中的i,jxk是刀触点i,jpk相对于三维预备轨迹标定坐标系op中的x坐标，将集合i,j,j+1yk中的i,j,j+1yk(t)赋值给投影方程中的得到投影方程：中的3个变量已知2个变量，此时投影方程是关于的2次函数，设定投影方程的两个解信息集为i,j,j+1zfk,1，i,j,j+1zsk,1，其中：设定i,j,j+1zfk,1(t)，i,j,j+1zsk,1(t)，i,j,j+1zfk,1(t)表示i,j,j+1zfk,1中的第t个值，i,j,j+1zsk,1(t)表示i,j,j+1zsk,1中的第t个值；

36、同理已知中的i,j+1xk是刀触点i,j+1pk相对于三维预备轨迹标定坐标系op中的x坐标，将集合i,j,j+1yk中的i,j,j+1yk(t)赋值给投影方程中的得到投影方程中的3个变量已知2个变量，此时投影方程是关于的2次函数，设定投影方程的两个解信息集为i,j,j+1zfk,2，i,j,j+1zsk,2，其中：设定i,j,j+1zfk,2(t)，i,j,j+1zsk,2(t)，i,j,j+1zfk,2(t)表示i,j,j+1zfk,2中的第t个值，i,j,j+1zsk,2(t)表示i,j,j+1zsk,2中的第t个值；

37、定义截面向量，用符号表示，截面向量表示由刀触点i,jpk指向i,j+1pk的单位向量，其三维坐标形式为其中：

38、定义倾斜度，用符号θ表示，倾斜度表示截面向量与三维预备轨迹标定坐标系op的xoy平面的单位法向量(0,0,1)之间的夹角，其中由数量积公式可得：其中代表(0,0,1)与进行数量积运算；设定预备体第i+1层理论高度信息集，用符号i,j,j+1zhk表示，其中：由于则有：

39、设定i,j,j+1zhk(t)，i,j,j+1zhk(t)表示i,j,j+1zhk中的第t个值；设定t的初值为1；

40、a)判断方程是否有解：

41、若方程有解，则将得到的两个解赋值给i,j,j+1zfk,1(t)，i,j,j+1zsk,1(t)，跳转至步骤三b)；

42、若方程无解，则i,j,j+1zfk,1(t)＝i,j,j+1zhk(t)，i,j,j+1zsk,1(t)＝i,j,j+1zhk(t)，跳转至步骤三b)；

43、b)判断方程是否有解：

44、若方程有解，则将得到的两个解赋值给i,j,j+1zfk,2(t)，i,j,j+1zsk,2(t)，将i,j,j+1zhk(t)中的最小值赋值给i,j,j+1zk(t)，跳转至步骤三c)；

45、若方程无解，则i,j,j+1zfk,2(t)＝i,j,j+1zhk(t)，i,j,j+1zsk,2(t)＝i,j,j+1zhk(t)，将i,j,j+1zhk(t)中的最小值赋值给i,j,j+1zk(t)，跳转至步骤三c)；

46、c)判断1≤t＜n+1是否成立：

47、若1≤t＜n+1成立，则t＝t+1跳转至步骤三a)；

48、若1≤t＜n+1不成立，则跳转至步骤三d)；

49、d)判断1≤k＜c是否成立：

50、若1≤k＜c成立，则k＝k+1，t＝1，跳转至步骤三a)；

51、若1≤k＜c不成立，则跳转至步骤三e)；

52、e)判断1≤j＜b-1是否成立：

53、若1≤j＜b-1成立，则j＝j+1，k＝1，t＝1，跳转至步骤三a)；

54、若1≤j＜b-1不成立，则跳转至步骤三f)；

55、f)判断1≤i＜a是否成立：

56、若1≤i＜a成立，则i＝i+1，j＝1，k＝1，t＝1，跳转至步骤三a)；

57、若1≤i＜a不成立，则i＝1，j＝1，k＝1，t＝1，跳转至步骤四；

58、步骤四、计算刀触点i,jpk与刀触点i,j+1pk之间的残余牙体在x＝i,jx评价截面上的轮廓的残差面积i,j,j+1sk：

59、定义残差面积，用符号表示i,j,j+1sk，残差面积表示进行牙体预备后，刀触点i,jpk与刀触点i,j+1pk之间的残余牙体在评价截面上的轮廓的面积，设定i,j,j+1sk的初始值为0；

60、a)判断1≤t＜n+1是否成立：

61、若1≤t＜n+1成立，则t＝t+1，跳转至步骤四a)；

62、若1≤t＜n+1不成立，则i,j,j+1sk＝i,j,j+1sk-i,j,j+1s梯，其中跳转至步骤四b)；

63、b)判断1≤k＜c是否成立：

64、若1≤k＜c成立，则k＝k+1，t＝1，跳转至步骤四a)；

65、若1≤k＜c不成立，则跳转至步骤四c)；

66、c)判断1≤j＜b-1是否成立：

67、若1≤j＜b-1成立，则j＝j+1，k＝1，t＝1，跳转至步骤四a)；

68、若1≤j＜b-1不成立，则跳转至步骤四d)；

69、d)判断1≤i＜a是否成立：

70、若1≤i＜a成立，则i＝i+1，j＝1，k＝1，t＝1，跳转至步骤四a)；

71、若1≤i＜a不成立，则i＝1，j＝1，k＝1，t＝1，跳转至步骤五；

72、步骤五、判断分层牙体预备轨迹是否合格：

73、设定残差面积最大允许限度smax；

74、a)判断i,j,j+1sk＞smax是否成立：

75、若i,j,j+1sk＞smax成立，则分层牙体预备轨迹不合格，b＝b+1，跳转至步骤一；

76、若i,j,j+1sk＞smax不成立，则继续进行优化，跳转至步骤五b)；

77、b)判断1≤k＜c是否成立：

78、若1≤k＜c成立，则k＝k+1，跳转至步骤五a)；

79、若1≤k＜c不成立，则跳转至步骤五c)；

80、c)判断1≤j＜b-1是否成立：

81、若1≤j＜b-1成立，则j＝j+1，k＝1，跳转至步骤五a)；

82、若1≤j＜b-1不成立，则跳转至步骤五d)；

83、d)判断1≤i＜a是否成立：

84、若1≤i＜a成立，则i＝i+1，j＝1，k＝1，跳转至步骤五a)；

85、若1≤i＜a不成立，则分层牙体预备轨迹合格，优化结束。

86、本发明的有益效果为：

87、1、本发明在对机器人辅助牙体预备分层预备轨迹进行优化时，提出根据不同车针姿态计算相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk的概念，对相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk进行了量化描述，并将残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk作为优化预备轨迹密度的参数，通过设定残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk的上限值smax，可以对相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积进行约束；

88、2、本发明针对判断相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk，以不同刀触点为原点，建立刀触点坐标系i,jhk，通过坐标变换，在评价平面x＝i,jxk中计算车针不同姿态的中轴面的外轮廓在三维预备轨迹标定坐标系op的表达式，方便计算相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk；

89、3、本发明针对判断相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk的大小，引入每层预备深度d的概念，在计算相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的材料残差面积i,j,j+1sk时考虑预备深度d，从而准确计算相邻预备轨迹上对应刀触点之间的残余牙体在评价截面上的轮廓的残差面积；

90、4、本发明针对判断相邻预备轨迹上同一评价截面刀触点i,jpk与i,j+1pk之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk，引入车针中轴面的外轮廓的表达式为一种常见的二次函数，当使用不同车针时，可以根据需要改变车针中轴面的外轮廓的表达式，进而计算相邻预备轨迹上对应刀触点之间的残余牙体在评价截面上的轮廓的残差面积i,j,j+1sk从而优化预备轨迹的密度；

91、5、同本发明人同日申报的发明专利《一种基于高宽比法的分层预备轨迹优化方法》相比，虽然两种方法都是用于优化机器人辅助牙体预备的预备轨迹，但是对于《一种基于高宽比法的分层预备轨迹优化方法》中提及的优化方法是针对相邻预备轨迹上对应刀触点之间的残余牙体在评价截面上的轮廓的尖锐程度，来优化预备轨迹的密度，此外，《一种基于高宽比法的分层预备轨迹优化方法》更加关注的是优化的过程中预备体的光滑程度，计算速度快，但其精准性没有本方法高；本方法是针对相邻预备轨迹上对应刀触点之间的残余牙体在评价截面上的轮廓的面积，来优化预备轨迹的密度，更加关注的是优化的过程中预备体残余牙体的面积，虽然计算速度慢，但精准性高，因此本方法的提出与《一种基于高宽比法的分层预备轨迹优化方法》互相补偿，进而完善了采用数值法优化机器人辅助牙体预备的轨迹的系列方法。