Querying points on a 3D spline at specific parametric values in Python
<p>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.</p>

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

<ol>
<li>subdivide spline into n subdivisions</li>
<li>add up subdivisions to compare agains queries</li>
<li>find the proper subdivision step for each query and extrapolate a position</li>
</ol>

<p>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.</p>

<p>Here's my function:</p>

<pre><code>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  
</code></pre>

<p>Here's a usage example:    </p>

<pre><code>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()
</code></pre>

<p>And the resulting plot:
<a href="https://i.stack.imgur.com/merun.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/merun.png" alt="enter image description here"></a></p>