www.pudn.com > Uniform-Rational-B-spline-fitting.rar > B样条拟合.m, change:2012-01-28,size:2238b


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