CSE 477 / G.Farin
B-spline Least Squares

Data points are data[[1]] thru data[[m+1]], associated parameter values are in tt[[1]] thru tt[[m+1]] where m>> n, n being the degree of the approximating B-spline curve.

  The B-splines are defined over a knot sequence knot[[1]] thru knot[[2n+ll+1]] where ll is the number of domain intervals.

The coefficient matrix vandermonde has m+1 rows (a lot) and   n+ll colums.
Solution: control points dd[[1]] thru dd[[n+ll]].

(SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];)

approx[n_,knot_,data_,tt_]:=
(m=Length[data]-1;
ll=Length[knot]-2*n-1; (* domain intervals *)

(* van is the Vandermonde matrix which typically has (many) more roes than columns *)
van=Table[BSplineBasis[{n,knot},j,tt[[i]]],{i,1,m+1},{j,0,n+ll-1 }];
dd=LinearSolve[Transpose[van].van,Transpose[van].data];
dataplot=ListPlot[data,PlotStyle-> Thick];
points=Graphics[{Red,PointSize[0.015],Point[dd]}];
ddplot=ListLinePlot[dd];
curve[t_]:=Sum[BSplineBasis[{n,knot},i,t]*dd[[i+1]],{i,0,n+ll-1}];
cplot=ParametricPlot[curve[t],{t,knot[[1]],knot[[2n+ll+1]]},PlotStyle->{Thick,Green}];
approxpic=Show[{ddplot,cplot, points,dataplot},Axes-> False]
);

n=3;  (* degree *)
m=70; (* counting from 0 to m *)
ll=10;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ RandomReal[], (2+t)Sin[t]+ RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]
(*Export["spiral-approx-cubic.eps",%]*)

We now use many more intervals (ll). The curve follows the data better, but it’s not a good curve anymore.

n=3;  (* degree *)
m=70; (* counting from 0 to m *)
ll=40;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

(*Increased noise level in data: *)
data1=Table[{(2+t)Cos[t]+ RandomReal[], (2+t)Sin[t]+ RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]

We now use a higher degree (n). The curve follows the data well, but it’s not a good curve either.

n=10;  (* degree *)
m=70; (* counting from 0 to m *)
ll=10;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

(*Increased noise level in data: *)
data1=Table[{(2+t)Cos[t]+ RandomReal[], (2+t)Sin[t]+ RandomReal[]},{t,0,10.0,10/m}]//N;
approx[n,knot,data1,tt1]

Next, we use only one domain interval. We get polynomial approximation (i.e., a Bezier curve),

n=7;  (* degree *)
m=70; (* counting from 0 to m *)
ll=1;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ 2RandomReal[], (2+t)Sin[t]+ 2RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]

Next, we again use only one domain interval, and the number of data points equals the number of unknowns. We get polynomial interpolation.

n=7;  (* degree *)
m=7; (* counting from 0 to m *)
ll=1;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ 2RandomReal[], (2+t)Sin[t]+ 2RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]

Next, we let the number of data equal the number of unknown control points. We get spline interpolation.

n=3;  (* degree *)
m=7; (* counting from 0 to m *)
ll=5;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ 2RandomReal[], (2+t)Sin[t]+ 2RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]

Next, we use only one interval and degree 1. We get Piecewise Linear Approximation.

n=1;  (* degree *)
m=70; (* counting from 0 to m *)
ll=10;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ 2RandomReal[], (2+t)Sin[t]+ 2RandomReal[]},{t,0,10.0,10/m}]//N;

approx[n,knot,data1,tt1]

Next, we use only one interval and degree 1. We get Linear Regression.

n=1;  (* degree *)
m=70; (* counting from 0 to m *)
ll=1;  (* domain intervals *)
tt1=Table[t,{t,0.0,1.0,1/m}]; (* data point parameter values *)

knot=Join[Table[0,{i,0,n}],Table[i/ll,{i,1,ll-1}],Table[1,{i,0,n}]]//N;

data1=Table[{(2+t)Cos[t]+ 2RandomReal[], (2+t)Sin[t]+ 2RandomReal[]},{t,0,10.0,10/m}]//N;

(* use a different order to display our geometry, simply for better viewing *)
approx1[n_,knot_,data_,tt_]:=
(m=Length[data]-1;
ll=Length[knot]-2*n-1; (* domain intervals *)

(* van is the Vandermonde matrix which typically has (many) more roes than columns *)
van=Table[BSplineBasis[{n,knot},j,tt[[i]]],{i,1,m+1},{j,0,n+ll-1 }];
dd=LinearSolve[Transpose[van].van,Transpose[van].data];
dataplot=ListPlot[data,PlotStyle-> Thick];
points=Graphics[{Red,PointSize[0.015],Point[dd]}];
ddplot=ListLinePlot[dd];
curve[t_]:=Sum[BSplineBasis[{n,knot},i,t]*dd[[i+1]],{i,0,n+ll-1}];
cplot=ParametricPlot[curve[t],{t,knot[[1]],knot[[2n+ll+1]]},PlotStyle->{Thick,Green}];
approxpic=Show[{dataplot,cplot, points,ddplot},Axes-> False]
);
approx1[n,knot,data1,tt1]