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

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 against 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 built-in 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 which 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,y))

# 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:

Continue reading...