Querying points on a 3D spline at specific parametric values in Python

Question

Given a list of control vertices defining a spline, and list of query values (0 = line start, 1= line end, 0.5= half way, 0.25= quarter of the way, etc...) I want to find (as fast and efficiently as possible) the 3D coordinates of those queries on the spline.

I tried to find something built-in scipy but failed. So i wrote the function to solve the problem with a brute force approach:

1. subdivide spline into n subdivisions
2. add up subdivisions to compare agains queries
3. find the proper subdivision step for each query and extrapolate a position

The code below works fine, but i'd love to know if there are faster/more efficient ways of calculating what i need, or better yet something already builtin scipy I may have missed.

Here's my function:

import numpy as np
import scipy.interpolate as interpolate

def uQuery(cv,u,steps=100,projection=True):
    ''' Brute force point query on spline
        cv     = list of spline control vertices
        u      = list of queries (0-1)
        steps  = number of curve subdivisions (higher value = more precise result)
        projection = method by wich we get the final result
                     - True : project a query onto closest spline segments.
                              this gives good results but requires a high step count
                     - False: modulates the parametric samples and recomputes new curve with splev.
                              this can give better results with fewer samples.
                              definitely works better (and cheaper) when dealing with b-splines (not in this examples)

    '''
    u = np.clip(u,0,1) # Clip u queries between 0 and 1

    # Create spline points
    samples = np.linspace(0,1,steps)
    tck,u_=interpolate.splprep(cv.T,s=0.0)
    p = np.array(interpolate.splev(samples,tck)).T  
    # at first i thought that passing my query list to splev instead
    # of np.linspace would do the trick, but apparently not.    

    # Approximate spline length by adding all the segments
    p_= np.diff(p,axis=0) # get distances between segments
    m = np.sqrt((p_*p_).sum(axis=1)) # segment magnitudes
    s = np.cumsum(m) # cumulative summation of magnitudes
    s/=s[-1] # normalize distances using its total length

    # Find closest index boundaries
    s = np.insert(s,0,0) # prepend with 0 for proper index matching
    i0 = (s.searchsorted(u,side='left')-1).clip(min=0) # Find closest lowest boundary position
    i1 = i0+1 # upper boundary will be the next up

    # Return projection on segments for each query
    if projection:
        return ((p[i1]-p[i0])*((u-s[i0])/(s[i1]-s[i0]))[:,None])+p[i0]

    # Else, modulate parametric samples and and pass back to splev
    mod = (((u-s[i0])/(s[i1]-s[i0]))/steps)+samples[i0]
    return np.array(interpolate.splev(mod,tck)).T  

Here's a usage example:

import matplotlib.pyplot as plt

cv = np.array([[ 50.,  25.,  0.],
   [ 59.,  12.,  0.],
   [ 50.,  10.,   0.],
   [ 57.,   2.,   0.],
   [ 40.,   4.,   0.],
   [ 40.,   14.,  0.]])


# Lets plot a few queries
u = [0.,0.2,0.3,0.5,1.0]
steps = 10000 # The more subdivisions the better
x,y,z = uQuery(cv,u,steps).T
fig, ax = plt.subplots()
ax.plot(x, y, 'bo')
for i, txt in enumerate(u):
    ax.annotate('  u=%s'%txt, (x[i],y[i]))

# Plot the curve we're sampling
tck,u_=interpolate.splprep(cv.T,s=0.0)
x,y,z = np.array(interpolate.splev(np.linspace(0,1,1000),tck))
plt.plot(x,y,'k-',label='Curve')

# Plot control points
p = cv.T
plt.scatter(p[0],p[1],s=80, facecolors='none', edgecolors='r',label='Control Points')

plt.minorticks_on()
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.xlim(35, 70)
plt.ylim(0, 30)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

And the resulting plot:

Answer 1

Sorry for my previous comment, I had misunderstood the question.

Mind you, I will be referring to your query as w = [0.,0.2,0.3,0.5,1.0], because I use u for something else.

Sadly, there is not a simple solution for your problem, since it implies calculating the length of a cubic spline, and this is not trivial. But there is a way to simplify your code using the integrate and optimize scipy libraries so you don't have to worry about precision.

First you have to understand that under the hood, splprep creates a cubic spline of the shape x=Fx(u) and y=Fy(u), where u is a parameter from 0 to 1, but is not linear related to the length of the spline, for example for this control points:

cv = np.array([[ 0.,  0.,  0.],
   [ 100,  25,   0.],
   [ 0.,  50.,   0.],
   [ 100,  75,   0.],
   [ 0.,   100.,  0.]])

You can see how the parameter u behaves. It is worth noticing that you can define which value of u you want for your control points, and this would have an impact on the shape of the spline.

Now, when you call splev, you are really asking for the coordinates of the spline for a given u parameter. So in order to do what you want to do, you need to find the u for a given fraction of the length of the spline.

First, in order to get the total length of the spline, there is not much you can do, but a numerical integration, like you did, but you can use the integrate library of scipy to do it easier.

import scipy.integrate as integrate

def foo(u):
     xx,yy,zz=interpolate.splev(u,tck,der=1)
     return (xx**2 + yy**2)**0.5

total_length=integrate.quad(foo,0,1)[0]

Once you have the total length of the spline, you can use the optimize library to find the value of u that would integrate to the fraction of the length that you want. And using this desired_u into splev will give you the coordinates you want.

import scipy.optimize as optimize

desired_u=optimize.fsolve(lambda uu:  integrate.quad(foo,0,uu)[0]-w*total_length,0)[0]

x,y,z = np.array(interpolate.splev(desired_u,tck))

---

Edit: I measured the performance of my method vs yours, and yours is way faster, and also very precise, the only metric in which is worst is in memory allocation. I found a way to speed up my method, still with low memory allocation, but it sacrifice precision.

I'll be using 100 query points as test.

My method as it is right now:

time = 21.686723 s
memory allocation = 122.880 kb

I'll be using the points given by my method as the true coordinates, and measure the average distance of the 100 points between each method and these ones

Your method as it is right now:

time = 0.008699 s
memory allocation = 1,187.840 kb
Average distance = 1.74857994144e-06

It is possible to improve the speed of my method by not using fsolve over the integral to each point, but by creating an interpolating function u=F(w) from a sample of points, and then use that function, which would be much faster.

import scipy.interpolate as interpolate
import scipy.integrate as integrate


def foo(u):
     xx,yy,zz=interpolate.splev(u,tck,der=1)
     return (xx**2 + yy**2)**0.5

total_length=integrate.quad(foo,0,1)[0]

yu=[integrate.quad(foo,0,uu)[0]/total for uu in np.linspace(0,1,50)]

find_u=interpolate.interp1d(yu,np.linspace(0,1,50))

x,y,z=interpolate.splev(find_u(w),tck)

With 50 samples I get:

time = 1.280629 s
memory allocation = 20.480 kb
Average distance = 0.226036973904

Which is much faster than before, but still not as fast as yours, precision also not as good as yours, but much better in terms of memory. But it depends as yours do, in the number of samples.

Your method with 1000, and 100 points:

1000 points
time = 0.002354 s
memory allocation = 167.936 kb
Average distance = 0.000176413655938

100 points
time = 0.001641 s
memory allocation = 61.440 kb
Average distance = 0.0179918600812

My method with 20 and 100 samples

20 samples
time = 0.514241 s
memory allocation = 14.384 kb
Average distance = 1.42356341648

100 samples
time = 2.45364 s
memory allocation = 24.576 kb
Average distance = 0.0506075927139

All things considered, I think your method is better with the right number of points for the precision needed, mine is only fewer lines of code.

Edit 2: I just realize something else, your method can give points outside of the spline, while mine is always in the spline, depending on what you are doing this might be important