# Find the closest parametric values corresponding to a BSpline's control points

On a curve generated by scipy.interpolate.BSpline I want to find the closest parameters relative to each control point, so that the given parametric range is monotonically increasing.

My first attempt was to naively sample the curve n times, find the index of the closest sample to each control point, and infer a parametric value from (closest index / number of samples) * max parameter.

import numpy as np
from numpy.core.umath_tests import inner1d
from scipy.interpolate import BSpline

def naive_subdivision(cv, degree=3, periodic=False, n=10000):
"""
Naively subdivide a BSpline and return the parametric
values of the closest point on curve to each control point.
"""
count = cv.shape[0]
max_param = count - (degree * (1-periodic))

# Closed curve
if periodic:
kv = np.arange(-degree,count+degree+1)
factor, fraction = divmod(count+degree+1, count)
cv = np.concatenate((cv,) * factor + (cv[:fraction],))
spl = BSpline(kv, cv, degree)

# Opened curve
else:
kv = np.clip(np.arange(count+degree+1)-degree,0,count-degree)
spl = BSpline(kv, cv, degree)

# Find the closest points
samples = spl(np.linspace(0,max_param,n)) # n samples on curve over parametric range
deltas = cv[:,np.newaxis]-samples # delta between control points and samples
mag = inner1d(deltas,deltas) # delta magnitudes, inner1d is faster than linalg.norm
ind = np.argmin(mag,axis=1) # closest index to each control point

# Return parameters expressed as a ratio of i*max_param/(n-1)
return ind * max_param / float(n-1)


Unfortunately this method didn't work as shows the animated gif below, where the parameters for points p1 and p3 jump under certain curve shapes:

# here are the 4 keyframes illustrated above, transposed to keep example formatting neat
keyframes = np.array([[[  3.,   3.,   3.,          10.54100958],
[ -3.,  -3.,  -2.15198582,  -2.15198582],
[ -2.,  -2.,  -2.,          -2.        ],
[  8.,   8.,   8.,           8.        ],
[  3.,   3.,   3.,           3.        ],
[-10.,  -7.,  -7.,          -7.        ]],
[[ -2.,  -2.,  -2.,          -1.2891395 ],
[ -4.,  -4.,  -1.33481258,  -1.33481258],
[  8.,   8.,   8.,           8.        ],
[  2.,   2.,   2.,           2.        ],
[-11., -11., -11.,         -11.        ],
[-11.,  -3.,  -3.,          -3.        ]]])
keyframes = keyframes.T

# keyframe[0]: all is well
print naive_subdivision(keyframes[0])
# [ 0., 0.33273327, 1.09810981, 1.88538854, 2.59135914, 3. ]
#
# keyframe[1]: eek!
print naive_subdivision(keyframes[1])
# [ 0., 2.88088809, 1.09810981, 1.88538854, 2.58475848, 3. ]
#       ^^^^^^^^^^
# keyframe[2]: back to normal
print naive_subdivision(keyframes[2])
# [ 0., 0.34713471, 1.04110411, 1.88508851, 2.58475848, 3. ]
#
# keyframe[3]: aak!
print naive_subdivision(keyframes[3])
# [ 0., 0.47824782, 1.04110411, 0.07770777, 2.58475848, 3. ]
#                               ^^^^^^^^^^


# My solution

To solve the problem I naively partition the BSpline into parametric segments for each control point, and apply the algorithm to each segment instead of the whole curve. The naive aspect of this approach is that it doesn't take into account the curve's actual shape, and just "hopes for the best" by keeping the subsamples within a reasonable range. The animated gif below shows each partition relative to it's control point:

Here is the full function:

import numpy as np
from numpy.core.umath_tests import inner1d
from scipy.interpolate import BSpline

def naive_partitioning(points, degree=3, periodic=False, n=100):
"""
Calculates the closest parameter on a BSpline curve relative
to each of the curve's control points.

points: BSpline curve's control points
degree: BSpline curve degree
periodic: True = curve is closed, False = curve is open
n: number of samples per partition
"""
points = np.asarray(points)
count = points.shape[0]

# On an open curve, clip degree so it doesn't exceed count-1
if not periodic:
degree = np.clip(degree,1,count-1)

# If the BSpline is linear (degree=1) stop right here
if degree <= 1:
return points, np.arange(count)

# Calculate BSpline's maximum parameter
max_param = count - (degree * (1-periodic))

#---BSpline function---#

# Periodic curve
if periodic:
# Create a periodic knot vector, ex: [-2, -1, 0, 1, 2, 3, 4, 5, 6]
kv = np.arange(-degree,count+degree+1)

# Create the BSpline function with "wrap around" control points
factor, fraction = divmod(count+degree+1, count)
spl = BSpline(kv, np.concatenate((points,) * factor + (points[:fraction],)), degree)

# Append the first control point to the end of the sequence to create a closed loop
points = np.concatenate((points,[points[0]]))
count = points.shape[0] # reset point count

# Opened curve
else:
# Create BSpline function with an open knot vector, ex [0, 0, 0, 0, 1, 2, 3, 3, 3, 3]
kv = np.clip(np.arange(count+degree+1)-degree,0,count-degree)
spl = BSpline(kv, points, degree)

#---Naive curve sub-sample partitioning---#
#
# For each control point, we'll look for a closest sample that lies inside a sample partition.
# For example, this will force the algorithm to find a closest point for p[1] somewhere between p[0] and p[2].
# This is to force the algotithm to produce an increasingly monotonic parameter sequence.
#
# Curve parameter for an opened curve will have always p[0]=0 and p[-1]=max_param
# For a periodic curve we offset the parametric range by "0.5*(degree+1)-1" to force the
# algorithm to keep p[0]=0 and p[-1]=max_param under an curve degree

# Periodic curve
if periodic:
# Calculate parametric offset according to curve degree.
# This will register parametric queries relative to the first control point.
offset =  0.5 * (degree+1) - 1

# Create a periodic parametric range. ex: [-2., -0.875, 0.25, 1.375, 2.5, 3.625, 4.75, 5.875, 7.]
u = np.linspace(-1,max_param+1,count+2) - offset

# Create partition matrix of n parametric subsamples per partition
# ex: p[0] will be searching between linspace(-2,0.25,n), p[1] between linspace(-0.875,1.375,n), etc...
_min = u[np.arange(0,u.shape[0]-2)]
_max = u[np.arange(2,u.shape[0])]
u = (_min[:, np.newaxis] + (_max-_min)[:, np.newaxis]/(n-1) * np.arange(n))%max_param

# Opened curve
else:
# Create an open parameric range. ex: [ 0. ,  0. ,  0.6,  1.2,  1.8,  2.4,  3. ,  3. ]
u = np.linspace(0,max_param,count)
u = np.insert(u,0,0)
u = np.insert(u,-1,u[-1])

# Create partition matrix of n parametric subsamples per partition
# ex: p[0] will be searching between linspace(0,0.6,n), p[1] between linspace(0,1.2,n), etc...
_min = u[np.arange(0,u.shape[0]-2)]
_max = u[np.arange(2,u.shape[0])]
u = _min[:, np.newaxis] + (_max-_min)[:, np.newaxis]/(n-1) * np.arange(n)

# Create samples for each partition
arange = np.arange(count)
samples = spl(u)

# Get the distances between control points to sample partition
deltas = points[arange,np.newaxis]-samples
mag = inner1d(deltas,deltas) # inner1d is faster than np.linalg.norm

# Get the index of the closest point per partition
index = np.argmin(mag,axis=1)

# Reverse modulus to maintain monotonicity on periodic curves
u = u[arange,index]
if periodic:
diff = np.diff(u)
if diff.min() < 0:
i = np.arange(0,np.where(diff < 0)[0][0] + 1)
u[i] = 0 - (count-1 - u[i])

# Return the closest points and corresponding parametric values
return samples[arange,index], u


This approach seems to work in both open and periodic forms:

# Same keyframes, better results:
print naive_partitioning(keyframes[0])[1]
# [ 0., 0.32727273, 1.0969697 , 1.89090909, 2.58787879, 3. ]
print naive_partitioning(keyframes[1])[1]
# [ 0., 0.32727273, 1.0969697 , 1.89090909, 2.58787879, 3. ]
print naive_partitioning(keyframes[2])[1]
# [ 0., 0.35151515, 1.03636364, 1.89090909, 2.58787879, 3. ]
print naive_partitioning(keyframes[3])[1]
# [ 0., 0.47272727, 1.03636364, 1.89090909, 2.58787879, 3. ]


I would love to know of better methods to get the parametric values I'm after. And should my approximation method be "the way to do it" I would welcome all suggestions to improve the function's execution speed and reliability.

• Could you supply some test cases? Perhaps the four key curves of the first gif? – Peter Taylor Oct 10 '17 at 10:17
• @PeterTaylor i added the original keyframe data, as well as the code to show how the original attempt failed. Hope this helps! – Fnord Oct 10 '17 at 22:46