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,:);

```