用matlab实现三次 NURBS插值曲线.docx
- 文档编号:23577126
- 上传时间:2023-05-18
- 格式:DOCX
- 页数:20
- 大小:148.97KB
用matlab实现三次 NURBS插值曲线.docx
《用matlab实现三次 NURBS插值曲线.docx》由会员分享,可在线阅读,更多相关《用matlab实现三次 NURBS插值曲线.docx(20页珍藏版)》请在冰豆网上搜索。
用matlab实现三次NURBS插值曲线
用matlab实现三次NURBS插值曲线
作者:
大漠孤狼 发表于matlab乐园()
作者:
这是我在大学时做大学生研究计划时写的,当时刚学会matlab,编写了这个程序,用了很多循环,效率不高.当时我并不清楚循环是matlab的弱点,等明白了,也不做这方面的工作了,也就懒的去改写了.如果谁需要用,就自己改吧.算法也有一些问题,我就不多说了,自己看吧
一、三次NURBS插值算法
给定平面控制顶点di(i=1,2,…n)及对应的权因子
(i=1,2,…n),可定义一条三次NURBS曲线。
先对控制顶点进行参数化,得一矢量:
U=[u0,u1,u2…,un+4]
则三次NURBS曲线的分段方程形式为:
u
[uk+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+aΔpi
其中,Δpi=pi+1-pi a=
A=
B=
在非周期情况下,首末各两数据点处的切矢可采用贝塞尔条件确定,即用如下差分矢量 Δ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
(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(d3(i+1))之间的曲线段是否有拐点出现:
如果矢量乘积ai-1×ai与ai×ai+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]上的曲线段
u∈[u3i+2,u3i+3]
令d3i=Vi,i=0,1,2,……,n
d3i+1=Vi+
(0<
i<
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)))*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+i))*u^3
1, 整体表示
为了便于编程,采用矩阵表示NURBS曲线。
三次NURBS的分段表示为
u
[uk+3,uk+4] k=0,1,2,…,n-3
设D=[ωidi]T i=0,1,2,……,n
Ω=[ω1,ω2,……,ωn] i=0,1,2,……,n
Qk=[N0,3(u),N1,3(u),……Nn,3(u)] i=0,1,2,……,n
U=[1 u u2 u3]
Qk表示u∈[uk+3,uk+4]内的基函数矢量,则NURBS的矩阵表示为
由此可见,NURBS曲线可由矩阵DΩQ表示。
进一步,可由分析矩阵而得到曲线的性质。
三、插值源程序
%给定点列V(i),随机给定权因子,构造三次NURBS插值点列V(i)
%载入数据 注意:
每一行必为矢量(三个分量),否则会出错
%参数化
V=[0.520;130;280;3100;5110;8300;10100];
s=size(V);
n=s(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;
fori=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; %求和
end
fori=3:
n
sumb=0;
forj=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
fori=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);
fori=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
fori=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)在0——1之间随机产生
%建立三维数组Q
fork=1:
(3*n-5)
fori=1:
3*n-2
ifi==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))/...
(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+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)*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+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/(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+i)-u(2+i))+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-...
u(2+i)))];
elseifi==(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))/...
(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));(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));(-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)))];
elseifi==(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;0;0;0];
end
end
end
clearaAbB;
%计算控制点
fori=2:
n
a(i,1:
3)=V(i,1:
3)-V(i-1,1:
3);
end
fori=2:
n-1
b(i,1:
3)=cross(a(i,1:
3),a(i+1,1:
3));
B(i)=b(i,3);
end
fori=2:
(n-2)
if(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 用matlab实现三次 NURBS插值曲线 matlab 实现 三次 NURBS 曲线