nurbs插值曲线（可运行）.m

用MATLAB编写的求曲线拟合的程序

%给定点列V( i ),随机给定权因子，构造三次NURBS插值点列V( i )
%载入数据       注意：每一行必为矢量（三个分量），否则会出错
 
%参数化
V=[0.5 2 0;1 3 0;2 8 0;3 10 0;5 11 0;8 30 0;10 10 0];
s=size(V);
n=s(1,1);
u(1)=0;
u(2)=0;
u(3)=0;
u(4)=0;
u(3*n-1)=1;
u(3*n)=1;
u(3*n+1)=1;
u(3*n+2)=1;
 
sum=0;
for i=1:(n-1)                               %计算边矢及求和
    a(i+2,1:3)=V(i+1,1:3)-V(i,1:3);
    A=sqrt(dot(a(i+2,1:3),a(i+2,1:3)));
    sum=sum+A;                              %求和
end
for i=3:n
    sumb=0;
    for j=2:(i-1)                           %???????
        b=V(j,1:3)-V(j-1,1:3);
        B=sqrt(dot(b,b));
        sumb=sumb+B;
    end
    u(3*i-2)=sumb/sum;
end
for i=1:(n-1)
    u(3*i-1)=2/3*u(3*i-2)+1/3*u(3*i+1);
    u(3*i)=1/3*u(3*i-2)+2/3*u(3*i+1);
end
%计算切矢
a(2,1:3)=2*a(3,1:3)-a(4,1:3);
a(1,1:3)=2*a(2,1:3)-a(3,1:3);
a(n+2,1:3)=2*a(n+1,1:3)-a(n,1:3);
a(n+3,1:3)=2*a(n+2,1:3)-a(n+1,1:3);
for i=1:n+2
    b(i,1:3)=cross(a(i,1:3),a(i+1,1:3));
    A(i)=sqrt(dot(b(i,1:3),b(i,1:3)));
end
for i=1:n
    M=A(i)/(A(i)+A(i+2));
    t(i,1:3)=(1-M)*a(i+1,1:3)+M*a(i+2,1:3);
end
%随机产生权因子
w=rand(3*n-2,1);                               %w(i)在０——１之间随机产生             
%建立三维数组Q
for k=1:(3*n-5)
    for i=1:3*n-2
        if i==k
            Q(1:4,i,k)=[1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)^3;...      
  -3/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i)^2;3/(u(4+i)-u(1+i))/...
  (u(4+i)-u(2+i))/(u(4+i)-u(3+i))*u(4+i);-1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/...
  (u(4+i)-u(3+i))];
        elseif i==(k+1)
                Q(1:4,i,k)=[-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-...
   u(2+i))*u(4+i)*u(1+i)*u(3+i)-1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/...
   (u(3+i)-u(2+i))*u(4+i)^2*u(2+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/...
   (u(3+i)-u(2+i))*u(i)*u(3+i)^2;(1/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-...
   u(2+i))*u(4+i)^2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)^2+1/...
   (u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(1+i)+1/(u(4+i)-...
   u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(4+i)*u(3+i)+2/(u(3+i)-u(i))/...
   (u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(i)*u(3+i)+2/(u(4+i)-u(1+i))/(u(4+i)-...
   u(2+i))/(u(3+i)-u(2+i))*u(4+i)*u(2+i)+1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/...
   (u(3+i)-u(2+i))*u(1+i)*u(3+i));(-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-...
   u(2+i))*u(4+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(3+i)-1/...
   (u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-u(2+i))*u(2+i)-1/(u(3+i)-u(i))/...
   (u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(i)-2/(u(3+i)-u(i))/(u(3+i)-u(1+i))/...
   (u(3+i)-u(2+i))*u(3+i)-2/(u(4+i)-u(1+i))/(u(4+i)-u(2+i))/(u(3+i)-...
   u(2+i))*u(4+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))*u(1+i));...
(1/(u(4+i)-u(1+i))/(u(3+i)-u(1+i))/(u(3+i)-u(2+i))+1/(u(4+i)-u(1+i))/...
(u(4+i)-u(2+i))/(u(3+i)-u(2+i))+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(3+i)-...
u(2+i)))];
            elseif i==(k+2)
                    Q(1:4,i,k)=[1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-...
   u(1+i))*u(i)^2*u(2+i)+1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
   u(1+i))*u(4+i)*u(1+i)^2+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(i)*u(3+i)*u(1+i);(-1/(u(3+i)-u(i))/(u(2+i)-u(i))/...
   (u(2+i)-u(1+i))*u(i)^2-2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-...
   u(1+i))*u(i)*u(2+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(i)*u(1+i)-1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
   u(1+i))*u(1+i)^2-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(i)*u(3+i)-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(3+i)*u(1+i)-2/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
   u(1+i))*u(4+i)*u(1+i));(1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
   u(1+i))*u(4+i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(1+i)+2/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-...
   u(1+i))*u(i)+1/(u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-...
   u(1+i))*u(2+i)+2/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-...
   u(1+i))*u(1+i)+1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))*u(3+i));(-1/(u(3+i)-u(i))/(u(3+i)-u(1+i))/(u(2+i)-...
   u(1+i))-1/(u(4+i)-u(1+i))/(u(3+i)-u(i))/(u(2+i)-u(1+i))-1/...
   (u(3+i)-u(i))/(u(2+i)-u(i))/(u(2+i)-u(1+i)))];
                elseif i==(k+3)
                        Q(1:4,i,k)=[-u(i)^3/(u(i+3)-u(i))/(u(i+2)-u(i))/...
   (u(i+1)-u(i));3*u(i)^2/(u(i+3)-u(i))/(u(i+2)-u(i))/(u(i+1)-u(i));...
   -3*u(i)/(u(i+3)-u(i))/(u(i+2)-u(i))/(u(i+1)-u(i));1/(u(i+3)-u(i))/...
   (u(i+2)-u(i))/(u(i+1)-u(i))];
                    else
                        Q(1:4,i,k)=[0;0;0;0];
                    end
                end
            end                                                             
clear a A b B;
%计算控制点
for i=2:n
    a(i,1:3)=V(i,1:3)-V(i-1,1:3);                                          
end
for i=2:n-1
    b(i,1:3)=cross(a(i,1:3),a(i+1,1:3));
    B(i)=b(i,3);
end
for i=2:(n-2)
    if (B(i)*B(i+1)<=0)
        G(i,1:3)=(V(i,1:3)+V(i+1,1:3))./2;
    else
        e=cross(a(i+1,1:3),t(i+1,1:3));
        E=sqrt(dot(e,e));
        c=cross(t(i,1:3),t(i+1,1:3));
        C=sqrt(dot(c,c));
        G(i,1:3)=V(i,1:3)+E/C.*t(i,1:3);
    end
end
clear e E b B;
G(1,1:3)=(V(1,1:3)+V(2,1:3))./2;
G(n-1,1:3)=(V(n-1,1:3)+V(n,1:3))./2;
b=G(1,1:3)-V(1,1:3);
T(1)=sqrt(dot(b,b));
b=V(n,1:3)-G(n-1,1:3);
T(n)=sqrt(dot(b,b));
for i=2:(n-1)
    b=G(i,1:3)-V(i,1:3);
    c=V(i,1:3)-G(i-1,1:3);
    B(1)=sqrt(dot(b,b));
    B(2)=sqrt(dot(c,c));
T(i)=min(B);
end
for i=1:n
    d(3*i-2,1:3)=V(i,1:3);
end
%确定辅助点d(3i+2),d(3i+3)
%计算d(3i+2)
for i=2:n-1
    s=sqrt(dot(t(i,1:3),t(i,1:3)));
    U=[1 u(3*i) u(3*i).^2 u(3*i).^3];
    H=U*Q(1:4,3*i-1,3*i-3);
    k=w(3*i).*s.*H;
    m=1/2*w(3*i)*H.*T(i);
    d(3*i-1,1:3)=V(i,1:3)+m/k.*t(i,1:3);                
end
%计算d(3i=3)
clear i
for i=1:(n-2)
    s=sqrt(dot(t(i+1,1:3),t(i+1,1:3)));
    U=[1 u(3*i+3) u(3*i+3).^2 u(3*i+3).^3];
    H=U*Q(1:4,3*i,3*i);
    k=w(3*i+3).*s.*H;
    m=1/2*w(3*i+3)*H.*T(i+1);
    d(3*i,1:3)=V(i+1,1:3)-m/k.*t(i+1,1:3);
end
%计算矩阵D
for i=1:3*n-2
    D(i,1:2)=w(i,1).*d(i,1:2);
end
%绘图
for i=4:(3*n-2)
    j=1;
    for m=u(i):0.001:u(i+1)
        
        U=[1 m m.^2 m.^3];
        H=U*Q(1:4,1:3*n-2,i-3);
        p(i-3,1:2)=H*D./(H*w);
        x(1,j)=p(i-3,1);
        y(1,j)=p(i-3,2);
        j=j+1;
    end
    plot(x,y)
    hold on
    clear x y
end
%可选择输出矩阵W,D,Q