使用 Python 和 Numpy 查找两个三次样条之间的最近点

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

我正在寻找一种方法来分析两个三次样条曲线并找到它们彼此最接近的点。我看过很多解决方案和帖子,但我一直无法实施建议的方法。我知道最近的点将是两条曲线的端点之一或两条曲线的一阶导数相等的点。检查端点很容易。很难找到一阶导数匹配的点。 给定:

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

最近点的候选点是:

B'(t) = C'(s)

每条曲线的一阶导数采用以下形式:

其中 a、b、c 系数由曲线的控制点形成:

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

为每条三次样条取 4 个控制点,我可以将每条曲线的参数部分转化为可以用 Numpy 表示的矩阵形式,代码如下 Python:

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)

这个 post 有一些很好的高级建议,但我不确定如何在 python 中实现一个解决方案,该解决方案可以找到具有匹配的一阶导数的 t 和 s 的正确值(连续下坡)。 B'(t) - C'(s) = 0 问题似乎是求根的问题。任何关于如何使用 python 和 Numpy 的建议将不胜感激。

使用 Numpy 假设问题应该用数字来解决。不失一般性,我们可以处理 0<s<=10<t<=1。您可以使用 SciPy 包来解决数字问题,例如

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)

这是两个圆(一个圆和一个圆弧)的示例。这可能适用于样条曲线。

您还可以使用函数 fmin

import numpy as np
import matplotlib.pylab as plt
from scipy.optimize import fmin

def BCubic(t, P0, P1, P2, P3):

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

    return a*3*(1-t)**2 + b*6*(1-t)*t + c*3*t**2

def B(t):
    return BCubic(t,4,2,3,1)

def C(t):
    return BCubic(t,1,4,3,4)

def f(t): 
    # L1 or manhattan distance
    return abs(B(t) - C(t))

init = 0 # 2
tmin = fmin(f,np.array([init]))
#Optimization terminated successfully.
#Current function value: 2.750000
#     Iterations: 23
#     Function evaluations: 46
print(tmin)
# [0.5833125]
tmin = tmin[0]

t = np.linspace(0, 2, 100)
plt.plot(t, B(t), label='B')
plt.plot(t, C(t), label='C')
plt.plot(t, abs(B(t)-C(t)), label='|B-C|')
plt.plot(tmin, B(tmin), 'r.', markersize=12, label='min')
plt.axvline(x=tmin, linestyle='--', color='k')
plt.legend()
plt.show()

您对 B'(t) = C'(s) 的假设太强了。

导数有方向和大小。方向必须与候选点重合,但幅度可能不同。

要找到具有相同导数斜率和最近距离的点可以求解方程组(当然,高幂:( )

 yb'(t) * xc'(u) - yc'(t) * xb'(u) = 0  //vector product of (anti)collinear vectors is zero
 ((xb(t) - xc(u))^2 + (xb(t) - xc(u))^2)' = 0   //distance derivative