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
1.1.0.0 | Added mex acceleration, removed separate uniform B-spline implementation |
Inspired by: mexme - write MEX files in no time
Create scripts with code, output, and formatted text in a single executable document.
johanna ganglbauer (view profile)
thanks 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)
gg (view profile)
I 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 Yeh (view profile)
Repeated knot (or knot with multiplicity larger than 1) doesn't seem to work using this class.
Jakob Ameres (view profile)
For 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 golk (view profile)
How to obtain for 3D spline tensor products? would make a sample? thanks
Kobye (view profile)
This 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!
Kobye (view profile)
Ben (view profile)
@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 Mineault (view profile)
These are 1d splines. You can obtain 2d/3d splines using tensor products of the basis given by the class.
Arso (view profile)
Hi,
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 Obler (view profile)
Patrick Mineault (view profile)
Richard 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 Obler (view profile)
Like Martijn I cannot compile on a 64bit machine. Same error messages.
Martijn (view profile)
I 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...