CSE 477 / G. Farin
Bezier and B - spline surface least squares approximation

 We are given an array data  of size numpts, containing data points to be approximated. The resulting surface has linearized control net ddlinear, which becomes the array ddapprox after packing to an (nu+lu)(nv+lv) rectangular control net net ddapprox. All data parameter values params are in the unit square. When there are none near the corners of the unit square, the process becomes unstable.

nu=8;  (*u-degree *)
nv=3; (*v-degree*)
llu=1; (* number of u-domain intervals *)
llv=1; (* number of v-domain intervals *)
(* for llu = llv = 1, we get Bezier surfaces *)
cols=nv+llv;
rows=nu+llu;
numpts=rows*cols*4; (* number of (linearized) data points *)
noise = 0.1; (* noise in data points *)
knotu = Join[Table[0,{i,0,nu}],Table[i/llu,{i,1,llu-1}],Table[1,{i,0,nu}]]//N;
knotv = Join[Table[0,{i,0,nv}],Table[i/llv,{i,1,llv-1}],Table[1,{i,0,nv}]]//N;

linearize[cols_,i_,j_]:=(cols)*(i-1) + j; (* produces a linear count in a rectangular array with cols columns *)
index[cols_,count_]:=(int=QuotientRemainder[count,cols]+{1,0};
If[int[[2]]==0,{int[[1]]-1,cols},int]);
(* Undoes linearize[] *)

ddd={Table[{i,j,i^2/5-i+Sin[j]},{i,nu+llu},{j,nv+llv}]//N,
Table[{j*Cos[i],j*Sin[i],.5Sqrt[(i^2+j^2)]},{i,1,nu+llu},{j,1,nv+llv}]//N};
dd=ddd[[2]];

surf[u_,v_]:=Sum[Sum[BSplineBasis[{nu,knotu},i,u]*BSplineBasis[{nv,knotv},j,v]*dd[[i+1,j+1]],{i,0,nu+llu-1}],{j,0,nv+llv-1}]; (* surface from which data points will be read off from *)

params=Join[Table[{RandomReal[],RandomReal[]},{i,1,numpts}],{{0,0},{1,0},{0,1},{1,1}}];
numpts = numpts+4; 
(*params=Table[{0.005(1+i)*Sin[2i]+0.5,0.005(1+i)*Cos[2i]+0.5},{i,10,numpts+10}];  *)
(* the range for random parameter values. These need to cover all of surf's domain.
 That's why the four corner parameters are added. *)

data=Table[surf[params[[i,1]],params[[i,2]]]+{0,0,noise*RandomReal[]},{i,1,numpts}]; (* data points, read off from an exising B-spline surface surf *)

prod[k_,u_,v_]:=BSplineBasis[{nu,knotu},index[cols,k][[1]]-1,u]*BSplineBasis[{nv,knotv},index[cols,k][[2]]-1,v]; (* product of two B-splines which have linearized index k *)

vand = Table[prod[k,params[[i,1]],params[[i,2]]],{i,1,numpts},{k,1,rows*cols}] ;
(* the nonsquare Vandermonde matrix *)
ddlinear=Inverse[Transpose[vand].vand ].Transpose[vand].data;
(* the initial, linearized, solution *)

ddapprox=Table[{0,0,0},{i,nu+llu},{j,nv+llv}];
Do[ddapprox[[index[cols,k][[1]],index[cols,k][[2]] ]] =ddlinear[[k]],{k,1,cols*rows}];
(* The initial solution packed into a rectangular array *)

appnetplot=Graphics3D[{PointSize[Large],Green,Map[Point,ddapprox],Gray,Line[ddapprox],Line[Transpose[ddapprox]]}];
dataplot=Graphics3D[{PointSize[Large],Map[Point,data]}];
appsurfplot=Graphics3D[BSplineSurface[ddapprox,SplineKnots->{knotu,knotv},SplineDegree-> {nu,nv}]];

Show[{dataplot(*appnetplot*),appsurfplot},Boxed->False,ViewPoint->{-2,0,1}]