NURBS插补快速求值算法*

孔祥洪　李迪　焦青松

(华南理工大学 机械与汽车工程学院, 广东 广州510640)



NURBS插补快速求值算法*

孔祥洪李迪焦青松

(华南理工大学 机械与汽车工程学院, 广东 广州510640)

摘要：当对非均匀有理B样条(NURBS)曲线进行高密度插值时,运用分段幂函数方法对基函数进行求值的效率远高于传统的de-Boor算法.为此,文中从NURBS插补计算的特点出发,结合de-Boor递推计算规律,设计了NURBS插补快速求值算法.首先采用该算法计算NURBS在各节点区间的基函数显式方程,再运用显式方程进行NURBS插补点求值,并设计相应的NURBS曲线插补器.复杂NURBS曲线的铣削加工实验结果表明,该算法能够有效地缩减NURBS曲线插补求值的计算耗时,提高插补计算的实时性.

关键词：NURBS;基函数;快速算法;参数插补

在手机工业、航空工业以及模具工业中,出于对特殊形状、特殊功能或高加工质量的需求,自由曲面零件的加工制造已成为常态.

加工自由曲面零件通常是一个复杂且费时费力的过程.越来越多的零件加工开始采用非均匀有理B样条(NURBS)插补来规避传统的计算机辅助制造(CAM)系统在自由曲面和自由曲线加工方面的不足.然而,目前大量的NURBS插补器的研究主要集中于三阶NURBS曲线插补[1-7].过于沉重的计算负荷成为限制高阶(四阶以上)NURBS曲线应用于加工工业的重要因素.高阶NURBS曲线能够提供更高的拟合精度、更优秀的连续性,适合有高加工精度需求的场合.目前的NURBS插补算法都只能在精度与实时性之间进行取舍[8].大部分的插补算法都倾向于为了获取更好的实时性而在插补精度方面做出适当的妥协.

高阶NURBS曲线插补的高负荷计算量主要来源于其求值求导计算,而在NURBS求值求导计算中最耗时的是NURBS基函数的求值计算.NURBS基函数的定义方法主要有采用截尾幂函数的差商定义、开花定义、递推公式定义等[9-11].文献[12-13]基于差商定义提出了NURBS曲线曲面的显示矩阵表示方法及其算法.文献[11]基于递推公式定义提出了基函数递推算法(de-Boor算法),由于该算法易于在计算机中实现而成为了NURBS基函数求值最为常用的算法.

运用分段幂函数对NURBS基函数进行求值的算法也渐渐引起人们的关注,当要对NURBS曲线进行高密度插值时,该算法要比传统de-Boor算法高效得多[13].NURBS曲面曲线矩阵表示正是源于该算法的研究.但有关系数矩阵的显式表示研究大多集中于推导计算,很少应用于插补算法中[14-15].王国勋等[16]提出的快速求值求导插补算法也局限于对NURBS曲线的某个节点区间进行插补研究.

运用系数矩阵推导计算NURBS基函数的分段幂函数,由于其计算过程过于复杂而难以应用于NURBS插补算法中.为了简化推导过程,提高计算效率,文中设计了基于前向路径增益求和方法的NURBS插补快速求值算法:通过计算基函数Ni,p(u)与零阶基函数Ni,0(u)之间的前向路径增益之和来快速计算Ni,p(u)在其对应的节点区间的分段幂函数表达式;应用该算法设计了相应的NURBS插补器,并通过加工实验验证该算法的效率.

1NURBS表达式和de-Boor算法

1.1　NURBS曲线表达式

根据NURBS曲线的数学定义[17],对于给定的节点向量U=(u0,u1,…,um)、权重向量W=(w0,w1,…,wn)以及控制点向量P=(P0,P1,…,Pn),p阶NURBS曲线公式可以表示为

(1)

k阶NURBS导矢计算公式为

(2)

式中,A(k)和W(i)可由B样条曲线的导矢公式

(3)

求得.

1.2　NURBS基函数的de-Boor递推算法

de-Boor算法的递推计算公式为

(4)

根据式(4),在参数区间u∈[ui,ui+1)内,零阶基函数有且仅有Ni,0(u)取值非0;根据式(4)中的递推关系,p阶基函数只有p+1个基函数取值为非0.该性质被称为基函数的局部支撑性.因此p阶NURBS曲线上任意点的求值都只需要计算p+1项非0的基函数的值.de-Boor算法是运用基函数局部支撑性的算法.在u∈[ui,ui+1)区间内,只有基函数Ni,p(u),Ni-1,p(u),…,Ni-p,p(u)取值非0,则运用de-Boor算法对以上基函数进行求值计算的过程形如金字塔的结构,如图1所示.

图1　de-Boor算法计算B样条基函数的过程Fig.1　Calculation process of B-spline basis function using de-Boor algorithm

2基函数的分段幂函数快速计算方法

如前所述,运用de-Boor算法计算NURBS基函数需要从零阶基函数开始依次向上进行递推计算.而运用式(1)计算NURBS曲线上的点只需要p阶基函数,即运用de-Boor算法计算得到的0~p-1阶基函数的值均为计算p阶基函数而产生的中间值.文中设计的算法是运用基函数的分段幂函数直接计算p阶基函数的值,避免了0~p-1阶基函数求值的计算,从而提高计算效率.

通常是使用差商法[13]计算NURBS基函数的分段幂函数,进而得到系数矩阵.但该方法的计算推导过程较为复杂且存在高阶误差.为了简化NURBS基函数的分段幂函数的计算过程,提高计算效率,文中设计了一种快速计算基函数的分段幂函数的方法,即将式(4)改写为

(5)

图2给出了一阶基函数Ni-1,1、二阶基函数Ni-1,2和三阶基函数Ni-1,3、Ni-2,3的递推计算过程.运用式(5)计算得到这4个基函数的取值分别为

图 2　B样条基函数的递推计算过程Fig.2　Recursive calculation process of B-spline basis function

根据图2所示的基函数递推过程,以上列举的基函数递推计算过程存在如下规律:假设图1中所有基函数值被看作计算节点,节点间的连线为通道,连线的箭头指示通道的前进方向,节点间的系数(fi,p或gi,p)为通道增益,将从初始节点Ni,0沿着通道前进方向连接到任意节点Ni,p所经过的通道定义为前向路径,所经过通道的增益的积定义为路径增益,则Ni,p的值等于所有不重复的路径增益之和,即

(6)

式中,L(j)为第j条前向路径的增益值,R为不重复的前向路径总数.

为了验证式(6)计算方法的可行性,以图2为例进行说明.如图2(a)所示,从Ni,0到Ni-1,1只有一条前向路径,其增益为gi,1,Ni-1,1的路径增益之和为gi,1.在图2(b)中,Ni,0到Ni-1,2有两条不重复的前向路径,其增益分别为fi,1gi,2和gi,1fi-1,2,Ni-1,2的路径增益之和为fi,1gi,2+gi,1fi-1,2.在图2(c)中,Ni,0到Ni-1,3有3条不重复的前向路径,它们的增益分别为fi,1fi,2gi,3、fi,1gi,2fi-1,3和gi,1fi-1,2gi-1,3,Ni-1,3的路径增益之和为fi,1fi,2gi,3+fi,1gi,2fi-1,3+gi,1fi-1,2gi-1,3.显然,以上结果与运用式(5)递推计算得到的结果完全一致.

现用数学归纳法证明式(6).

证明当p=1时,根据图1可以知道,Ni,1(u)和Ni-1,1(u)的前向路径数为1,其通道增益分别为fi,1和gi,1,根据式(6)计算可以得到Ni,1(u)=fi,1,Ni-1,1(u)=gi,1.

运用de-Boor递推公式(式(4))计算Ni,1(u)和Ni-1,1(u)的值:

显然，当p=1时式(6)成立.

假设当p=n-1时式(6)成立,则有

根据式(4)有

Ni,n(u)=fi,pNi,n-1(u)+gi+1,pNi+1,n-1(u)=

(7)

式中,∑fi,pL1(j)为Ni,0节点经由Ni,n-1节点连接到Ni,n节点的前向路径的增益之和,∑gi+1,pL2( j )为Ni,0节点经由Ni+1,n-1节点连接到Ni,n节点的前向路径的增益之和,两者的路径不重复且所有从Ni,0节点连接至Ni,n节点的前向路径必然经过以上两节点,故∑fi,pL1(j)与∑gi+1,pL2(j)的和等于所有Ni,0节点连接到Ni,n节点的不重复路径增益之和,即当p=n时式(6)的结论也成立,证毕.

根据式(5)可以确定任意阶NURBS基函数与Ni,0的递推关系.图3给出了一阶到五阶的NURBS基函数递推关系.结合图3所示的NURBS基函数递推关系,运用式(6)能够简单、直观地计算得到五阶以内的NURBS各参数节点区间的基函数的分段幂函数,继之能够快速、高效地计算得到基函数的值,完成NURBS插补求值计算.若需要计算五阶以上的NURBS的分段幂函数,只需要运用式(5),把图3中的递推关系依次向上往高阶递推即可.

图3　五阶B样条基函数的递推计算过程Fig.3　Recursive calculation process of 5-order B-spline basisfunction

3NURBS曲线插补器设计

应用于NURBS曲线的插补器设计流程如图3所示.该NURBS曲线插补器主要由速度规划算法、插补算法及插补点求值算法组成.速度规划算法主要是根据设定的动力学模型规划插补过程中每个插补周期的插补速度.插补算法主要是根据当前的规划速度和当前样条参数值(uk)计算出下一插补点的参数值(uk+1).求值算法则要根据计算得到的样条参数值(u)计算出与之相对应的插补点的值.

运用式(6)计算得到的基函数分段幂函数表达式在NURBS的参数节点区间内是分段连续的.在NURBS曲线插补的过程中,样条参数逐渐增大并依次跨越NURBS曲线的各个参数节点.基于以上幂函数得到分段连续性,当样条参数跨越参数节点时,需要定位样条参数当前所处的节点区间,并重新计算基函数对应于该参数节点区间的显式方程.在样条参数跨越下一个参数节点前,计算NURBS插补点所需要的基函数值，均可以运用以上显式方程进行快速地计算.因此，在图4的NURBS曲线插补器中设计了判断环节,用于判断插补过程中样条参数的值是否跨越参数节点,并依次判断是否需要重新计算基函数在该节点区间的显式方程.

图4　NURBS曲线插补器流程图Fig.4　Flowchart of NURBS curve interpilator

4算法效率分析与实验

4.1　算法效率分析

对各阶基函数的求值计算在NURBS的求值求导计算过程中占用了绝大部分的计算耗时.通过分析每个NURBS插补点各阶基函数求值的计算规模，即能评估不同求值求导算法的计算效率.

为评估文中提出的快速求值求导算法的计算效率,本节设计了实例仿真对算法进行分析.用于仿真实验的NURBS曲线参数如下:阶数p=3;节点序列U={0,0,0,0,1,2,3,4,5,6,7,8,9,10,10,10,10};控制顶点集P={(0,0),(0.1,20.01),(0.2,20.02),(0.3,20.03),(20,40),(25,33),(30,20),(35,7),

(40,0),(54.8,18.98),(54.9,18.99),(55,19),(60,20)};权重向量w=(1,1,1,1,1,1,1,1,1,1,1,1,1).由以上参数确定的NURBS曲线如图5所示.

文中采用的仿真实例为三阶NURBS曲线,故在曲线插补计算中需要计算三阶基函数的显式方程.假设当前所处的参数节点区间为[ui,ui+1],运用式(6)表示的算法并结合图3所示的基函数递推关系,得到基函数的分段幂函数显式表达式的计算过程,如图6所示.根据NURBS的局部支撑性性质,

图5　NURBS仿真曲线及控制多边形Fig.5　NURBS simulation curve and control ploygon

图6　B样条基函数的分段幂函数的计算过程Fig.6　Calculation process of piecewise power function of B-spline basis function

在该节点区间内只有Ni,3、Ni-1,3、Ni-2,3和Ni-3,3的值为非0.图6中Ni,3和Ni-3,3分别只有一条前向路径,其增益分别为fi,1fi,2fi,3和gi,1gi-1,2gi-2,3;Ni-1,3和Ni-2,3各有3条前向路径,它们的增益分别为fi,1fi,2gi,3、fi,1gi,2fi-1,3、gi,1fi-1,2fi-1,3和fi,1fi,2gi,3、fi,1gi,2fi-1,3、fi,1gi,2fi-1,3、gi,1fi-1,2fi-1,3.根据式(6)可以求得基函数Ni,3、Ni-1,3、Ni-2,3和Ni-3,3的前向路径增益之和:

(8)

把式(5)中的fi,p和gi,p代入式(7),则有

(ui+1-ui-1)-1(ui+1-ui)-1,

(ui+1-ui-1)-1(ui+1-ui)-1,

以上等式即为在u∈[ui,ui+1)区间内的基函数分段幂函数的显式表达式.取图5所示曲线的u∈[1,2)计算该节点区间的基函数显式方程,该节点区间的节点索引值i=4,根据NURBS局部支撑性,在该节点区间不为0的基函数为N1,3、N2,3、N3,3和N4,3.计算显式方程所需参数节点值为ui-2=ui-1=0,ui=1,ui+1=2,ui+2=3,ui+3=4.

把以上参数节点的值代入基函数的分段幂函数方程,可以得到u∈[1,2)的基函数显式方程,即

运用以上基函数显式方程,可以通过样条参数u代入直接计算基函数的值.在NURBS插补点求值计算过程中,计算三阶基函数N1,3、N2,3、N3,3和N4,3的取值只需要进行12次乘法运算和12次加法运算.

为比较,按图1所示的计算方法,基于de-Boor算法进行基函数求值是一个逐阶向上递推的计算过程.根据递推公式(4),对于NURBS曲线任意一个p阶基函数的求值,运用de-Boor算法需要进行5p(p+1)/2次加法计算、p(p+1)次乘法计算和p(p+1)次除法计算,故计算p阶所有基函数时需要进行5p(p+1)2/2次加法计算、 p(p+1)2次乘法计算和p(p+1)2次除法计算.

文献[17]给出了优化的de-Boor算法.运用该优化算法能够减少传统de-Boor算法中存在的冗余计算,明显提高计算效率.3种算法在一次NURBS插补中计算基函数取值需要耗费的计算量比较如表1所示.

表13种算法的计算量对比

Table1Comparisonofarithmeticoperationamongthreealgorithms

运算传统de-Boor算法优化de-Boor算法文中算法加法1204512乘法481812除法48180

为了简化比较结果,假设乘法运算和除法运算的计算耗时与加法运算的计算耗时相等,则运用传统de-Boor算法在一次NURBS插补中需要进行216次运算,运用优化de-Boor算法需要进行81次运算,而文中算法只需要进行24次运算.文中算法所需计算耗时为传统de-Boor算法的1/9,相当于优化de-Boor算法的1/3.基于以上比较结果,文中算法的计算效率明显优于de-Boor算法.若是四阶以上的NURBS基函数求值,文中算法的优势将会更加明显.

4.2　加工实验及结果分析

如图7所示,加工实验由基于PC平台构建的运动控制系统完成.运动控制系统主要是由PC机、上位机编程软件、运动控制卡以及插补算法构成.上位机编程软件主要用于实现人机数据交互功能,方便地导入NURBS曲线、设置参数和监视加工状态.运动控制卡采用PCI接口与PC机实现数据交换,在板卡中采用了高性能DSP芯片实现插补算法,再通过板卡上的FPGA芯片采集IO信号和实现控制脉冲输出.

图7　加工实验Fig.7　Machining experiment

加工实验采用了一段较为复杂的三阶NURBS曲线(形如一只飞翔中的鸽子).NURBS插补器被构建于板卡的DSP芯片内.上位机软件采集的NURBS曲线参数直接传输至板卡上的DSP芯片中进行插补运算,并结合FPGA芯片完成NURBS的加工控制.加工实验采用的运动控制参数如下:最大速度为50 mm/s,最大加速度为100 mm/s2,最大加加速度为5 000 mm/s3,插补周期为1 ms,前瞻周期数为2 000,弦高容差为0.001 mm,速度浮动容差为0.005%.

如图4所示,完整的NURBS插补算法主要由运动规划、参数插补、插补点求值3个步骤组成.在运动规划方面,加工试验采用了具有前瞻功能的自适应运动规划算法及S型加减速策略.在参数插补求值方面,加工试验中采用了重构映射插补算法[1].

运用DSP内部的定时器能够监测插补过程的总耗时.传统de-Boor算法[11]、差商算法[12]、系数矩阵算法[14]、优化de-Boor算法[17]及文中算法的插补总耗时分别为35.636、30.691、26.706、26.592、20.628 s.如前所述,插补器中除了插补求值算法外,还有运动规划算法、前瞻算法、曲率自适应算法等.插补总耗时是插补求值算法耗时加上其他算法耗时的结果.由于实验条件所限，无法在加工实验过程中分别检测得到各算法在插补过程中的耗时，但以前面的总插补耗时结果进行粗略评估,文中算法的耗时依然明显少于其他算法,算法效率以及实时性都有明显的优势.若将前面的插补总耗时减去其他算法的耗时,则文中算法的计算效率相对于其他算法具有更加明显的优势.

5结论

文中从NURBS插补计算的特点出发,结合de-Boor递推计算规律,设计了快速获取NURBS基函数的分段多项式解析式的算法.该算法获得的多项式解析式能直接用于样条函数求值,避免了繁杂的递推计算,大幅度削减所需的计算耗时.文中算法通过有效控制样条基函数的计算量来实现NURBS曲线曲面上插补点及其导数的快速求值,增强了插补器的实时性.文中基于所提快速求值算法设计开发了NURBS插补运动控制系统,并完成复杂NURBS曲线的铣削加工实验,结果表明,该算法能够有效地缩减NURBS曲线插补求值的计算耗时,提高插补计算的实时性.文中所提算法不仅适用于提高样条基函数的计算效率,还可以进一步导出NURBS曲线曲面的基函数系数矩阵,拓宽NURBS在理论和应用中的研究空间.

参考文献:

[1]LIU Min,HUANG Yu,YIN Ling.Development and implementation of a NURBS interpolator with smooth feedrate scheduling for CNC machine tools [J].International Journal of Machine Tools & Manufacture,2014,87:1-15.

[2]ZHAO Huan,ZHU LiMin,DING Han.A parametric interpolator with minimal feed fluctuation for CNC machine tools using arc-length compensation and feedback correction [J].International Journal of Machine Tools & Manufacture,2013,75:1-8.

[3]CHENG C W,TSAI M C.Real-time variable feedrate NURBS curve interpolator for CNC machining [J].International Journal of Advanced Manufacturing Technology,2004,23(11):865- 873.

[4]SUN Yuwen,ZHAO Yang,BAO Yurong,et al.A novel adaptive-feedrate interpolation method for NURBS tool path with drive constraints [J].International Journal of Machine Tools & Manufacture,2014,77:74-81.

[5]FANG Jing-Jing,HUNG Chia-Lien.An improved parameterization method for B-spline curve and surface interpolation [J].Computer-Aided Design,2013,45(6):1005-1028.

[6]GRECO L,CUOMO M.B-spline interpolation of Kirchhoff-love space rods [J].Computer Methods in Applied Mechanics and Engineering,2013,56:251-269.

[7]王世勇,李迪.NURBS图形激光雕刻算法及其嵌入式实现 [J].华南理工大学学报(自然科学版),2010,38(6):112-117.

WANG Shi-yong,LI Di.Algorithm and embedded implementation of laser marking of NURBS images [J].Journal of South China University of Technology(Natural Science Edition),2010,38(6):112-117.

[8]LEI W T,WANG S B.Robust real-time NURBS path interpolators [J].International Journal of Machine Tools & Manufacture,2009,49:625- 633.

[9]CURRY H B,SCHOENBERG I J.On spline distributions and their limits:the Pólya distribution functions [J].Bulletin of the American Mathematical Society,1947,53:109.

[10]RANSHAW L.Blossoming:a connect-the-dots approach to splines [R].Palo Alto:Digital,Systems Research Center,

1987.

[11]DE BOOR C.On calculating with B-splines [M].New York:Springer-Verlag,1978.

[12]SCHUMKER L.Spline functions:basic theory [M].New York:John Wiley & Sons,1981.

[13]潘日晶.NURBS曲线曲面的显式矩阵表示及其算法 [J].计算机学报,2001,24(4):358-366.

PAN Ri-jing.Explicit matrix representation for NURBS curves and surfaces and its algorithm [J].Chinese Journal of Computers,2001,24(4):358-366.

[14]CHOI B K,YOO W S,LEE C S.Matrix representation for NURBS curves and surfaces [J].Computer-Aided Design,1990,22(4):235-240.

[15]王学福,孙家广,秦开怀.NURBS的符号矩阵表示及其应用 [J].计算机学报,1993,16(1):29-34.

WANG Xue-fu,SUN Jia-guang,QIN Kai-huai.Symbolic matrix representation of NURBS and its applications [J].Chinese Journal of Computers,1993,16(1):29-34.

[16]王国勋,舒启林,王军,等.NURBS 直接插补技术中快速求值求导算法 [J].东北大学学报(自然科学版),2012,33(7):1021-1024.

WANG Guo-xun,SHU Qi-lin,WANG Jun,et al.Algorithm of fast evaluation and derivation for the technique of NURBS direct interpolation [J].Journal of Northeas-tern University(Natural Science),2012,33(7):1021-1024.

[17]PIEGL L,TILLER W.The NURBS book [M].2nd ed.Berlin:Springer,1997.

An Efficient Evaluation Algorithm for NURBS Interpolation

KONGXiang-hongLIDiJIAOQing-song

(School of Mechanical and Automotive Engineering, South China University of Technology, Guangzhou 510640, Guangdong, China)

Abstract:In the high-density interpolation of non-uniform rational B-spline (NURBS) curves, using the piecewise power function method to evaluate the B-spline basis function consumes much less computing time than the traditional de-Boor algorithm. Therefore, on the basis of the characteristic of the NURBS interpolation, an efficient eva-luation algorithm for the NURBS interpolation is proposed by drawing on the recursive calculation laws of the de-Boor algorithm. First, the proposed algorithm is used to deduce the explicit equations of the B-spline basis function in each spline parameter knot interval. Then, NURBS interpolation points are evaluated by using explicit equations, and a corresponding NURBS curve interpolator is designed. The results of the milling experiment with complex NURBS curves show that the proposed algorithm can effectively reduce the computing time of the NURBS curve interpolation and can improve the real-time performance of NURBS interpolators.

Key words:NURBS; spline function; efficient algorithm; parametric interpolation

doi：10.3969/j.issn.1000-565X.2016.01.013

中图分类号：TP391

作者简介：孔祥洪(1984-),男,博士生,主要从事嵌入式运动控制算法研究.E-mail:6382026@qq.com

*基金项目：国家“863”计划项目(2012AA040909)

收稿日期：2015-05-07

文章编号：1000-565X(2016)01- 0085- 08

Foundation item： Supported by the National High-Tech R & D Program of China(2012AA040909)