Determining Value of B-Spline at an Arbitrary X Coordinate

Question

I have created an open, clamped, cubic, b-spline on the set of data points xp and yp.

The spline is parametarized by the vector u which spans the domain of xp.

My goal is to determine the "y" coordinate of the b-spline at a given "x" coordinate in the domain of xp.

As is expected behavior when generating parametric curves, when I pass the value "4" into splev after calculating tck, the value for both the x and y coordinate corresponding to the parameter 4 is returned.

I am able to use Newton's method to determine the value of the parameter u at a given "x" coordinate; This is however indirect and requires more computation time than my final application can permit.

Can anyone suggest a more direct way to determine the "y" coordinate on the b-spline for a given "x"?

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

xp = [0., 0.71428571, 1.42857143, 2.14285714, 2.85714286, 3.57142857, 4.28571429, 5.]
yp = [0., -0.86217009, -2.4457478, -2.19354839, -2.32844575, -0.48680352, -0.41055718, -3.]

length = len(xp)
t = np.linspace(0., xp[-1], length - 2, endpoint=True)
t = np.append([0, 0, 0], t)
t = np.append(t, [xp[-1], xp[-1], xp[-1]])

tck = [t, [xp, yp], 3]
u = np.linspace(0, 5., 1000, endpoint=True)
out = interpolate.splev(u, tck)

x_value_in_xp_domain = 4.
y_value_out = interpolate.splev(x_value_in_xp_domain, tck)

plt.plot(xp, yp, linestyle='--', marker='o', color='purple')
plt.plot(out[0], out[1], color = 'teal')
plt.plot(x_value_in_xp_domain, y_value_out[1], marker='o', color = 'orangered')
plt.plot(y_value_out[0], y_value_out[1], marker='o', color = 'black')
plt.axvline(x=x_value_in_xp_domain, color = 'orangered')
plt.show()

The below image shows the guide polygon and the b-spline generated by the above code. The orange point at x=4 corresponds the point I wish to directly determine the y value of the b-spline at. The black point is the value of the b-spline when the value of 4 is passed as a parameter.

providing a few useful references:

Fast b-spline algorithm with numpy/scipy

https://github.com/kawache/Python-B-spline-examples

https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html

http://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node17.html

Answer 1

The roots method of the PPoly class may be the solution for your problem. This will be faster than Newton, and it will give you all the solutions in case there is more than one.

sx = interpolate.BSpline(t, xp, 3)
sy = interpolate.BSpline(t, yp, 3)

x0 = 4
u0 = interpolate.PPoly.from_spline((sx.t, sx.c - x0, 3)).roots()
sy(u0)

plt.plot(xp, yp, linestyle='--', marker='o', color='purple')
plt.plot(out[0], out[1], color = 'teal')
plt.plot(sx(u0), sy(u0), 'o')
plt.show()

Answer 2

Let me suggest a frame change. As you know, there are two ways to create a spline in the plane: as a function and as a parametric curve. You have chosen the latter but nowhere in the question do I see any compelling reason to do so; whereas taking the former would make a solution to your question almost trivial.

We start with some input data.

xdata  = [0.,  0.60626102,  1.36904762,  2.14285714,  2.85714286,  3.63095238,  4.39373898,  5.]
ydata  = [0., -0.82509685, -2.00782014, -2.25806452, -1.99902249, -0.77468231, -1.19975743, -3.]

Note that unlike your xp and yp, which are called control points in my circles, I collect a few points xdata, ydata that we will interpolate. I did that by sampling your function at the Greville abscissae to have something similar to your input. However, the exact type of input depends on your application; you have mentioned function approximation in one of the comments, so I guess this might be a way.

Now, we can use UnivariateSpline to create the spline function interpolating the data:

spline = interpolate.InterpolatedUnivariateSpline(xdata, ydata)

There are some alternatives; I recommend having a look at the interpolation tutorial and the references therein.

And this is more or less it. You can now sample the values at a lot of points to plot spline:

xfunc  = np.linspace(0, 5., 1000, endpoint=True)
yfunc  = spline(xfunc)
plt.plot(xfunc, yfunc,  color = 'green')

However, the reason we go through all this is to answer your question: how to find the y-coordinate corresponding to a given x-coordinate? This is quite simple:

x_value_in_xp_domain = 4.
y_value_out = spline(x_value_in_xp_domain)

Here is a picture showing my input data (blue; again, they are not identical to your inputs), the interpolating spline function (green), x_value_in_xp_domain (orange curve) and y_value_out (orange point).

Finally, here is the complete source code.

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

xdata  = [0.,  0.60626102,  1.36904762,  2.14285714,  2.85714286,  3.63095238,  4.39373898,  5.]
ydata  = [0., -0.82509685, -2.00782014, -2.25806452, -1.99902249, -0.77468231, -1.19975743, -3.]

spline = interpolate.InterpolatedUnivariateSpline(xdata, ydata)
xfunc  = np.linspace(0, 5., 1000, endpoint=True)
yfunc  = spline(xfunc)

x_value_in_xp_domain = 4.
y_value_out = spline(x_value_in_xp_domain)

plt.plot(xfunc, yfunc,  color = 'green')
plt.plot(xdata, ydata, 'o', color = 'blue')
plt.plot(x_value_in_xp_domain, y_value_out, marker='o', color = 'orangered')
plt.axvline(x=x_value_in_xp_domain, color = 'orangered')
plt.show()