4
\$\begingroup\$

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:

non-monotonic parametric values

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

my naive partitioning of the curve

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:

looks pretty good

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.