matlab实现三次插值.docx
《matlab实现三次插值.docx》由会员分享，可在线阅读，更多相关《matlab实现三次插值.docx（17页珍藏版）》请在冰豆网上搜索。
matlab实现三次插值
用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）在０——１之间随机产生
%建立三维数组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（B（
升级会员 