B-splines are already included in the commonly used libraries of most scientific scripting languages. If this is not the case for your favorite language, then it is not difficult to implement the necessary formulas such as the recurrence relation by Cox and de Boor.

For each of the following languages, not all relevant B-spline functions are listed. However, the given functions are linked to the respective documentation and should serve you as a reasonable starting point.

## Python

SciPy contains quite extensive (B-)spline functionality in its two modules `scipy.signal`

and `scipy.interpolate`

. *Caution:* When evaluating (B-)splines, always give the evaluation points \(x\) as `float`

s, not as `int`

s. Otherwise you may get wrong results.

Functions for directly evaluating B-splines are located in `scipy.signal`

, for example:

`y = bspline(x, p)`

evaluates the centralized B-spline \(y_i = b^p(x_i + \frac{p+1}{2})\) of degree \(p\).`y = quadratic(x)`

,`y = cubic(x)`

are special cases of`bspline`

for \(p = 2, 3\).

Functions for calculating with splines as linear combinations of B-splines are in `scipy.interpolate`

:

`s = BSpline(xi, c, p)`

constructs the spline \(s = \sum_{k=0}^{m-1} c_k b_{k,\vec{\xi}}^p\) of degree \(p\) with knots \(\vec{\xi} = (\xi_0, \dotsc, \xi_{m+p})\) and coefficients \(\vec{c} = (c_0, \dotsc, c_{m-1})\).`BSpline`

was recently introduced and requires at least version 0.19.0 of SciPy.`y = s(x, nu=q)`

evaluates the spline (default without`nu`

`) or its \(q\)-th derivative, i.e., \(y_i = s^{(q)}(x_i)\).`

`r = s.derivative(nu=q)`

constructs the \(q\)-th derivative (default: first derivative) of \(s\) as a`BSpline`

object \(r = s’\).`S = s.antiderivative(nu=q)`

constructs a \(q\)-th antiderivative (default: first antiderivative) of \(s\) as a`BSpline`

object \(S = \int s(x) \mathsf{d}x\).`I = s.integrate(a, b)`

computes the integral \(I = \int_a^b s(x) \mathsf{d}x\).

`s = UnivariateSpline(x, y, k=p, w=w, s=s)`

constructs a spline \(s\) of degree \(p \le 5\) trying to fit the data \((x_i, y_i)\). The knots are automatically chosen and determined by the smoothing parameter \(s \ge 0\); the larger \(s\) is, the smoother the spline will be. Setting \(s = 0\) will interpolate the data points. Optionally, the data points can be weighted with the vector \(\vec{w}\).`s = InterpolatedUnivariateSpline(x, y, k=p, w=w)`

is a shortcut for`UnivariateSpline`

with \(s = 0\).`s = LSQUnivariateSpline(x, y, xi, k=p, w=w)`

is like`UnivariateSpline`

, but lets you choose the knots \(\vec{\xi}\) yourself.

`scipy.interpolate`

contains analogous classes and functions for bivariate splines. Use `RectBivariateSpline`

if your data is on a regular grid. For scattered data, use `SmoothBivariateSpline`

and `LSQBivariateSpline`

instead.

As an example, the following code first plots the cubic cardinal B-spline with knots \(0, 1, 2, 3, 4\) and then a cubic non-uniform spline as a linear combination of \(m = 5\) B-splines with coefficients \(c\).

importnumpyasnpimportmatplotlib.pyplotaspltimportscipy.signalasspsimportscipy.interpolateasspi# plot cubic cardinal B-spline (knots 0, 1, 2, 3, 4)p = 3 xx = np.linspace(0, p+1, 100) yy = sps.bspline(xx - (p+1)/2, p) plt.plot(xx, yy) plt.show()# plot cubic non-uniform spline (m=5 DOFs)xi = [0, 1, 3, 4, 6, 7, 8, 10, 11] c = [2, -1, 1, 0, 1] s = spi.BSpline(xi, c, p) m = len(c) xx = np.linspace(xi[p], xi[m]) yy = s(xx) plt.plot(xx, yy) plt.show()

It should be mentioned that SciPy is not the only Python library with B-splines in it. If you prefer to calculate your results symbolically using a CAS (computer algebra system), you can do so with SymPy (and, thus, SageMath) or its fork Diofant. For example, use the function `b = bspline_basis(p, xi, k, x)`

in `diofant.functions.special.bsplines`

to obtain \(b_{k,\vec{\xi}}^p(x)\) (where \(x\) is a symbolic variable), i.e., for the standard cardinal cubic B-spline:

>>> importdiofant>>> importdiofant.functions.special.bsplinesasbs>>>x = diofant.symbols("x")>>>bs.bspline_basis(3, range(5), 0, x) Piecewise((x**3/6, And(x < 1, x >= 0)), (-x**3/2 + 2*x**2 - 2*x + 2/3, And(x < 2, x >= 1)), (x**3/2 - 4*x**2 + 10*x - 22/3, And(x < 3, x >= 2)), (-x**3/6 + 2*x**2 - 8*x + 32/3, And(x <= 4, x >= 3)), (0, true))

## R

R ships with the `splines`

package, which features functions to generate and evaluate (“predict” in R terminology) B-splines:

`b <- bs(x, knots=zeta, degree=p, intercept=TRUE, Boundary.knots=sigma)`

constructs the \(m\) B-splines of degree \(p \ge 1\) with knots $$\vec{\xi} = (\overbrace{\sigma_0, \dotsc, \sigma_0}^{\text{$(p+1)$ times}}, \zeta_0, \dotsc, \zeta_{m-p-2}, \overbrace{\sigma_1, \dotsc, \sigma_1}^{\text{$(p+1)$ times}}),$$ where \(m\) is determined by the dimension of the vector \(\vec{\zeta} \in \mathbb{R}^{m-p-1}\) of the interior knots. The pair \(\vec{\sigma} \in \mathbb{R}^2\) contains the boundary knots, defaulting to the range \(\vec{\sigma} = (\min(x), \max(x))\). The function`bs`

returns a matrix of size \(N \times m\), where \(\vec{x} \in \mathbb{R}^N\) and the \((i,k)\)-th entry equals \(b_{k,\vec{\xi}}^p(x_i)\). If`intercept`

is set to`FALSE`

, then the first basis function (column) is omitted from the matrix.`b <- bs(x, df=m, degree=p, intercept=TRUE, Boundary.knots=sigma)`

sets \(m\) explicitly and chooses the interior knots automatically as appropriate quantiles of \(x\).`y <- predict(b, x)`

re-evaluates the B-splines of \(b\) at \(\vec{x}\), returning a \((N \times m)\)-matrix as before with \(\vec{x} \in \mathbb{R}^N\).

The following example shows how to plot the seven cubic B-splines corresponding to the knot vector \(\vec{\xi} = (0, 0, 0, 0, 1, 2, 3, 4, 4, 4, 4)\). The “middle” B-spline with index \(k = p = 3\) is the cardinal B-spline with knots \(0, 1, 2, 3, 4\).

require("splines") require("graphics")# construct cubic B-splines# including the cardinal one with knots 0, 1, 2, 3, 4p <- 3 x <- c(0, p+1) xi <- 0:(p+1) b <- bs(x, knots=xi, degree=p, intercept=TRUE)# plot B-splinesxx <- seq(0, p+1, length=100) yy <- predict(b, xx) matplot(xx, yy, type="l")

If you want to calculate derivatives or integrals of B-splines, you may want to have a look at the `splines2`

package in CRAN. The included `bSpline`

function also generalizes the `bs`

function to the case \(p = 0\) of piecewise constant B-splines.

## MATLAB

MATLAB contains a lot of (B-)spline functionality. If you only want to work with cubic splines and you do not need B-splines or splines of other degrees, then the MATLAB core (without any toolboxes) suffices. The function `mkpp`

creates spline objects, given knots and coefficients of the polynomial pieces. The objects can be evaluated with `ppval`

or the underlying data can be read out by `unmkpp`

. For 1D cubic spline interpolation, the spline objects can be calculated via `spline`

. This is also possible using the more general function `interp1`

, or, for 2D and 3D data, with `interp2`

and `interp3`

. The most general method is using `griddedInterpolant`

for \(d\)-dimensional data on a regular cartesian grid. However, other newer functions such as `fillmissing`

can employ spline interpolation, too.

If you want to work directly on the B-spline basis or want to employ sophisticated B-spline methods, you will need the Curve Fitting Toolbox. Historically, the functions were contained in the Spline Toolbox, which was written by Carl de Boor himself. The code got merged into the Curve Fitting Toolbox in 2008.

`b = bspline(xi)`

returns the spline object corresponding to the B-spline with knots \(\vec{\xi} \in \mathbb{R}^{p+2}\) (degree \(p\)). Omit the output argument \(b\) to plot the B-spline together with its polynomial pieces.`y = fnval(b, x)`

evaluates the B-spline \(y_i = b(x_i)\).`a = fnder(b)`

calculates the spline object for the derivative \(a = b’\).`B = fnint(b)`

calculates the spline object for an antiderivative \(B = \int b(x) \mathsf{d}x\).

`bspligui()`

brings up an interactive GUI to experiment with the effect of the knots on the corresponding B-spline.`splinetool(x, y)`

starts an interactive GUI for experimenting with different 1D spline approximation methods for the data \(\vec{x}, \vec{y}\). Omit \(\vec{x}, \vec{y}\) to use included example data.

As a last remark: If you want to use these more advanced B-spline methods, but do not have access to the Curve Fitting Toolbox and you do not want to write the code yourself, the MATLAB File Exchange could help you out, for example this Spline Toolbox.

## Julia

Have a look at the Interpolations.jl package, which seems to be currently the most recently updated Julia package for B-splines. There are also numerous other packages (ApproXD.jl, BSplines.jl, Dierckx.jl, Splines.jl, …) which are not quite up-to-date, but their implementation could serve you as starting point for writing your own B-spline code.

## Rust

Rust is a relatively new programming language, but that does not mean that no B-spline libraries would be available. There is a simple B-spline curve library named bspline available to be downloaded via Cargo. Together with the gnuplot package, it is possible to write a little example program to plot the cardinal B-spline \(b^p\) of degree \(p\).

*Disclaimer: I am a Rust novice, which means that the following code may be a bit complicated unnecessarily.*

externcratebspline;externcrategnuplot;fnmain() {// set B-spline degreeletp = 3;// construct coefficient and knot vectors// to be c = [0, ..., 0, 1, 0, ..., 0] and// xi = [0, ..., 0, 0, 1, ..., p+1, p+1, ..., p+1]// (all of the "0, ..., 0" and "p+1, ..., p+1"// mean "repeat this number p times")letmutc: std::vec::Vec<f32> = std::vec::Vec::with_capacity(2*p+1);letmutxi: std::vec::Vec<f32> = std::vec::Vec::with_capacity(3*p+2);forkin0..3*p+2 {ifk < p { c.push(0.0); xi.push(0.0); }elseifk < 2*p+1 { c.push((k == p)asi32asf32); xi.push((k-p)asf32); }else{ xi.push((p+1)asf32); } }// construct B-splineletb = bspline::BSpline::new(p, c, xi);// evaluate B-spline for plottingletnn = 100;letmutxx: std::vec::Vec<f32> = std::vec::Vec::with_capacity(nnasusize);letmutyy: std::vec::Vec<f32> = std::vec::Vec::with_capacity(nnasusize);foriin0..nn {letx = ((p+1)asf32) * (iasf32) / ((nn-1)asf32);lety = b.point(x); xx.push(x); yy.push(y); }// plot via gnuplotletmutfig = gnuplot::Figure::new(); fig.axes2d().lines(&xx, &yy, &[]); fig.show(); }

## C

This content is to be added. Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat. Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur. Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.