Search code examples
pythonscipyinterpolationspline

What is the difference between the various spline interpolators from Scipy?


My goal is to calculate a smooth trajectory passing through a set of points as shown below. I have had a look at the available methods of scipy.interpolate, and also the scipy user guide. However, the choice of the right method is not quite clear to me.

What is the difference between BSpline, splprep, splrep, UnivariateSpline, interp1d, make_interp_spline and CubicSpline?

According to the documentation, all functions compute some polynomial spline function through a sequence of input points. Which function should I select? A) What is the difference between a cubic spline and a B-spline of order 3? B) What is the difference between splprep() and splrep()? C) Why is interp1d() going to be deprecated?

I know I'm asking different questions here, but I don't see the point in splitting the questions up as I assume the answers will be related.

All in all, I find that the scipy.interpolate module is organized a little confusingly. I thought maybe I'm not the only one who has this impression, which is why I'm reaching out to SO.


Here's how far I've come. Below is some code that runs the different spline functions for some test data. It creates the figure below.

  1. I've read somewhere: "All cubic splines can be represented as B-splines of order 3", and that "it's a matter of perspective which representation is more convenient". But why do I end up in different results if I use CublicSpline and any of the B-spline methods?
  2. I found that it is possible to construct a BSpline object from the output of splprep and splrep, such that the results of splev() and BSpline() are equivalent. That way, we can convert the output of splrep and splprep into the object-oriented interface of scipy.interpolate().
  3. splrep, UnivariateSpline and make_interp_spline lead to the same result. In my 2D data, I need to apply the interpolation independently per data dimension that it works. The convenience function interp1d yields the same result, too. Related SO-question: Link
  4. splprep and splrep seem unrelated. Even if I compute splprep twice for every data axis independently (see p0_new), the result looks differently. I see in the docs that splprep computes the B-spline representation of an n-D curve. But should splrep and splprep not be related?
  5. splprep, splrep and UnivariateSpline have a smoothing parameter, while other interpolators have no such parameter.
  6. splrep pairs with UnivariateSpline. However, I couldn't find a matching object-oriented counterpart for splprep. Is there one?

Comparison of different spline methods

import numpy as np
from scipy.interpolate import *
import matplotlib.pyplot as plt

points = [[0, 0], [4, 4], [-1, 9], [-4, -1], [-1, -9], [4, -4], [0, 0]]
points = np.asarray(points)

n = 50 
ts = np.linspace(0, 1, len(points))
ts_new = np.linspace(0, 1, n)

(t0_0,c0_0,k0_0), u = splprep(points[:,[0]].T, s=0, k=3)
(t0_1,c0_1,k0_1), u = splprep(points[:,[1]].T, s=0, k=3)
p0_new = np.r_[np.asarray(splev(ts_new, (t0_0,c0_0,k0_0))),
               np.asarray(splev(ts_new, (t0_1,c0_1,k0_1))),
                ].T

# splprep/splev
(t1,c1,k1), u = splprep(points.T, s=0, k=3)
p1_new = splev(ts_new, (t1,c1,k1))
# BSpline from splprep
p2_new = BSpline(t1, np.asarray(c1).T, k=k1)(ts_new)
# splrep/splev (per dimension)
(t3_0,c3_0,k3_0) = splrep(ts, points[:,0].T, s=0, k=3)
(t3_1,c3_1,k3_1) = splrep(ts, points[:,1].T, s=0, k=3)
p3_new = np.c_[splev(ts_new, (t3_0,c3_0,k3_0)),
               splev(ts_new, (t3_1,c3_1,k3_1)),
               ]
# Bspline from splrep
p4_new = np.c_[BSpline(t3_0, np.asarray(c3_0), k=k3_0)(ts_new),
               BSpline(t3_1, np.asarray(c3_1), k=k3_1)(ts_new),
               ]
# UnivariateSpline
p5_new = np.c_[UnivariateSpline(ts, points[:,0], s=0, k=3)(ts_new),
               UnivariateSpline(ts, points[:,1], s=0, k=3)(ts_new),]
# make_interp_spline
p6_new = make_interp_spline(ts, points, k=3)(ts_new)
# CubicSpline
p7_new = CubicSpline(ts, points, bc_type="clamped")(ts_new)
# interp1d
p8_new = interp1d(ts, points.T, kind="cubic")(ts_new).T

fig, ax = plt.subplots()
ax.plot(*points.T, "o-", label="Original points")
ax.plot(*p1_new,   "o-", label="1: splprep/splev")
ax.plot(*p2_new.T, "x-", label="1: BSpline from splprep")
ax.plot(*p3_new.T, "o-", label="2: splrep/splev")
ax.plot(*p4_new.T, "x-", label="2: BSpline from splrep")
ax.plot(*p5_new.T, "*-", label="2: UnivariateSpline")
ax.plot(*p6_new.T, "+-", label="2: make_interp_spline")
ax.plot(*p7_new.T, "x-", label="3: CubicSpline")
#ax.plot(*p8_new.T, "k+-", label="3: interp1d")
#ax.plot(*p0_new.T, "k+-", label="3: CubicSpline")
ax.set_aspect("equal")
ax.grid("on")
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

Solution

  • What is the difference between BSpline, splprep, splrep, UnivariateSpline, interp1d, make_interp_spline and CubicSpline?

    • a BSpline object represents a spline function in terms of knots t, coefficients c and degree k. It does not know anything about the data x, y. It is a low-level implementation object, on par with PPoly --- think of PPoly vs BSpline as a change of basis.

    • make_interp_spline constructs a spline (BSpline) which passes through the data x and y -- it's an interpolating spline, so that spl(x) == y. It can handle batches of data: if y is two-dimensional and the second dimension has length n, make_interp_spline(x, y) represents a stack of functions y_1(x), y_2(x), ..., y_n(x).

    • CubicSpline is similar but more limited: it only allows cubic splines, k=3. If you ever need to switch from cubics to e.g. linear interpolation, with make_interp_spline you just change k=3 to k=1 at the call site. OTOH, CubicSpline has some features BSpline does not: e.g. root-finding. It is a PPoly instance, not BSpline.

    • interp1d is not deprecated, it is legacy. A useful part of interp1d is nearest/previous/next modes. The rest just delegates to make_interp_spline, so better use that directly.

    • splrep constructs a smoothing spline function given data. The amount of smoothing is controlled by the s parameter, with s=0 being interpolation. It returns not a BSpline but a tck tuple (knots, coefficients and degree).

    • splprep constructs a smoothing spline curve in a parametric form. That is (x(u), y(u)) not y(x). Also returns tck-tuples.

    • UnivariateSpline is equivalent to splrep.

    Note that scipy.interpolate has grown organically. With at least four generations of developers over almost quarter century. And the FITPACK library which splrep, splprep and UnivariateSpline delegate to is from 1980s. CubicSpline, BSpline and make_interp_spline are not using FITPACK.

    So in new code I personally would recommend make_interp_spline + BSpline; if you only need cubics, CubicSpline works, too.

    But why do I end up in different results if I use CublicSpline and any of the B-spline methods?

    Boundary conditions. Both make_interp_spline and CubicSpline default to not-a-knot, and you can change it. splrep et al only use not-a-knot.