0

I am trying to fit a smoothing B-spline to some data and I found this very helpful post on here. However, I not only need the spline, but also its derivatives, so I tried to add the following code to the example:

tck_der = interpolate.splder(tck, n=1)
x_der, y_der, z_der = interpolate.splev(u_fine, tck_der)

For some reason this does not seem to work due to some data type issues. I get the following traceback:

Traceback (most recent call last):
  File "interpolate_point_trace.py", line 31, in spline_example
    tck_der = interpolate.splder(tck, n=1)
  File "/home/user/anaconda3/lib/python3.7/site-packages/scipy/interpolate/fitpack.py", line 657, in splder
     return _impl.splder(tck, n)
   File "/home/user/anaconda3/lib/python3.7/site-packages/scipy/interpolate/_fitpack_impl.py", line 1206, in splder
     sh = (slice(None),) + ((None,)*len(c.shape[1:]))
 AttributeError: 'list' object has no attribute 'shape'

The reason for this seems to be that the second argument of the tck tuple contains a list of numpy arrays. I thought turning the input data to be a numpy array as well would help, but it does not change the data types of tck.

Does this behavior reflect an error in scipy, or is the input malformed? I tried manually turning the list into an array:

tck[1] = np.array(tck[1])

but this (which didn't surprise me) also gave an error: 

ValueError: operands could not be broadcast together with shapes (0,8) (7,1)

Any ideas of what the problem could be? I have used scipy before and on 1D splines the splder function works just fine, so I assume it has something to do with the spline being a line in 3D.

------- edit --------

Here is a minimum working example:

import numpy as np
import matplotlib.pyplot as plt
from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D

total_rad = 10
z_factor = 3
noise = 0.1

num_true_pts = 200
s_true = np.linspace(0, total_rad, num_true_pts)
x_true = np.cos(s_true)
y_true = np.sin(s_true)
z_true = s_true / z_factor

num_sample_pts = 80
s_sample = np.linspace(0, total_rad, num_sample_pts)
x_sample = np.cos(s_sample) + noise * np.random.randn(num_sample_pts)
y_sample = np.sin(s_sample) + noise * np.random.randn(num_sample_pts)
z_sample = s_sample / z_factor + noise * np.random.randn(num_sample_pts)

tck, u = interpolate.splprep([x_sample, y_sample, z_sample], s=2)
x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
u_fine = np.linspace(0, 1, num_true_pts)
x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck)

# this is the part of the code I inserted: the line under this causes the crash
tck_der = interpolate.splder(tck, n=1)
x_der, y_der, z_der = interpolate.splev(u_fine, tck_der)

# end of the inserted code

fig2 = plt.figure(2)
ax3d = fig2.add_subplot(111, projection='3d')
ax3d.plot(x_true, y_true, z_true, 'b')
ax3d.plot(x_sample, y_sample, z_sample, 'r*')
ax3d.plot(x_knots, y_knots, z_knots, 'go')
ax3d.plot(x_fine, y_fine, z_fine, 'g')
fig2.show()
plt.show()
meetaig
  • 913
  • 10
  • 26

2 Answers2

2

Stumbled into the same problem...

I circumvented the error by using interpolate.splder(tck, n=1) and instead used interpolate.splev(spline_ev, tck, der=1) which returns the derivatives at the points spline_ev (see Scipy Doku).

If you need the spline I think you can then use interpolate.splprep() again.

In total something like:

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

points = np.random.rand(10,2) * 10

(tck, u), fp, ier, msg = interpolate.splprep(points.T, s=0, k=3, full_output=True)

spline_ev = np.linspace(0.0, 1.0, 100, endpoint=True)

spline_points = interpolate.splev(spline_ev, tck)

# Calculate derivative
spline_der_points = interpolate.splev(spline_ev, tck, der=1)
spline_der = interpolate.splprep(spline_der_points.T, s=0, k=3, full_output=True)


# Plot the data and derivative
fig = plt.figure()

plt.plot(points[:,0], points[:,1], '.-', label="points")
plt.plot(spline_points[0], spline_points[1], '.-', label="tck")
plt.plot(spline_der_points[0], spline_der_points[1], '.-', label="tck_der")

#   Show tangent
plt.arrow(spline_points[0][23]-spline_der_points[0][23], spline_points[1][23]-spline_der_points[1][23], 2.0*spline_der_points[0][23], 2.0*spline_der_points[1][23])

plt.legend()

plt.show()

EDIT:

I also opened an Issue on Github and according to ev-br the usage of interpolate.splprep is depreciated and one should use make_interp_spline / BSpline instead.

qwertz
  • 46
  • 8
0

As noted in other answers, splprep output is incompatible with splder, but is compatible with splev. And the latter can evaluate the derivatives.

However, for interpolation, there is an alternative approach, which avoids splprep altogether. I'm basically copying a reply on the SciPy issue tracker (https://github.com/scipy/scipy/issues/10389):

Here's an example of replicating the splprep outputs. First let's make sense out of the splprep output:

# start with the OP example
import numpy as np
from scipy import interpolate

points = np.random.rand(10,2) * 10

(tck, u), fp, ier, msg = interpolate.splprep(points.T, s=0, k=3, full_output=True)

# check the meaning of the `u` array: evaluation of the spline at `u`
# gives back the original points (up to a list/transpose)

xy = interpolate.splev(u, tck)
xy = np.asarray(xy)
np.allclose(xy.T, points)

Next, let's replicate it without splprep. First, build the u array: the curve is represented parametrically, and u is essentially an approximation for the arc length. Other parametrizations are possible, but here let's stick to what splprep does. Translating the pseudocode from the doc page, https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.splprep.html

vv = np.sum((points[1:, :] - points[:-1, :])**2, axis=1)
vv = np.sqrt(vv).cumsum()
vv/= vv[-1]
vv = np.r_[0, vv]

# check: 
np.allclose(u, vv)

Now, interpolate along the parametric curve: points vs vv:

spl = interpolate.make_interp_spline(vv, points)

# check spl.t vs knots from splPrep
spl.t - tck[0]

The result, spl, is a BSpline object which you can evaluate, differentiate etc in a usual way:

np.allclose(points, spl(vv))

# differentiate
spl_derivative = spl.derivative(vv)
ev-br
  • 24,968
  • 9
  • 65
  • 78