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.
- 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
CublicSplineand any of the B-spline methods? - I found that it is possible to construct a
BSplineobject from the output ofsplprepandsplrep, such that the results ofsplev()andBSpline()are equivalent. That way, we can convert the output ofsplrepandsplprepinto the object-oriented interface ofscipy.interpolate(). splrep,UnivariateSplineandmake_interp_splinelead to the same result. In my 2D data, I need to apply the interpolation independently per data dimension that it works. The convenience functioninterp1dyields the same result, too. Related SO-question: Linksplprepandsplrepseem unrelated. Even if I computesplpreptwice 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 shouldsplrepandsplprepnot be related?splprep,splrepandUnivariateSplinehave a smoothing parameter, while other interpolators have no such parameter.splreppairs withUnivariateSpline. However, I couldn't find a matching object-oriented counterpart forsplprep. Is there one?
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()

a
BSplineobject represents a spline function in terms of knotst, coefficientscand degreek. It does not know anything about the datax, y. It is a low-level implementation object, on par withPPoly--- think of PPoly vs BSpline as a change of basis.make_interp_splineconstructs a spline (BSpline) which passes through the dataxandy-- it's an interpolating spline, so thatspl(x) == y. It can handle batches of data: if y is two-dimensional and the second dimension has lengthn,make_interp_spline(x, y)represents a stack of functionsy_1(x), y_2(x), ..., y_n(x).CubicSplineis similar but more limited: it only allows cubic splines,k=3. If you ever need to switch from cubics to e.g. linear interpolation, withmake_interp_splineyou just changek=3tok=1at the call site. OTOH,CubicSplinehas 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.
splrepconstructs a smoothing spline function given data. The amount of smoothing is controlled by thesparameter, withs=0being interpolation. It returns not a BSpline but a tck tuple (knots, coefficients and degree).splprepconstructs a smoothing spline curve in a parametric form. That is(x(u), y(u))noty(x). Also returns tck-tuples.UnivariateSplineis equivalent tosplrep.Note that
scipy.interpolatehas grown organically. With at least four generations of developers over almost quarter century. And the FITPACK library whichsplrep,splprepandUnivariateSplinedelegate to is from 1980s.CubicSpline,BSplineandmake_interp_splineare not using FITPACK.So in new code I personally would recommend
make_interp_spline+BSpline; if you only need cubics,CubicSplineworks, too.Boundary conditions. Both
make_interp_splineandCubicSplinedefault tonot-a-knot, and you can change it.splrepet al only usenot-a-knot.