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


function [Px,Py,Pz,U,V] = BsplineSurfApp(Dx,Dy,Dz,ku,kv,nu,nv,iCtrlu,iCtrlv) 
% B-spline surface approximation, 
% this is a global approximation routine 
% input: Dx, Dy (data point coordinates x,y correspondingly),  
%        nu,nv (n=# of data point, u,v direction correspondingly) 
%        ku,kv (order of the curve, u,v direction correspondingly)  
%        iCtrlx,iCtrly (iCtrl=# control points, u,v direction 
%        correspondingly), 
% output: Px,Py,Pz (control point of the resulting B-spline surface, 
% coordinates x,y correspondingly), 
 
% data points 
% Dx=[0 3 4 5 8 10 8 0;2 3 3 5 3 6 8 1;-1 2 3 5 3 9 7 0;-2 3 4 5 6 2 5 7];  
% Dy=[2 3 5 4 10 1 3 0;1 2 3 6 4 5 2 9;-2 3 4 5 6 2 5 7;-1 2 3 5 3 9 7 0]; 
% Dz=[2 3 5 4 10 1 3 0;1 2 3 6 4 5 2 9;-2 3 4 5 6 2 5 7;-1 2 3 5 3 9 7 0]; 
% ku=2; 
% kv=2; 
% nu=8; 
% nv=4; 
% iCtrlu = 8; 
% iCtrlv = 4; 
 
%%%%%%%stop here 
 
%Number of data points should be no less than 1 in both directions 
%Number of control points should be no less than that of the data points in 
%both directions 
%Number of control points should be no more than the Order of both 
%directions 
 
if nu<=1|nv<=1, 
    return; 
end 
if ku>nu,  
    ku = nu; 
end 
if kv>nv,  
    kv = nv; 
end 
if iCtrlu>nu, 
    iCtrlu = nu; 
end 
if iCtrlv>nv, 
    iCtrlv = nv; 
end 
if ku > iCtrlu, 
    ku = iCtrlu; 
end 
if kv > iCtrlv, 
    kv = iCtrlv; 
end 
 
% Step 1: Computing parameters: ubu(), ubv() 
ubu = zeros(1,nu); 
for l=0:nv-1, 
    totalChord = 0; 
    for k=1:nu-1, 
        cds(k) = sqrt((Dx(l+1,k+1)-Dx(l+1,k))^2+(Dy(l+1,k+1)-Dy(l+1,k))^2); 
        totalChord = totalChord + cds(k); 
    end 
     
    dt = 0;    
    for k=1:nu-1, 
        dt = dt + cds(k); 
        ubu(k+1) = ubu(k+1) + dt/totalChord; 
    end     
end 
 
ubu(1) = 0.; 
for k=1:nu-2, 
    ubu(k+1) = ubu(k+1)/nv; 
end 
ubu(nu) = 1.0 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%% Para v direction 
ubv = zeros(1,nv); 
for k=0:nu-1, 
    totalChord = 0; 
    for l=1:nv-1, 
        cds(l) = sqrt((Dx(l+1,k+1)-Dx(l,k+1))^2+(Dy(l+1,k+1)-Dy(l,k+1))^2); 
        totalChord = totalChord + cds(l); 
    end 
     
    dt = 0;    
    for l=1:nv-1, 
        dt = dt + cds(l); 
        ubv(l+1) = ubv(l+1) + dt/totalChord; 
    end     
end 
 
ubv(1) = 0.; 
for l=1:nv-2, 
    ubv(l+1) = ubv(l+1)/nu; 
end 
ubv(nv) = 1.0 ; 
 
%Step 2: Computing knot vector: U,V 
 
% 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 
 
for j=1:iCtrlu-ku, 
    U(j+ku) = 0.; 
    for i=j:j+ku-2 % averaging k-1 parameters 
        U(j+ku) = U(j+ku)+ubu(i+1); 
    end 
    U(j+ku) = U(j+ku)/(ku-1); 
end 
 
for j=1:ku, 
    U(j) = 0.; % repeating 1st k knots 
end 
 
for j=iCtrlu+1:iCtrlu+ku, 
    U(j) = 1.; % repeating last k knots 
end 
%%%%%%%%%%%%%%%%%%%%%%Calculate V 
for j=1:iCtrlv-kv, 
    V(j+kv) = 0.; 
    for i=j:j+kv-2 % averaging k-1 parameters 
        V(j+kv) = V(j+kv)+ubv(i+1); 
    end 
    V(j+kv) = V(j+kv)/(kv-1); 
end 
 
for j=1:kv, 
    V(j) = 0.; % repeating 1st k knots 
end 
 
for j=iCtrlv+1:iCtrlv+kv, 
    V(j) = 1.; % repeating last k knots 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Step 3: Get the intermediate control points Q 
 
for i=0:iCtrlv-1, 
    Qx(i+1,:) = (AppForSurf(Dx(i+1,:)',ku,nu,iCtrlu,ubu,U))'; 
    Qy(i+1,:) = (AppForSurf(Dy(i+1,:)',ku,nu,iCtrlu,ubu,U))'; 
    Qz(i+1,:) = (AppForSurf(Dz(i+1,:)',ku,nu,iCtrlu,ubu,U))'; 
end 
 
%Step 4: Get the final Control points P 
 
for j=0:iCtrlu-1, 
    Px(:,j+1) = AppForSurf(Qx(:,j+1),kv,nv,iCtrlv,ubv,V); 
    Py(:,j+1) = AppForSurf(Qy(:,j+1),kv,nv,iCtrlv,ubv,V); 
    Pz(:,j+1) = AppForSurf(Qz(:,j+1),kv,nv,iCtrlv,ubv,V); 
end 
 
%Step 5: End of program