version 1.1.0.0 (11.6 KB) by
Patrick Mineault

Fit, evaluate, differentiate non-uniform B-splines of any order - fast

fastBSpline - A fast, lightweight class that implements

non-uniform B splines of any order

Matlab's spline functions are very general. This generality comes at

the price of speed. For large-scale applications, including model

fitting where some components of the model are defined in terms of

splines, such as generalized additive models, a faster solution is

desirable.

The fastBSpline class implements a lightweight set of B-spline

features, including evaluation, differentiation, and parameter fitting.

The hard work is done by C code, resulting in up to 10x acceleration

for evaluating splines and up to 50x acceleration when evaluating

of spline derivatives.

Nevertheless, fastBSplines are manipulated using an intuitive, high-

level object-oriented interface, thus allowing C-level performance

without the messiness. Use CompileMexFiles to compile the required

files. If mex files are not available, evaluation will be done in .m

code, so you may still use the code if you can't use a compiler for your

platform.

B splines are defined in terms of basis functions:

y(x) = sum_i B_i(x,knots)*weights_i

B (the basis) is defined in terms of knots, a non-decreasing sequence

of values. Each basis function is a piecewise polynomial of order

length(knots)-length(weights)-1. The most commonly used B-spline is

the cubic B-spline. In that case there are 4 more knots than there

are weights. Another commonly used B-spline is the linear B-spline,

whose basis function are shaped like tents, and whose application

results in piecewise linear interpolation.

The class offers two static functions to fit the weights of a spline:

lsqspline and pspline. It includes facilities for computing the basis

B and the derivatives of the spline at all points.

Constructor:

sp = fastBSpline(knots,weights);

Example use:

%Fit a noisy measurement with a smoothness-penalized spline (p-spline)

x = (0:.5:10)';

y = sin(x*pi*.41-.9)+randn(size(x))*.2;

knots = [0,0,0,0:.5:10,10,10,10];

%Notice there are as many knots as observations

%Because there are so many knots, this is an exact interpolant

sp1 = fastBSpline.lsqspline(knots,3,x,y);

%Fit penalized on the smoothness of the spline

sp2 = fastBSpline.pspline(knots,3,x,y,.7);

clf;

rg = -2:.005:12;

plot(x,y,'o',rg,sp1.evalAt(rg),rg,sp2.evalAt(rg));

legend('measured','interpolant','smoothed');

fastBSpline properties:

outOfRange - Determines how the spline is extrapolated outside the

range of the knots

knots - The knots of the spline (read only)

weights - The weights of the spline (read only)

fastBSpline Methods:

fastBSpline - Construct a B spline from weights at the knots

lsqspline - Construct a least-squares spline from noisy measurements

pspline - Construct a smoothness-penalized spline from noisy

measurements

evalAt - Evaluate a spline at the given points

getBasis - Get the values of the underlying basis at the given points

Btimesy - Evaluate the product getBasis(x)'*y

dx - Returns another fastBSpline object which computes the derivative of

the original spline

Disclaimer: fastBSpline is not meant to replace Matlab's spline functions;

it does not include any code from the Mathworks

Patrick Mineault (2021). Fast B-spline class (https://www.mathworks.com/matlabcentral/fileexchange/32509-fast-b-spline-class), MATLAB Central File Exchange. Retrieved .

Created with
R2009b

Compatible with any release

**Inspired by:**
mexme - write MEX files in no time

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Daniel DurhamNice, I'm seeing ~2x faster evaluation and the getBasis functionality is MUCH faster, thank you!

johanna ganglbauerthanks for the great contribution!

A few comments on how to use the tool:

1) boundary conditions: i am not quite sure how the boundary conditions are set in this code, but the problem of oscillations on the edges can be circumvented by providing many data points close to the edges and using lsqspline which makes a good approximation using least squares. The knot sequence should not be denser on the edges.

2) extension to 2 dimensions: is well explained by Les Piegl and Wayne Tiller in 'The nurbs book' in chapter 9.2.5

basically the interpolation Q_k,l=N_i(x)*N_j(y)*P_i,j (sum over i and j) with Q_k,l being the control points, N the basis functions and P_i,j the weight matrix has to be split into two parts:

first of all coefficients R_i,l are found by interpolation in x direction for each fixed y-value y_l:

Q_k,l=N_i(x)*R_i,l (sum over i)

Then coefficients are interpolated again in y direction leading to the coefficient P_i,j with

R_i,l=N_j(y)*P_i,j (sum over j)

ggI would expect the variable tt in the example below to be a vector of all ones but the last entry is zero. Is that a bug or a feature?

knots = [0, 0, 0, 0, 1, 1, 1, 1];

weights = [1, 1, 1, 1];

sp = fastBSpline(knots,weights);

testx = linspace(0,1,100);

tt = sum(getBasis(sp, testx), 2);

all(tt==1) % I would expect TRUE but this is FALSE

tt(end) == 0 % Last value is the culprit! All others are fine.

Raine YehRepeated knot (or knot with multiplicity larger than 1) doesn't seem to work using this class.

Jakob AmeresFor Compilation under Linux it might be useful:

function [] = CompileMexFiles()

%Compiles mex files that accelerate B-spline evaluation

mex -v CFLAGS="$CFLAGS -std=c99" evalBin.c

mex -v CFLAGS="$CFLAGS -std=c99" evalBSpline.c

mex -v CFLAGS="$CFLAGS -std=c99" evalBinTimesY.c

end

Ehsan golkHow to obtain for 3D spline tensor products? would make a sample? thanks

KobyeThis class is excellent. The compilation on my 64 bit machine running Matlab R2012b worked after (1) downloading Microsoft SDK 7.1 (2) renaming all files ending in .c to .ccp, both in the directory and all references contained within the .c files. Great work!

KobyeBen@Patrick

Could you kindly provide a sample code for 2D/3D spline using tensor product. I believe that will further promote your code and enlarge your contributions to the society. Thank you!

Patrick MineaultThese are 1d splines. You can obtain 2d/3d splines using tensor products of the basis given by the class.

ArsoHi,

This is 2D case only?

Is there parametrized interpolation like evalAt(t=0...1)?

It fails in general, for messy points x=[0; 6 ; 3 ; 12; 8] and y=[0;-3; 6; 12; 0]?

If I am not wrong, it works similar as function provided with Matlab...

If You could extend such possibilities in 3D - it would be great, because Matlab is lacking with CAD features.

Richard OblerPatrick MineaultRichard has sent me an email with this solution for his problems which stem from using a different compiler (VC++ on Windows rather than gcc on Linux):

"Under Windows it is important that all files have the ending .cpp instead of .c so that the right compiler is taken. After correcting that (*.c à *.cpp) and stepping through all calls of .c files in the code which also have to be corrected, of course, everything works perfect!"

Richard OblerLike Martijn I cannot compile on a 64bit machine. Same error messages.

MartijnI get a long list of syntax errors when trying to compile on x64. I have not tried 32.

evalBin.c

evalBin.c(98) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(99) : error C2065: 'x_ptr' : undeclared identifier

evalBin.c(99) : warning C4047: 'function' : 'const mxArray *' differs in levels of indirection from 'int'

evalBin.c(99) : warning C4024: 'mexmetypecheck' : different types for formal and actual parameter 1

evalBin.c(100) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(101) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(102) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(103) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(104) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(105) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(106) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(107) : error C2143: syntax error : missing ';' before 'const'

evalBin.c(108) : error C2065: 'firstknot_ptr' : undeclared identifier

etc...