clc;clear;
X=load('data.txt');
n=length(X);
%控制点系数矩阵A
A=zeros(n);
A(1,1:2)=[9,-3];A(n,n-1:n)=[-3,9];A(2,1:3)=[1/4 7/12 1/6];A(n-1,n-2:n)=[1/6 7/12 1/4];
for i=3:1:n-2
A(i,i-1:i+1)=[1/6 2/3 1/6];
end
%方程右边
e=[0,0];
e(1,:)=6*X(1,:);e(n,:)=6*X(n,:);
for i=2:n-1
e(i,:)=X(i,:);
end
%得到控制点,
d=inv(A)*e;
d(n+2,:)=X(n,:);
for i=n+1:-1:2
d(i,:)=d(i-1,:);
end
d(1,:)=X(1,:);
%画出图形
hold on
%原始数据,红色,点
plot(X(:,1),X(:,2),'r.');
%控制多边形,蓝色,线
plot(d(:,1),d(:,2),'b');
%插值B样条曲线
syms u;
syms v;
if n==4
B=[-1 3 -3 1
3 -6 3 0
-3 3 0 0
1 0 0 0]';
end
if n==5
B(:,:,1)=[-1 3 -3 1
7/4 -9/2 3 0
-1 3/2 0 0
1/4 0 0 0]';
B(:,:,2)=[-1/4 3/4 -3/4 1/4
1 -3/2 0 1/2
-7/4 3/4 3/4 1/4
1 0 0 0]';
end
if n==6
B(:,:,1)=[-1 3 -3 1
7/4 -9/2 3 0
-11/12 3/2 0 0
1/6 0 0 0]';
B(:,:,2)=[-1/4 3/4 -3/4 1/4
7/12 -9/2 3 0
-7/12 1/2 1/2 1/6
1 0 0 0]';
B(:,:,3)=[-1/6 1/2 -1/2 1/6
11/12 -5/4 -1/4 7/12
-7/4 3/4 3/4 1/4
1 0 0 0]';
end
if length(d)-3>=5
B(:,:,1)=[-1 3 -3 1
7/4 -9/2 3 0
-11/12 3/2 0 0
1/6 0 0 0]';
B(:,:,2)=[-1/4 3/4 -3/4 1/4
7/12 -5/4 1/4 7/12
-1/2 1/2 1/2 1/6
1/6 0 0 0]';
for a=3:n-3
B(:,:,a)=[-1/6 1/2 -1/2 1/6
1/2 -1 0 2/3
-1/2 1/2 1/2 1/6
1/6 0 0 0]';
end
B(:,:,(n-2))=[-1/6 1/2 -1/2 1/6
1/2 -1 0 2/3
-7/12 1/2 1/2 1/6
1/4 0 0 0]';
B(:,:,n-1)=[-1/6 1/2 -1/2 1/6
11/12 -5/4 -1/4 7/12
-7/4 3/4 3/4 1/4
1 0 0 0]';
end
for i=1:(n-1)
P(i,1)=[u^3 u^2 u 1]*B(:,:,i)*d(i:i+3,1);
P(i,2)=[u^3 u^2 u 1]*B(:,:,i)*d(i:i+3,2);
uu=0:0.01:1;
for j=1:length(uu)
PX(j,i)=subs(P(i,1),'u',uu(j));PX=PX;
PY(j,i)=subs(P(i,2),'u',uu(j));PY=PY;
end
%画出图形
hold on
%拟合数据,黄色
plot(PX(:,i),PY(:,i),'y');
end
评论0