www.pudn.com > B-Interpolation.rar > BsplineSurfInter.m, change:2006-04-10,size:3041b

```function [Px,Py,Pz,U,V] = BsplineSurfInter(Dx,Dy,Dz,ku,kv,nu,nv)
%B-spline surface interpolation,
% this is a global interpolation 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),
% output: Px,Py,Pz (control point of the resulting B-spline surface, coordinates x,y,z correspondingly)

%test data;
% 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]; % data points
% 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];
% ku=2;
% kv=2;
% nu=8;
% nv=4;

% Step 0:
% Number of data points should be no less than the Order of the surface
if(ku>nu)
ku = nu;
end

if(kv>nv)
kv = nv;
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+(Dz(l+1,k+1)-Dz(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+(Dz(l+1,k+1)-Dz(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:nu-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=nu+1:nu+ku,
U(j) = 1.; % repeating last k knots
end
%%%%%%%%%%%%%%%%%%%%%%Calculate V
for j=1:nv-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=nv+1:nv+kv,
V(j) = 1.; % repeating last k knots
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Step 3: Get the intermediate control points Q

for i=0:nv-1,
Qx(i+1,:) = (InterForSurf(Dx(i+1,:)',ku,nu,ubu,U))';
Qy(i+1,:) = (InterForSurf(Dy(i+1,:)',ku,nu,ubu,U))';
Qz(i+1,:) = (InterForSurf(Dz(i+1,:)',ku,nu,ubu,U))';
end

%Step 4: Get the final Control points P

for j=0:nu-1,
Px(:,j+1) = InterForSurf(Qx(:,j+1),kv,nv,ubv,V);
Py(:,j+1) = InterForSurf(Qy(:,j+1),kv,nv,ubv,V);
Pz(:,j+1) = InterForSurf(Qz(:,j+1),kv,nv,ubv,V);
end

%Step 5: End of program

```