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

1、用matlab实现三次 NURBS插值曲线作者:大漠孤狼 发表于matlab乐园()作者:这是我在大学时做大学生研究计划时写的,当时刚学会matlab ,编写了这个程序,用了很多循环,效率不高.当时我并不清楚循环是matlab的弱点,等明白了,也不做这方面的工作了,也就懒的去改写了.如果谁需要用,就自己改吧.算法也有一些问题,我就不多说了,自己看吧 一 、三次NURBS插值算法 给定平面控制顶点di(i=1,2,n)及对应的权因子(i=1,2,n)，可定义一条三次NURBS曲线。先对控制顶点进行参数化，得一矢量： U=u0,u1,u2,un+4 则三次NURBS曲线的分段方程形

2、式为： uuk+3,uk+4 k=0,1,2,n-3 (1) 首先证明一条性质： 若三点dk,dk+1,dk+2共线，且满足 则三次NURBS曲线插值点dk+1. 证明：由（1）式可得： 由以知可得 故性质得证。 下面构造三次 NURBS插值曲线； 取一型值点列V0,V1,V2,Vn,用下列方法决定各点Vi处的切矢Ti: 设顺序五个数据点Vi-2,Vi-1,Vi+2,在中间数据点Vi的切矢ti被局部定义为 ti=(1-a)pi-1+api 其中，pi=pi+1-pi a=, A= B= 在非周期情况下，首末各两数据点处的切矢可采用贝塞尔条件确定，即用如下差分矢量 p-1=2p0p1 p-2=2

3、p-1p0 pn=2pn1pn2 pn+1=2pnpn1 仍按前述公式确定t0,t1,tn-1,tn。 下面，我们令d3i=pi (i=1n) 然后对数据点d3i进行参数化 u0=u1=u2=u3=0 (i=2,3,，n) u3n+1=u3n+2=u3n+3=u3n+4=1 在区间u3i,u3(i+1)插入节点u3i+1 ,u3i+2 u3i+1= 即得节点矢量。 下面在点d3i ,d3(i+1) (i=0 , 1 , ,n-1) 之间插入两个辅助控制顶点 d3i+1 , d3i+2 ,构成新控制顶点组，得分段三次NURBS方程。 我们按下列原则判断两型值点 pi (d3i) , pi+1(

4、d3(i+1) ) 之间的曲线段是否有拐点出现：如果矢量乘积ai-1ai 与aiai+1同向，则pi 与pi+1之间无拐点，此时两切线Li , Li+1之间有一交点 Vi*=Vi+ 否则，Vi 与Vi+1 之间有一拐点，记 Vi*=（Vi+Vi+1）2 令 T0 = Ti = min , (i = 1 , 2 , 3 , n-1) Tn = 我们按下列方法决定 Vi 与 Vi+1 之间的辅助点 d3i+1 , d3i+2 讨论在区间 u3i+2 , u3i+3上的曲线段 uu3i+2 , u3i+3 令 d3i = Vi , i = 0 , 1, 2 , n d3i+1 = Vi + (0 i

5、 Ti ) d3i+2 = Vi - (0 i+1 Ti+1 ) 由此三式得 由性质可得 P( u3i+2 ) = d3i , 即此曲线插值点 d3i ( Vi ) .二、基函数计算及整体表示1， 基函数根据伯恩斯坦多项式可计算Ni , 3( u )=+ + + 整理后得：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)*u3+(1/(u(4+i)-u(1+

6、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)-

7、u(i)/(u(3+i)-u(1+i)/(u(2+i)-u(1+i)*u(3+i)*u2+(-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

8、(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

9、(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)*u3+(-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)/(

10、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)*u2+(1/(u(4+i)-u(1+i)/(u(4+i)-u(2+i)/(u(3+i)-u(2+i)*u(

11、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

12、/(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)2Ni+3,0(u) 的系数1/(u(4+i)-u(1+i)/(u(4+i)-u(2+i)/(u(4+i

13、)-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)*u2-1/(u(4+i)-u(1+i)/(u(4+i)-u(2+i)/(u(4+i)-u(3+i)*u31， 整体表示为了便于编程，采用矩阵表示NURBS曲线。三次NURBS的分段表示为 uuk+3,uk+4 k=0,1,2,n-3 设D=idi T i= 0 , 1 , 2 , , n = 1 , 2 , , n i= 0 , 1 , 2

14、, , n Qk =N0 , 3(u) , N1 , 3(u) , Nn , 3(u) i= 0 , 1 , 2 , , n U =1 u u2 u3 Qk 表示uuk+3 ,uk+4 内的基函数矢量,则NURBS的矩阵表示为 由此可见，NURBS曲线可由矩阵 D Q 表示。进一步，可由分析矩阵而得到曲线的性质。三、插值源程序%给定点列V( i ),随机给定权因子，构造三次NURBS插值点列V( i )%载入数据 注意：每一行必为矢量（三个分量），否则会出错%参数化V=0.5 2 0;1 3 0;2 8 0;3 10 0;5 11 0;8 30 0;10 10 0;s=size(V);n=s(

15、1,1);u(1)=0;u(2)=0;u(3)=0;u(4)=0;u(3*n-1)=1;u(3*n)=1;u(3*n+1)=1;u(3*n+2)=1;sum=0;for i=1:(n-1) %计算边矢及求和 a(i+2,1:3)=V(i+1,1:3)-V(i,1:3); A=sqrt(dot(a(i+2,1:3),a(i+2,1:3); sum=sum+A; %求和endfor i=3:n sumb=0; for j=2:(i-1) %? b=V(j,1:3)-V(j-1,1:3); B=sqrt(dot(b,b); sumb=sumb+B; end u(3*i-2)=sumb/sum;end

16、for i=1:(n-1) u(3*i-1)=2/3*u(3*i-2)+1/3*u(3*i+1); u(3*i)=1/3*u(3*i-2)+2/3*u(3*i+1);end%计算切矢a(2,1:3)=2*a(3,1:3)-a(4,1:3);a(1,1:3)=2*a(2,1:3)-a(3,1:3);a(n+2,1:3)=2*a(n+1,1:3)-a(n,1:3);a(n+3,1:3)=2*a(n+2,1:3)-a(n+1,1:3);for i=1:n+2 b(i,1:3)=cross(a(i,1:3),a(i+1,1:3); A(i)=sqrt(dot(b(i,1:3),b(i,1:3);end

17、for i=1:n M=A(i)/(A(i)+A(i+2); t(i,1:3)=(1-M)*a(i+1,1:3)+M*a(i+2,1:3);end%随机产生权因子w=rand(3*n-2,1); %w(i)在之间随机产生 %建立三维数组Qfor k=1:(3*n-5) for i=1:3*n-2 if i=k Q(1:4,i,k)=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;3/(u(4+i)-u(1+i)/. (

18、u(4+i)-u(2+i)/(u(4+i)-u(3+i)*u(4+i);-1/(u(4+i)-u(1+i)/(u(4+i)-u(2+i)/. (u(4+i)-u(3+i); elseif i=(k+1) Q(1:4,i,k)=-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)*

19、u(3+i)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

20、+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);(-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/

21、(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);.(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

22、+i)-u(2+i)+1/(u(3+i)-u(i)/(u(3+i)-u(1+i)/(u(3+i)-.u(2+i); elseif i=(k+2) Q(1:4,i,k)=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);(-1/(u(3+i)-u(i)/(u(2+i)-u(i)/

23、. (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)*

24、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);(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)

25、/(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);(-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); elseif i=(k+3) Q(1:4,i,k)=-u(i)3/(u(i+3)-u(i)/(u(i+2)-u(i)/. (u(i+1)-u(i);3*u(i)2/(u(i+3)-u(i)/(u(i+2)-u(i)/(u(i+1)-u(i);. -3*u(i)/(u(i+3)-u(i)/(u(i+2)-u(i)/(u(i+1)-u(i);1/(u(i+3)-u(i)/. (u(i+2)-u(i)/(u(i+1)-u(i); else Q(1:4,i,k)=0;