用matlab实现三次NURBS插值曲线.docx
- 文档编号:13298155
- 上传时间:2023-06-12
- 格式:DOCX
- 页数:20
- 大小:77.06KB
用matlab实现三次NURBS插值曲线.docx
《用matlab实现三次NURBS插值曲线.docx》由会员分享,可在线阅读,更多相关《用matlab实现三次NURBS插值曲线.docx(20页珍藏版)》请在冰点文库上搜索。
用matlab实现三次NURBS插值曲线
作者:
这是我在大学时做大学生研究计划时写的,当时刚学会matlab,编写了这个程序,用了很多循环,效率不高.当时我并不清楚循环是matlab的弱点,等明白了,也不做这方面的工作了,也就懒的去改写了.如果谁需要用,就自己改吧•算法也有一些问题,我就不多说了,自己看吧
一、三次NURBS插值算法
给定平面控制顶点di(i=l,2r-n)及对应的权因子©(i=l,2,-n)t可定义一条三次
NURBS曲线。
先对控制顶点进行参数化,得一矢量:
U=[u0,Ui,U2--.Un-l]
则三次NURBS曲线的分段方程形式为:
恥)=
Iu€[uk%uz]k=0,1,2,
…,n-3
(1)
首先证明一条性质:
若三点dk,d“d“2共线,且满足
7_©N匕3("女+3)久+级+3弘+2・3(你+3)〃奸2
线M.3(以+3)+畋+3弘+2・3他+3)
则三次NURBS曲线插值点dm
证明:
由
(1)式可得:
&(你+3)=
皱九・3(你+3)久+©+1Nr+].3("奸3)〃k+l+臥十2"知2・3(叫+3)〃女+2
©M・3“+3)+畋利N奸i・3(%3)+畋+2“奸2・3(昭3)
3如\N如\3(你.4)久+1+d+2‘V知2.3("—4)心+2+^A+3^+2.3(“A+4)久+3
©+1M+1.3(你*4)+3gNRd.、(你+4)+3gN"3(巾如4)
由以知可得
叭N匕3(你+3)久+Q却M+1.3(%3)久+1+叭+2M+2.3(%3)〃知2
3rM・3("人+3)+d+1八'如]3(心+3)+3如2"知2.3(你+3)
故性质得证。
下面构造三次NURBS插值曲线;
取一型值点列V0,VuV2,-Vn,用下列方法决定各点£处的切矢Ti:
A
a=A+B9
设顺序五个数据点Vi-2,S,•••,%,在中间数据点Vi的切矢仁被局部定义为ti=(l-a)Apji+aApi
其中,ApFpxu-pi
A」辱2X%|
在非周期情况下,首末各两数据点处的切矢可釆用贝塞尔条件确定,即用如下差分矢量△p1=2Apo—Ap,Ap-o=2ap.|
—Apo
△Pn=2Apni—Apn2APn・|=2Apn—A
Pn1
仍按前述公式确定to.tuUutno
下面,我们令d3i=pi(i=l—n)然后对数据点山进行参数化
Uo二U:
二112二山二。
工|心-如-1)
>-1
(i=2・3,…,n)
在区间[U3i•U3(iH)]插入节点U:
h41,1131*2
2
12
=§“3f+JW3(r+l)
1
~U3i+"W3(/+1)
二上》
即得节点矢量。
下面在点dai,ds(iM)(i=0,1•…,nT)之间插入两个辅助控制顶点d:
小]•chw•构成新控制顶点组,得分段三次NURBS方程。
我们按下列原则判断两型值点p.(砥),卩边(&(3)之间的曲线段是否有拐点出现:
如果矢量乘枳a»iXaj与aiXaiu同向,则Pi与Pi“之间无拐点,此时两切线Lt,LiM之间有一交点
廊+|X/;+i|f
V>Vi+l/fx/;+,l*
否则,Vi与Vi“之间有一拐点,记
Vi=(Vi+ViH>/2
令To=K-V0|
■罰H}
3,n-1)
T匕-必
我们按下列方法决定V,与Vw之间的辅助点d^2
讨论在区间血2,叽]上的曲线段
3r+2
P(“)=
上・3i-l
2>如(“)
A-3/-1
令chi=Vi•i=0.1.2•
h
(0
d3iM=Vi+吗・+1“3也3(如2)#・|
d3i.2=Vi-色+2冷+2.3("3:
+5)|如「(0
轴+2心+2.3("3,+5片,)
由此三式得
J_力31NZ3(“3i+2Mi+幺\+1弘田〕(“3i+2)〃3i+i
力31“31.3(“3i+2)+”3i+l弘田」(叫+2)
由性质可得P(lbw)=chi,即此曲线插值点ck(Vi)•
二.基函数计算及整体表示
1,基函数
根据伯恩斯坦多项式可计算
(”一曾)(%・+3-II)U-
仏+3一终)(陷3一陷1)|_陷2一山+1
%(")+
终十3一如
Ni+2.O(H)
(%+4一")("一%+|)灯一"田
(心+4一"田)(陷3一如)比+2一如
二的)+企二心的)均+3一均+2陷4一⑷+3.
整理后得:
([)NiT.O的系数
(-1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-u(1+i))-1/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(l+i))-l/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(l+i)))*u"3+(l/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))*u(4+i)+l/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(l+i))*u(i)+l/(u(3+i)-u(i))/(u(3+i)-u(l+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(l+i))*u(i)+l/(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(l+i))*u(l+i)+l/(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)-ii(i))/(u(2+i)-u(1+i))*ii(i)"2-2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i))♦u(i)*u(2+i)-l/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(l+i))*u(i)*u(l+i)-l/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(l+i))*u(l+i)*2-1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-u(1+i))⑴(3+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(1+i))*u(3+i)*u(l+i)-2/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(l+i))*u(4+i)*u(l+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)+l/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(l+i))*u(4+i)*u(1+i)"2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-u(l+i))*u(i)*u(3+i)*u(1+i)
②Ni,2.0(U)的系数
(1/(u(4+i)-u(l+i))/(ii(3+i)-u(1+i))/(u(3+i)-u(2+i))+1/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))+1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(3+i)-u(2+i)))*u*3+(-1/(u(4+i)-u(l+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)T/(u(4+i)-u(l+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)-1/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(2+i)-l/(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(l+i))/(u(3+i)-u(2+i))(3+i)-2/(u(4+i)-u(l+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(l+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))(4+i)2+1/(u(3+i)-u(i))/(u(3+i)-u(l+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(l+i)+l/(u(4+i)-u(1+i))/(u(3+i)-u(l+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(l+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(4+i)*u(2+i)+l/(u(4+i)-u(l+i))/(u(3+i)-u(l+i))/(u(3+i)-u(2+i))*u(1+i)*u(3+i))*u-l/(u(4+i)-u(l+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)-l/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(3+i)-u(2+i))*u(i)*u(3+i)「2
③Ngo(u)的系数
1/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)3-3/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)"2*u+3/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)*u*2-l/(u(4+i)-u(l+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u3
1>整体表不
为了便于编程,采用矩阵表示NURBS曲线。
三次NL1RBS的分段表示为
设D=[o)idi]
上・3
如3
工诃3仪)
u€[uu*3,urh]k=O,1,2,
i.k
i=0,1,2,
i=0,1,2,
Qk=[No.3(U)•Na.3(U),……N„.3(U)]
U=[1uu2u3]
Qk表示uG[Uk43,umi]内的基函数矢量,则NURBS的矩阵表示为
由此可见,NURBS曲线可由矩阵DQQ表示。
进一步,可由分析矩阵而得到曲线的性质。
三.插值源程序
%给定点列口i)•随机给定权因子,构造三次NURBS插值点列\,(i)
%参数化
V=[20;130;280;3100;5110:
8300;10100];s=size(V);
n=s(lt1);
u(l)=0;
u
(2)=0;
u⑶=0;
u(4)=0;
u(3*n-l)=l;
u(3*n)=l;
u(3*n+l)=l;
u(3*n+2)=l;sum=0;
fori=l:
(n-1)%计算
边矢及求和
a(i+2fl:
3)=V(i+l,l:
3)-V(ifl:
3);
A=sqrt(dot(a(i+2,1:
3),a(i+2,1:
3)));
sum=sum+A;%求
和
end
fori=3:
n
sumb=0;forj=2:
(i-1)
b=V(jtl:
3)-V(j-lJ:
3);
B二sqrt(dot(b,b));
sumb=sumb^B;
end
u(3*i-2)=sumb/sum;
end
fori=l:
(n-1)
u(3*i-l)=2/3*u(3*i-2)+l/3*u(3*i+l):
u(3*i)=l/3*u(3*1-2)+2/3*u(3*i+l);
end
%计算切矢
a(2,l:
3)=2*a(3tl:
3)-a(4,l:
3);
a(l,l:
3)=2*a(2tl:
3)-a(3,l:
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+1t1:
3);
fori=l:
n+2
b(i,1:
3)=cross(a(i,1:
3)ta(i+l,1:
3));
A(i)=sqrt(dot(b(i,1:
3),b(i,1:
3)));
end
fori=l:
n
M=A(i)/(A(i)+A(i+2));
t(i,l:
3)=(l-M)*a(i+l,l:
3)+M*a(i+2tl:
3);
end
%随机产生权因子
%w
w=rand(3*n-2t1);
(i)在0——1之间随机产生
%建立三维数组Q
fork=l:
(3*n-5)
fori=l:
3*n-2
ifi==k
Q(l:
4,itk)=[l/(u(4+i)-u(l+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(l+i))/...
(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))];
elseifi==(k+l)
Q(l:
4,i,k)=[-l/(u(4+i)-u(1+i))/(u(3+i)-u(l+i))/(
u(3+i)-・..
u(2+i))*u(4+i)*u(l+i)*u(3+i)-l/(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)-l/(u(3+i)-u(i))/(u(3+i)-u(l+i))/...
(u(3+i)-u(2+i))*u(i)*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(l+i))/(u(3+i)-u(2+i))*u(3+i)*2+
1/…
(u(4+i)-u(l+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(1+i)+1/(u(4+i)-・•
u(l+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(l+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)和(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));(T/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)
u(2+i))*u(4+i)T/(u(4+i)-u(1+i))/(u(3+i)-u(l+i))/(u(3+i)-u(2+i))*u(3+i)T/
•••
(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))*11⑴-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(l+i))/(u(4+i)-u(2+i))/(u(3+i)-...u(2+i))*u(4+i)T/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(l+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(l+i))/・・.
(u(4+i)-u(2+i))/(u(3+i)-u(2+i))+1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(3+i)-..・u(2+i)))];
elseifi==(k+2)
Q(l:
4,i,k)=[l/(u(3+i)-u(i))/(u(2+i)-u(i)
)/(u(2+i)-..・
u(l+i))*u(i)2*u(2+i)+l/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-..•u(l+i))*u(4+i)*u(l+i)2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i).
u(l+i))*u(i)*u(3+i)*u(1+i);(-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(l+i))*u(i)*u(2+i)T/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-..・u(l+i))*u(i)*u(l+i)-l/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-..・u(l+i))*u(1+i)2-1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-.・.
u(l+i))*u(i)*u(3+i)T/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-..・u(l+i))*u(3+i)*u(1+i)-2/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
u(l+i))*u(4+i)*u(1+i));(1/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-...
u(l+i))*u(4+i)+l/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-…・
u(l+i))*u(i)+l/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-..・
u(l+i))*u(l+i)+2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-.•.
u(l+i))*u(i)+l/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-..•
u仃+i))*u(2+i)+2/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-..・
u(l+i))*u(l+i)+l/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-...
u(l+i))*u(3+i)):
(-1/(u(3+i)-u(i))/(u(3+i)-u(l+i))/(u(2+i)-….
u(l+i))-l/(u(4+i)-u(l+i))/(u(3+i)-u(i))/(u(2+i)-u(l+i))-1/..・
(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(l+i)))];
elseifi==(k+3)
Q(l:
4,iFk)=[-u(i)*3/(u(i+3)-u(i))/(u(i+2)-u(i))/...
(u(i+l)-u(i));3*u(i)*2/(u(i+3)-u(i))/(u(i+2)-u(i))/(u(i+l)-u⑴);...-3*u(i)/(u(i+3)-u(i))/(u(i+2)-u(i))/(u(i+l)-u(i));1/(u(i+3)-u(i))/...(u(i⑵-u(i))/(u(i+l)-u(i))];
else
Q(l:
4,i,k)=[0;0;0;0];
end
end
end
clearaAbB;
%计算控制点
fori=2:
n
a(itl:
3)=V(i,l:
3)-V(i-l,l:
3);
endfori=2:
nT
b(i,l:
3)=cross(a(i,1:
3)fa(i+l,1:
3));
B(i)=b(i,3);
end
fori=2:
(n-2)
if(B(i)*B(i+l)<=0)
G(i,l:
3)=(V(iJ:
3)+V(i+1J:
3))./2;
else
e=cross(a(i+lt1:
3),t(i+lt1:
3));
E=sqrt(dot(e,e));
c=cross(t(i,1:
3),t(i+l,1:
3));
C=sqrt(dot(c,c));
G(i,l:
3)=V(i,l:
3)+E/C.*t(i,l:
3);
end
end
cleareEbB;
G(l,l:
3)=(V(lJ:
3)+V(2,l:
3))./2;
G(n-lJ:
3)=(V(n-l,l:
3)+V(n,l:
3))./2;
b=G(lJ:
3)-V(lJ:
3);
T(l)=sqrt(dot(b,b));
b=V(nt1:
3)-G(n~l,1:
3);
T(n)=sqrt(dot(b,b));
fori=2:
(n-1)
b=G(i,l:
3)-V(i,l:
3);
c=V(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 实现 三次 NURBS 曲线