B-Spline basis functions

The Main function

(*B-SPline function of arbitrary degree*)
(*The func[x_]:=func[x]=.. structure turns the function into a lookup function (it remembers the values once there are evaluated)*)
B[n_][x_]:=B[n][x]=B[n][-x]=Simplify[(1/n!)*Sum[Binomial[n+1,k]*(x+(n+1)/2-k)^n*(-1)^k*UnitStep[x+(n+1)/2-k],{k,0,n+1}]]

Interval function (on which the B-spline at a certain scale and position is non-zero)

(*The interval on which the B-spline is non-zero*)
BXMin[n_,xcenter_,sx_]:=-sx((n+1)/2)+xcenter
BXMax[n_,xcenter_,sx_]:=+sx((n+1)/2)+xcenter
BNonZeroInterval[n_,xcenter_,sx_]:={BXMin[n,xcenter,sx],BXMax[n,xcenter,sx]}

Sample a N-dimensional single B-spline on a spatial grid

(*Sample N-D spline function*)
(*centers: a list of length N specifying the location/index around which the spline is centered*)
(*scales: the scaling in each dimension*)
(*size: the size of the output array*)
BSplineSampled[n_,centers_,scales_,size_,sparseQ_:False]:=Module[
{
dim=Length[size],(*The N in N-D*)
index,
bases1D,
basisND
},
bases1D=Table[
(*Each 1D basis is highly sparse with only values around its center*)
SparseArray[{{index_}/;IntervalMemberQ[Interval[BNonZeroInterval[n,centers[[i]],scales[[i]]]],index]->N[B[n][((index-centers[[i]])/(scales[[i]]))]]},{size[[i]]}]
,{i,1,dim,1}];
(*This is to make sure that the code also works for 1D signals (s.t. the Kronecker product can be used)*)
bases1D=Prepend[bases1D,{1}];
(*The final basis function is the Kronecker product of the 1D basis functions*)
basisND=ArrayReshape[KroneckerProduct@@bases1D,size];
(*Return*)
If[sparseQ,
Return[basisND]
,
Return[Normal[basisND]]
]
]

Compute the spatial and center of the B-Spline corresponding to a certain index

BSplineIndexToCenter[ijk_,basisSize_,dataSize_]:=Module[
{
scales,
shifts
},
scales=dataSize/basisSize;
shifts=(scales-1)/2;
Return[ijk*scales-shifts]
]

Generate the complete basis, the full set of sampled B-spline functions

BSplineBasisUniform[n_,basisSize_,dataSize_]:=Module[{center,scales=dataSize/basisSize,ijk},
SparseArray@ParallelArray[BSplineSampled[n,BSplineIndexToCenter[{##},basisSize,dataSize],scales,dataSize,True]&,basisSize]
]
BSplineBasisUniformFlat[n_,basisSize_,dataSize_]:=ArrayReshape[Flatten[BSplineBasisUniform[n,basisSize,dataSize]],{Times@@basisSize,Times@@dataSize}]

BSplineBasisUniform2[n_,basisSize_,dataSize_]:=Module[{center,scales=dataSize/basisSize,shifts,ijk},
shifts=(scales-1)/2;
SparseArray@ParallelArray[BSplineSampled[n,{##}*scales-shifts,scales,dataSize,True]&,basisSize]
]
BSplineBasisUniformFlat2[n_,basisSize_,dataSize_]:=ArrayReshape[Flatten[BSplineBasisUniform2[n,basisSize,dataSize]],{Times@@basisSize,Times@@dataSize}]

Example 1D

Fit spline basis of scale 1 (same nr of basis functions as there are "pixels")

The random signal

signal=RandomReal[{-1,1},40];
ListLinePlot[signal,InterpolationOrder->0,PlotStyle->Blue]

Settings and construct basis (the code for N-D basis construction is not yet optimized, this may take a while for large data)

n=2
dataSize=Dimensions[signal](*Number of data elements*)
basisSize=Dimensions[signal](*Number of basis functions to use*)
scale=First[dataSize/basisSize](*Scale the basisfunctions to make sure the whole domain is covered*)
bFlat=BSplineBasisUniformFlat[n,basisSize,dataSize];//AbsoluteTiming(*This generates a list of b-splines sampled on the same grid as the data*)

Each row in bFlat represents a sampled B-spline:

ListLinePlot[bFlat]

Do least squares fitting

(*Manually*)
c=Inverse[bFlat.Transpose[bFlat]].bFlat.signal

(*Or via the least squares function of Mathematica*)
c=LeastSquares[Transpose[bFlat],signal]

Show the fitted signal (it coincides with the data)

(*Reconstruction*)
signalRecon=Transpose[bFlat].c
ListLinePlot[{signal,signalRecon},InterpolationOrder->0,PlotStyle->{Blue,{Orange,Dashed}}]

Show the interpolation function (the spline coincides with the sampled data at the grid points, check yourself what it looks like for other order splines)

Here the shift by 0.5 is used for plotting purposes (by default a “pixel” or data element corresponds to the all x \in [i,i+1] which we now change to [i-0.5,i+0.5] (i.e., i is the center of the pixel).

Show[
ListLinePlot[{signal,signalRecon},InterpolationOrder->0,PlotStyle->{Blue,{Orange,Dashed}}],
Plot[Sum[c[[i]]*B[n][((x-0.5)/scale)-i],{i,1,Length[c],1}],{x,1,Length[signal]},PlotStyle->Orange],
ImageSize->800
]

Fit coarse scale basis

Settings and construct basis

n=3
dataSize=Dimensions[signal]
basisSize={20}(*Now we fit with a lower number of basis functions!*)
scale=N@First[dataSize/basisSize]
bFlat=BSplineBasisUniformFlat[n,basisSize,dataSize];//AbsoluteTiming

Do the fitting

(*Manually*)
c=Inverse[bFlat.Transpose[bFlat]].bFlat.signal

(*Or via the least squares function of Mathematica*)
c=LeastSquares[Transpose[bFlat],signal]

Show the fitted signal (it coincides with the data)

(*Reconstruction*)
signalRecon=Transpose[bFlat].c
ListLinePlot[{signal,signalRecon},InterpolationOrder->0,PlotStyle->{Blue,{Orange,Dashed}}]

Show the interpolation function

Show[
ListLinePlot[{signal,signalRecon},InterpolationOrder->0,PlotStyle->{Blue,{Orange,Dashed}}],
Plot[Sum[c[[i]]*B[n][(x/scale)-i],{i,1,Length[c],1}],{x,1,Length[signal]},PlotStyle->Orange],
ImageSize->800
]

Example 2D

Load an input image

image=ImageResize[ColorConvert[ExampleData[{"TestImage","Man"}],"Grayscale"],128]
f=ImageData[image];

Expand the image in a B-spline basis (with the same nr of B-splines as there are pixels)

Settings and construction of the basis

(*Settings*)
n=3(*Spline degree*)
dataSize=Dimensions[f];
basisSize=dataSize;
(*The basis functions*)
bFlat=BSplineBasisUniformFlat[n,basisSize,dataSize];//AbsoluteTiming

Find the least squares fit via gradient descent

(*Initialize c*)
c=ConstantArray[0.,Times@@basisSize];
(*Precompute S\[Transpose]S*)
SS=bFlat.Transpose[bFlat];
(*Precompute S.y*)
Sy=bFlat.Flatten[f];
(*Compute the initial gradient*)
grad=SS.c-Sy;
(*The stepsize in gradient direction*)
Δ=0.05;
(*Iterate gradient decent a couple of times*)
Do[
c=c-Δ*grad;
grad=SS.c-Sy;
,{10000}];//AbsoluteTiming

Or directly using mathematica’s LeastSquares function

(*c=LeastSquares[Transpose[bFlat],Flatten[f]];*)

Show the reconstructed image

fRecon=ArrayReshape[Transpose[bFlat].c,dataSize];
Image[fRecon]

Up sample the image

factor=4;
bFlatUp=BSplineBasisUniformFlat[n,basisSize,factor*dataSize];//AbsoluteTiming
fReconUp=ArrayReshape[Transpose[bFlatUp].c,factor*dataSize];
Image[fReconUp]

Compare with nearest neighbor up-sampling

Row[{Image[fRecon,ImageSize->500],Image[fReconUp,ImageSize->500]}]

Expand the image in a coarse scale B-spline basis (with a lower nr of B-splines than there are pixels)

Settings and construction of the basis

(*Settings*)
n=3;(*Spline degree*)
dataSize=Dimensions[f];
basisSize=Round[dataSize/8];
(*The basis functions*)
bFlat=BSplineBasisUniformFlat[n,basisSize,dataSize];//AbsoluteTiming

Find the least squares fit via gradient descent (this is much faster than inverting the huge basis matrix)

(*Initialize c*)
c=ConstantArray[0.,Times@@basisSize];
(*Precompute S\[Transpose]S*)
SS=bFlat.Transpose[bFlat];
(*Precompute S.y*)
Sy=bFlat.Flatten[f];
(*Compute the initial gradient*)
grad=SS.c-Sy;
(*The stepsize in gradient direction*)
Δ=0.01;
(*Iterate gradient decent a couple of times*)
Do[
c=c-Δ*grad;
grad=SS.c-Sy;
,{1000}];//AbsoluteTiming

Show the reconstructed image

fRecon=ArrayReshape[Transpose[bFlat].c,dataSize];
Image[fRecon,ImageSize->500]