www.pudn.com > B-Interpolation.rar > BsplineApproximation.m, change:2006-03-28,size:2595b


%function [P] = BsplineApproximation(D,k,n,iCtrl) 
%B-spline curve approximation, 
% this is a global approximation routine 
% input: D (data point), 
%       n (n=# of data point) 
%       k (order of the curve), 
%       iCtrl (iCtrl=# control points) 
% output: P (control point of the resulting B-spline curve) 
 
D=[0 0; 3 4; -1 4; -4 0; -4 -3]; % data points 
k=2; 
n=5; 
iCtrl = 4; 
 
%Number of data points should be no less than 1 
%Number of control points should be no less than that of the data points 
%Number of control points should be no more than the Order of the curve 
 
if n<=1 
    return; 
end 
if k>n 
    k = n; 
end 
if iCtrl>n 
    iCtrl = n; 
end 
if k > iCtrl 
    k = iCtrl; 
end 
 
% Step 1: Computing parameters: ub() 
totalChord = 0; 
for i=2:n 
    totalChord = totalChord + sqrt((D(i-1,1)-D(i,1))^2+(D(i-1,2)-D(i,2))^2); 
end 
 
ub(1) = 0.; 
for i=2:n-1 
    ub(i) = ub(i-1) + sqrt((D(i-1,1)-D(i,1))^2+(D(i-1,2)-D(i,2))^2)/totalChord; 
end 
ub(n) = 1.0 ; 
 
%Step 2: Computing knot vector: U 
factor = n/iCtrl; 
 
for j=1:iCtrl-k, 
    U(j+k) = 0.; 
    for i=j:j+k-2, % averaging k-1 parameters 
        ii = round(i*factor); 
        sub = i*factor - ii; 
        U(j+k) = U(j+k) + sub*ub(round((i-1)*factor)+1) + (1-sub)*ub(ii+1); 
    end 
    U(j+k) = U(j+k)/(k-1); 
end 
 
for j=1:k, 
    U(j) = 0.; % repeating 1st k knots 
end 
 
for j=iCtrl+1:iCtrl+k, 
    U(j) = 1.; % repeating last k knots 
end 
 
% step 3: compute Qk & N 
% for(i=2:n-1) 
%     Q(i,:)=D(i,:)-N(0)*D(1,:)-N(h,k)*D(n,:); 
% end; 
 
N(1,1) = 1.0; 
N(n,iCtrl) = 1.0; 
for i=0:n-1 
    %Find span 
    if  ub(i+1)>=U(iCtrl+1), 
        span = iCtrl - 1; 
    else if ub(i+1)<=U(k), 
           span = k-1; 
        else 
            span = 0; 
            while ub(i+1)>U(span+1) 
                span = span+1; 
            end 
            span = span - 1; 
        end 
    end 
    % Define the basis function 
    funs = Basisfuns(ub(i+1),U,span,k); 
    for j=0:k-1 
        N(i+1,span-k+2+j) = funs(j+1); 
    end 
    Qk(i+1,:) = D(i+1,:) - N(i+1,1)*D(1,:) - N(i+1,iCtrl)*D(n,:); 
end 
 
% step 4: compute Q 
 
for i=0:iCtrl-1, 
    Q(i+1,1:2) = 0.; 
    for j=0:n-1, 
        Q(i+1,1:2) = Q(i+1,1:2) + N(j+1,i+1)*Qk(j+1,:);   
    end 
end 
 
% step 5: solve the function NT*N*P=Q 
%         notice that P0 & Pn-1 are not included 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
Q = Q(2:iCtrl-1,:); 
N = N(2:n-1,2:iCtrl-1); 
 
P(2:iCtrl-1,1) = (N'*N)\Q(:,1); 
P(2:iCtrl-1,2) = (N'*N)\Q(:,2); 
 
P(1,:) = D(1,:); 
P(iCtrl,:) = D(n,:);