Scipy Spline Interpolation Parameter

Question

Documentation in scipy.interpolate (found at https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) states:

"The parameter variable is given with the keyword argument, u, which defaults to an equally-spaced monotonic sequence between 0 and 1."

From experimenting with the function in 2D, it appears that the parameter is the length along the curve. Here's a snippet my test script:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splprep, splev

n = 20

iota = np.arange(n)
gamma = 5
theta = np.pi/2 * (1.0 * iota / (n-1)) ** gamma

r = np.cos(theta)
z = np.sin(theta)

spl, u = splprep([r,z])
spline = splev(u, spl)

# Sampled point at half-length?
test = splev(0.5, spl)

plt.plot(r, z, '*')
plt.plot(spline[0], spline[1], '*', color = 'red')
plt.plot(test[0], test[1], '*', color = 'cyan')

# Exact mid-point.
plt.plot(np.sqrt(2)/2, np.sqrt(2)/2, '*', color = 'purple')
plt.axis('equal')
plt.show()

The above statement does not (to me) seem to clearly convey this behavior. Would anyone have input on this?

Answer 1

The parameter variable is given with the keyword argument, $u$, which defaults to an equally-spaced monotonic sequence between...

When calling the function splprep you may pass an argument of $u$. This argument $u$ should be the same length as your data, and the algorithm will treat it as the parameter of the curve, sometimes called the 'time' along the curve. Later, when using splev, the argument x you supply to it, is the 'time' for which you want the interpolated value.

I edited your code a bit to try experimenting with this feature. You see that upon changing $beta$, you would get overshoots or undershoots in your interpolation.

I also think that the documentation you linked to is incorrect, because setting beta to 1 does not yield the same result as omitting the $u$ argument to splprep. In fact, I tested it, and the default $u$ is cumulative sum of the segment lengths. The code below prints False:

_, u_default = splprep([r,z])
C = np.cumsum((np.diff(np.array([r,z]).transpose(),axis=0)**2).sum(axis=1)**0.5)
C /= C[-1]
C = np.concatenate(([0],C))
print(np.any(u_default-C))

Some more experiments with the $u$ argument:

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splprep, splev
plt.figure()

n = 20

iota = np.arange(n)
gamma = 5
beta = 2
theta = np.pi/2 * (1.0 * iota / (n-1)) ** gamma

r = np.cos(theta)
z = np.sin(theta)

u_arg = np.linspace(0,1,n)**beta
spl, u = splprep([r,z],u=u_arg)
print(np.any(u_arg!=u))
spline = splev(u, spl)

# Sampled point at half-length?
test = splev(0.5, spl)

plt.plot(r, z, 'v','b')
plt.plot(spline[0], spline[1], '^', color = 'red')
plt.plot(test[0], test[1], 'x', color = 'cyan')

# Exact mid-point.
plt.plot(np.sqrt(2)/2, np.sqrt(2)/2, '+', color = 'purple')
plt.axis('equal')
plt.show()