查看: 4993|回复: 8

非均匀三次B样条插值和插入节点算法

字体大小: 正常放大

2620 主题	162 听众	1万积分

升级 0%

TA的每日心情

	开心 2015-3-12 15:35

签到天数: 207 天

[LV.7]常住居民III

群组: 第六届国赛赛前冲刺培

群组: 国赛讨论

群组: 2014美赛讨论

群组: 2014研究生数学建模竞

群组: 数学中国试看培训视频

电梯直达

1^#

发表于 2015-1-15 10:25 |只看该作者 |倒序浏览

|招呼Ta 关注Ta

%X:原始资料,d:控制顶点
%n:数据条数,k:B样条的次数
%
%see also http://www.matlabsky.com
%
X=load('data.txt');
[n,numy]=size(X); %得数据维数;
k=3;
%弦长参数化
u(k+1)=0;
for i=1

n-1)
u(k+i+1)=u(k+i)+abs(X(i+1,1)-X(i,1));
end
%====规范参数化=========================
temp=u(n+k);
for i=(k)

n+k) %(k+1)

n-1+k)足够了
u(i)=u(i)/temp;
end
%首末节点重复k+1个
for i=1

k+1)
u(i)=0;
u(n-1+k+i)=1;
end
%% ===============================
%A:方程系数------采用递推算法-----------------------
A=zeros(n+2);
%控制多形与曲线相切,切曲线图端点为第1个控制点重合
%A(1,1)=1;A(1,2)=-2;A(1,3)=1;A(n+2,n)=1;A(n+2,n+1)=-2;A(n+2,n+2)=1;
%控制多形与曲线相切,切曲线图端点和第1,2个控制点重合
A(1,1)=1;A(1,2)=-1;A(n+2,n+1)=-1;A(n+2,n+2)=1;
A(2,2)=1;A(n+1,n+1)=1;
for i=3:n
for j=0:2
A(i,i+j-1)=Bbase(i+j-1,3,u,u(2+i));
end
end
%e:方程右边.
e=0;
for i=1:numy
e(n+2,i)=0;
end
for i=1:n
e(i+1,

=X(i,

;
end
%得到控制点,
d=inv(A)*e;
%clear A;
%clear e;
%===画出图形================================
hold on
%原始数据,红色,点
plot(X(:,1),X(:,2),'r.');
%控制多边形,蓝色,线
plot(d(:,1),d(:,2),'b');
%===插值B样条曲线============================
x=0;y=0;down=0;
for j=1

n-1)
uu=(u(j+3)):0.005:u(j+4);
for kk=1:length(uu)
      temp=j+4;
      down=down+1;
      x(down)=d(j,1)*Bbase(temp-4,3,u,uu(kk))+d(j+1,1)*Bbase(temp-3,3,u,uu(kk))+d(j+2,1)*Bbase(temp-2,3,u,uu(kk))+d(j+3,1)*Bbase(temp-1,3,u,uu(kk));
      y(down)=d(j,2)*Bbase(temp-4,3,u,uu(kk))+d(j+1,2)*Bbase(temp-3,3,u,uu(kk))+d(j+2,2)*Bbase(temp-2,3,u,uu(kk))+d(j+3,2)*Bbase(temp-1,3,u,uu(kk));
end
end
plot(x,y,'g');
%clear uu x y;
%===输入增加的节点并得到该点数据==============
fprintf('请输入一个%d-%d之间的',min(X(:,1)),max(X(:,1)) );
Knot=input('一个补充节点:');
%Knot=25;
if Knot<min(X(:,1)) || Knot>max(X(:,1))
fprintf('你的输入有误,超过了范围.\n');
return;
end
%转换为节点
knot=(Knot-min(X(:,1)))/(max(X(:,1))-min(X(:,1)));
kn=find(u>knot,1,'first');
newu=u;
newu(kn)=knot;
newu(kn+1:length(u)+1)=u(kn:length(u));
%===得到插入节点的数据======================
inknot=DeBoor(d,u,knot,kn-1,3);
plot(inknot(1),inknot(2),'m.');

%==生成新的控制序列dd,Deboor算法=============
dd=d;
dd(n+3,

=0;
dd(kn:n+3,

=dd(kn-1:n+2,

;
for i=kn-k

kn-1)
dd(i,

=DeBoor(d,u,knot,i,1);
end
%新控制多边形,蓝色,线
plot(dd(:,1),dd(:,2),'b--');
%===画出修正后的图=========================
x=0;y=0;down=0;
for j=1:n
uu=(newu(j+3)):0.005:newu(j+4);
for kk=1:length(uu)
      temp=j+4;
      down=down+1;
      x(down)=dd(j,1)*Bbase(temp-4,3,newu,uu(kk))+dd(j+1,1)*Bbase(temp-3,3,newu,uu(kk))+dd(j+2,1)*Bbase(temp-2,3,newu,uu(kk))+dd(j+3,1)*Bbase(temp-1,3,newu,uu(kk));
      y(down)=dd(j,2)*Bbase(temp-4,3,newu,uu(kk))+dd(j+1,2)*Bbase(temp-3,3,newu,uu(kk))+dd(j+2,2)*Bbase(temp-2,3,newu,uu(kk))+dd(j+3,2)*Bbase(temp-1,3,newu,uu(kk));
end
end
plot(x,y,'r-.');
hold off

%****************第i段k次B样条基,Deboor递推递归算法******************************
function result = Bbase(i,k,u,t)
%第i段k次B样条基,Deboor递推递归算法
%t为变量,u(i)<=t<u(i+1),k=0时result=1;
if k==0
if (u(i)<=t && t<u(i+1)) %注意t=u(i+1)) 时的情况,要用t<=u(i+1);
      result=1;

游客，如果您要查看本帖隐藏内容请回复

zan