Search code examples
pythonnumpymathbezierspline

Finding the Closest Points between Two Cubic Splines with Python and Numpy

I'm looking for a way to analyze two cubic splines and find the point where they come the closest to each other. I've seen a lot of solutions and posts but I've been unable to implement the methods suggested. I know that the closest point will be one of the end-points of the two curves or a point where the first derivative of both curves is equal. Checking the end points is easy. Finding the points where the first derivatives match is hard. Given:

Curve 0 is B(t)   (red)
Curve 1 is C(s)   (blue)

enter image description here

A candidate for closest point is where:

B'(t) = C'(s)

The first derivative of each curve takes the following form: enter image description here

Where the a, b, c coefficients are formed from the control points of the curves:

a=P1-P0
b=P2-P1
c=P3-P2

Taking the 4 control points for each cubic spline I can get each curve's parametric sections into a matrix form that can be expressed with Numpy with the following Python code:

def test_closest_points():
    # Control Points for the two qubic splines.
    spline_0 = [(1,28), (58,93), (113,95), (239,32)]
    spline_1 = [(58, 241), (26,76), (225,83), (211,205)]

    first_derivative_matrix = np.array([[3, -6, 3], [-6, 6, 0], [3, 0, 0]])

    spline_0_x_A = spline_0[1][0] - spline_0[0][0]
    spline_0_x_B = spline_0[2][0] - spline_0[1][0]
    spline_0_x_C = spline_0[3][0] - spline_0[2][0]

    spline_0_y_A = spline_0[1][1] - spline_0[0][1]
    spline_0_y_B = spline_0[2][1] - spline_0[1][1]
    spline_0_y_C = spline_0[3][1] - spline_0[2][1]

    spline_1_x_A = spline_1[1][0] - spline_1[0][0]
    spline_1_x_B = spline_1[2][0] - spline_1[1][0]
    spline_1_x_C = spline_1[3][0] - spline_1[2][0]

    spline_1_y_A = spline_1[1][1] - spline_1[0][1]
    spline_1_y_B = spline_1[2][1] - spline_1[1][1]
    spline_1_y_C = spline_1[3][1] - spline_1[2][1]

    spline_0_first_derivative_x_coefficients = np.array([[spline_0_x_A], [spline_0_x_B], [spline_0_x_C]])
    spline_0_first_derivative_y_coefficients = np.array([[spline_0_y_A], [spline_0_y_B], [spline_0_y_C]])

    spline_1_first_derivative_x_coefficients = np.array([[spline_1_x_A], [spline_1_x_B], [spline_1_x_C]])
    spline_1_first_derivative_y_coefficients = np.array([[spline_1_y_A], [spline_1_y_B], [spline_1_y_C]])

    # Show All te matrix values
    print 'first_derivative_matrix:'
    print first_derivative_matrix
    print
    print 'spline_0_first_derivative_x_coefficients:'
    print spline_0_first_derivative_x_coefficients
    print
    print 'spline_0_first_derivative_y_coefficients:'
    print spline_0_first_derivative_y_coefficients
    print
    print 'spline_1_first_derivative_x_coefficients:'
    print spline_1_first_derivative_x_coefficients
    print
    print 'spline_1_first_derivative_y_coefficients:'
    print spline_1_first_derivative_y_coefficients
    print

# Now taking B(t) as spline_0 and C(s) as spline_1, I need to find the values of t and s where B'(t) = C'(s)

This post has some good high-level advice but I'm unsure how to implement a solution in python that can find the correct values for t and s that have matching first derivatives (slopes). The B'(t) - C'(s) = 0 problem seems like a matter of finding roots. Any advice on how to do it with python and Numpy would be greatly appreciated.

Solution

  • Using Numpy assumes that the problem should be solved numerically. Without loss of generality we can treat that 0<s<=1 and 0<t<=1. You can use SciPy package to solve the problem numerically, e.g.

    from scipy.optimize import minimize
import numpy as np

def B(t):
    """Assumed for simplicity: 0 < t <= 1
    """
    return np.sin(6.28 * t), np.cos(6.28 * t)

def C(s):
    """0 < s <= 1
    """
    return 10 + np.sin(3.14 * s), 10 + np.cos(3.14 * s)



def Q(x):
    """Distance function to be minimized
    """
    b = B(x[0])
    c = C(x[1])
    return (b[0] - c[0]) ** 2 + (b[1] - c[1]) ** 2

res = minimize(Q, (0.5, 0.5))


print("B-Point: ", B(res.x[0]))
print("C-Point: ", C(res.x[1]))

    B-Point: (0.7071067518175205, 0.7071068105555733)
    C-Point: (9.292893243165555, 9.29289319446135)

    This is example for two circles (one circle and arc). This will likely work with splines.