Bezier  bicubic surface interpolation

 We are given an array pp  of size 4x4, containing data points to be interpolated. Data points have uniform parameter values.

knotu={0,0.333,0.6667,1};
knotv={0,0.333,0.6667,1};
(* Vandermonde matrices vandu, vandv: *)
vandu=Transpose[Table[BernsteinBasis[3,i,knotu[[j]]],{i,0,3},{j,1,4}]];
vandv=Table[BernsteinBasis[3,j,knotv[[i]]],{j,0,3},{i,1,4}];

(*  array ppp of data sets: *)
ppp={
Table[{j*Cos[i],j*Sin[i],.5Sqrt[(i^2+j^2)]},{i,1,4},{j,1,4}]//N,
Table[{i,j,i^2/5-i+Sin[j]},{i,1,4},{j,1,4}]//N};

pp=ppp[[2]];

(* Resulting contol net:*) 
bb=Inverse[vandu].Transpose[Inverse[Transpose[vandv]].Transpose[pp]];
(* We need this because Mathematica's way of matrix multiplication if
one matrix holds column vectors. It is equavialent to 
Inverse[vandu].pp.Inverse[vandv] *)
center=Sum[Sum[bb[[i,j]],{i,1,4}],{j,1,4}]/16;
dists=Table[Norm[bb[[i,j]]-center],{i,1,4},{j,1,4}];
radius=Max[dists];
Do[bb[[i,j]]=bb[[i,j]]-center,{i,1,4},{j,1,4}];
Do[pp[[i,j]]=pp[[i,j]]-center,{i,1,4},{j,1,4}];

netplot=Graphics3D[{PointSize[Large],Map[Point,pp],Gray,Line[bb],Line[Transpose[bb]]}];
bbplot=Graphics3D[{PointSize[Large],Red,Map[Point,bb]}];
surfplot=Graphics3D[BSplineSurface[bb]];

Show[{bbplot,surfplot,netplot},Boxed->False,PlotRange->0.8radius]

We are given an array pp  of size(nu+1)(nv+1), containing data points to be interpolated. Data points have uniform parameter values.

nu = 13;
nv = 5;
knotu=Table[i/(nu),{i,0,nu}]//N;
knotv=Table[j/(nv),{j,0,nv}]//N;

(* Vandermonde matrices vandu, vandv: *)
vandu=Transpose[Table[BernsteinBasis[nu,i,knotu[[j]]],{i,0,nu},{j,1,nu+1}]];
vandv=Table[BernsteinBasis[nv,j,knotv[[i]]],{j,0,nv},{i,1,nv+1}];

(*  array ppp of data sets: *)
ppp={
Table[{Sqrt[j]*Cos[i],Sqrt[j]*Sin[i],.5Sqrt[(i^2+j^2)]},{i,1,nu+1},{j,1,nv+1}]//N,
Table[{i,j,i^2/5-i+Sin[j]},{i,1,nu+1},{j,1,nv+1}]//N};

pp=ppp[[2]];
(* Resulting contol net:*) 
bb=Inverse[vandu].Transpose[Inverse[Transpose[vandv]].Transpose[pp]];
(* We need this because Mathematica's way of matrix multiplication if
one matrix holds column vectors. It is equavialent to 
Inverse[vandu].pp.Inverse[vandv] *)
center=Sum[Sum[pp[[i,j]],{i,1,nu+1}],{j,1,nv+1}]/((nu+1)*(nv+1));
Do[bb[[i,j]]=bb[[i,j]]-center,{i,1,nu+1},{j,1,nv+1}];
dists=Table[Norm[pp[[i,j]]-center],{i,1,nu+1},{j,1,nv+1}];
radius=Max[dists];
Do[pp[[i,j]]=pp[[i,j]]-center,{i,1,nu+1},{j,1,nv+1}];

netplot=Graphics3D[{Gray,Line[bb],Line[Transpose[bb]]}];
dataplot=Graphics3D[{PointSize[Large],Map[Point,pp]}];
bbplot=Graphics3D[{PointSize[Large],Red,Map[Point,bb]}];
surfplot=Graphics3D[{BSplineSurface[bb,SplineDegree->{nu,nv}]}];

Show[{surfplot,dataplot},Boxed->False,PlotRange->0.9radius]