- 积分
- 0
- 注册时间
- 2011-7-1
- 仿真币
-
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2012-6-9 20:54:47
|
显示全部楼层
来自 天津
本帖最后由 5爱football 于 2012-6-9 20:55 编辑
此程序是曲线的,好像有点错误,请大牛们指教,多谢!最好给个曲面重构的
%X:原始资料,d:控制顶点
%n:数据条数,k:B样条的次数
%
%
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+1)=u(k+i)+abs((X(i+1,1)-X(i,1))^2+(X(i+1,2)-X(i,2))^2+(X(i+1,3)-X(i,3))^2)^(1/2);
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;
return;
else
result=0;
return;
end
end
if u(i+k)-u(i)==0
alpha=0;
else
alpha=(t-u(i))/(u(i+k)-u(i));
end
if u(i+k+1)-u(i+1)==0
beta=0;
else
beta=(u(i+k+1)-t)/(u(i+k+1)-u(i+1));
end
result=alpha*Bbase(i,k-1,u,t)+beta*Bbase(i+1,k-1,u,t);
end
%************递归,德布尔递推算法***************************************
function y=DeBoor(d,u,kn,i,L)
%控制点序列,u节点序列,kn插入的节点,u(i)<kn<u(i+1),L递推层次
% k B样条次数
k=3;
if L==0
y=d(i,:);
return;
end
temp=(kn-u(i))/(u(i+k+1-L)-u(i));
y=(1-temp)*DeBoor(d,u,kn,i-1,L-1)+temp*DeBoor(d,u,kn,i,L-1);
|
|