Scientific Programming with B-Splines

bsplines.org September 16, 2019

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 floats, not as ints. Otherwise you may get wrong results.

Functions for directly evaluating B-splines are located in scipy.signal, for example:

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}\).

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\).

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sps
import scipy.interpolate as spi

# 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:

>>> import diofant
>>> import diofant.functions.special.bsplines as bs
>>> 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.

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, 4
p <- 3
x <- c(0, p+1)
xi <- 0:(p+1)
b <- bs(x, knots=xi, degree=p, intercept=TRUE)

# plot B-splines
xx <- 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.

extern crate bspline;
extern crate gnuplot;

fn main() {
  // set B-spline degree
  let p = 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")
  let mut c:  std::vec::Vec<f32> =
      std::vec::Vec::with_capacity(2*p+1);
  let mut xi: std::vec::Vec<f32> =
      std::vec::Vec::with_capacity(3*p+2);

  for k in 0..3*p+2 {
    if k < p {
      c.push(0.0);
      xi.push(0.0);
    } else if k < 2*p+1 {
      c.push((k == p) as i32 as f32);
      xi.push((k-p) as f32);
    } else {
      xi.push((p+1) as f32);
    }
  }

  // construct B-spline
  let b = bspline::BSpline::new(p, c, xi);

  // evaluate B-spline for plotting
  let nn = 100;
  let mut xx: std::vec::Vec<f32> =
      std::vec::Vec::with_capacity(nn as usize);
  let mut yy: std::vec::Vec<f32> =
      std::vec::Vec::with_capacity(nn as usize);

  for i in 0..nn {
    let x = ((p+1) as f32) * (i as f32) / ((nn-1) as f32);
    let y = b.point(x);
    xx.push(x);
    yy.push(y);
  }

  // plot via gnuplot
  let mut fig = 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.