非均匀B样条拟合MATLAB程序

update on : 20.6.14

直接上代码，多的不再说了。

1、写一个Base函数

文件保存为Base.m文件

function result = Base(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))%注意1=u(i)<=t<u(i+1)=1时的情况,这里要用t<=u(i+1);
        result=1;
        return;
    else
        result=0;
        return;
    end
else
    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
end
result=alpha*Base(i,k-1,u,t)+beta*Base(i+1,k-1,u,t);

　　

2、B样条程序

文件随便保存为一个脚本文件即可，如demo.test.m。

%------------------非均匀B样条拟合MATLAB程序-----------------
clear
k=3;
% x=load('data.txt');
data_x = [-1;-1;1;1;0.2];
data_y = [1;-0.5;-1;1;0.8];
x = [data_x,data_y];



[n,m]=size(x);
%-----------弦长参数化--------------------------------------
u(k+n)=0;
for i=1:n-1
 u(k+i+1)=u(k+i)+sqrt((x(i+1,1)-x(i,1))^2+(x(i+1,2)-x(i,2))^2);
end
L=u(n+k);
for i=1:n
u(k+i)=u(k+i)/L;
end
for i=1:3
u(k+i+n)=1;
end
%控制多边线
plot(x(:,1),x(:,2),'o');
hold on
%------------反求n+2个控制点--------------------
%首位重节点v1=v2
%首位与控制多边形相切
A=zeros(n+2);
A(1,1)=1;A(1,2)=-1;
A(2,2)=1;
A(n+2,n+1)=-1;A(n+2,n+2)=1;
A(n+1,n+1)=1;
for i=3:n
  for j=0:2
    A(i,i+j-1)=Base(i+j-1,k,u,u(i+2));
   end
end
%e:方程右边.
e=0;
for i=1:m
    e(n+2,i)=0;
end
for i=1:n
    e(i+1,:)=x(i,:);
end
%求出控制点d
d=inv(A)*e;
plot(d(:,1),d(:,2),'g');


hold on
%------------插值并作出样条曲线-----------------
x=0;y=0;down=0;
for j=1:(n-1)
    uu=(u(j+3)):0.0005:u(j+4);
    for kk=1:length(uu)
       down=down+1;
        x(down)=d(j,1)*Base(j,3,u,uu(kk))+d(j+1,1)*Base(j+1,3,u,uu(kk))+d(j+2,1)*Base(j+2,3,u,uu(kk))+d(j+3,1)*Base(j+3,3,u,uu(kk));
        y(down)=d(j,2)*Base(j,3,u,uu(kk))+d(j+1,2)*Base(j+1,3,u,uu(kk))+d(j+2,2)*Base(j+2,3,u,uu(kk))+d(j+3,2)*Base(j+3,3,u,uu(kk));
    end
end
axis('equal');

plot(x,y,'red');

xlabel('x');ylabel('y');
grid on

效果

 参考文章

3次非均匀B样条拟合程序