用matlab实现三次-NURBS插值曲线

用matlab实现三次NURBS插值曲线作者:大漠孤狼发表于matlab乐园(kyxt.com)作者:这是我在大学时做大学生研究计划时写的,当时刚学会matlab,编写了这个程序,用了很多循环,效率不高.当时我并不清楚循环是matlab的弱点,等明白了,也不做这方面的工作了,也就懒的去改写了.如果谁需要用,就自己改吧.算法也有一些问题,我就不多说了,自己看吧一、三次NURBS插值算法给定平面控制顶点di(i=1,2,…n)及对应的权因子i(i=1,2,…n)，可定义一条三次NURBS曲线。先对控制顶点进行参数化，得一矢量：U=[u0,u1,u2…,un+4]则三次NURBS曲线的分段方程形式为：33,33,1)()()(kkiiikkiiikuNuNduPu[uk+3,uk+4]k=0,1,2,…,n-3(1)首先证明一条性质：若三点dk,dk+1,dk+2共线，且满足)()()()(33,2333,233,2333,1kkkkkkkkkkkkkkkuNuNduNduNd则三次NURBS曲线插值点dk+1.证明：由（1）式可得：)()()()()()()(33,2233,1133,233,22133,1133,3kkkkkkkkkkkkkkkkkkkkkkkuNuNuNduNduNduNuP)()()()()()()(43,3343,2243,11343,23243,22143,114kkkkkkkkkkkkkkkkkkkkkkkuNuNuNduNduNduNuP由以知可得)()()()()()(33,2233,1133,233,22133,1133,1kkkkkkkkkkkkkkkkkkkkkkuNuNuNduNduNduNd故性质得证。下面构造三次NURBS插值曲线；取一型值点列V0,V1,V2,…Vn,用下列方法决定各点Vi处的切矢Ti:设顺序五个数据点Vi-2,Vi-1,…,Vi+2,在中间数据点Vi的切矢ti被局部定义为ti=(1-a)Δpi-1+aΔpi其中，Δpi=pi+1-pia=BAA,A=12iippB=1iipp在非周期情况下，首末各两数据点处的切矢可采用贝塞尔条件确定，即用如下差分矢量Δp-1=2Δp0－Δp1Δp-2=2Δp-1－Δp0Δpn=2Δpn－1－Δpn－2Δpn+1=2Δpn－Δpn－1仍按前述公式确定t0,t1,tn-1,tn。下面，我们令d3i=pi(i=1…n)然后对数据点d3i进行参数化u0=u1=u2=u3=0njjjijjjiddddu1)1(3311)1(333(i=2,3,…，n)u3n+1=u3n+2=u3n+3=u3n+4=1在区间[u3i,u3(i+1)]插入节点u3i+1,u3i+2u3i+1=)1(333132iiuu)1(33233231iiiuuu即得节点矢量。下面在点d3i,d3(i+1)(i=0,1,…,n-1)之间插入两个辅助控制顶点d3i+1,d3i+2,构成新控制顶点组，得分段三次NURBS方程。我们按下列原则判断两型值点pi(d3i),pi+1(d3(i+1))之间的曲线段是否有拐点出现：如果矢量乘积ai-1×ai与ai×ai+1同向，则pi与pi+1之间无拐点，此时两切线Li,Li+1之间有一交点Vi*=Vi+iiiiitttta111否则，Vi与Vi+1之间有一拐点，记Vi*=（Vi+Vi+1）／2令T0=0*0VVTi=min{iiVV*,*1iiVV}(i=1,2,3……,n-1)Tn=*1nnVV我们按下列方法决定Vi与Vi+1之间的辅助点d3i+1,d3i+2讨论在区间[u3i+2,u3i+3]上的曲线段23133,23133,)()(iikkkiikkkkuNuNduPu∈[u3i+2,u3i+3]令d3i=Vi,i=0,1,2,……,nd3i+1=Vi+iiiiiittuN)(233,1313(0i)(233,1313iiiuNTi)d3i+2=Vi-11533,23231)(iiiiiittuN(0i+1)(533,2323iiiuNTi+1)由此三式得)()()()(233,1313233,131313233,13133233,13133iiiiiiiiiiiiiiiuNuNduNduNd由性质可得P(u3i+2)=d3i,即此曲线插值点d3i(Vi).二、基函数计算及整体表示1，基函数根据伯恩斯坦多项式可计算Ni,3(u)=)()())(()(0,11220,1232uNuuuuuNuuuuuuuuuuiiiiiiiiiiiii+)()())(())((0,21330,11211333uNuuuuuNuuuuuuuuuuuuiiiiiiiiiiiiii+)()())(())((0,22330,1121131414uNuuuuuNuuuuuuuuuuuuiiiiiiiiiiiiii+)()())(()(0,33440,2232241424uNuuuuuNuuuuuuuuuuiiiiiiiiiiiii整理后得：①Ni+1,0的系数(-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))-1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))-1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i)))*u^3+(1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(4+i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(1+i)+2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))*u(i)+1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))*u(2+i)+2/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(1+i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(3+i))*u^2+(-1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))*u(i)^2-2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))*u(i)*u(2+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(i)*u(1+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(1+i)^2-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(i)*u(3+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(3+i)*u(1+i)-2/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(4+i)*u(1+i))*u+1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))*u(i)^2*u(2+i)+1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(4+i)*u(1+i)^2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(i)*u(3+i)*u(1+i)②Ni+2,0(u)的系数(1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))+1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i)))*u^3+(-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)-1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(2+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(i)-2/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)-2/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(4+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(1+i))*u^2+(1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(4+i)^2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)^2+1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(1+i)+1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(3+i)+2/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(i)*u(3+i)+2/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(4+i)*u(2+i)+1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(1+i)*u(3+i))*u-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(1+i)*u(3+i)-1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(4+i)^2*u(2+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(i)*u(3+i)^2③Ni+3,0(u)的系数1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)^3-3/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)^2*u+3/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)*u^2-1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3