使用 Numpy 找到两个线性方程的交点

Find the point of intersection of two linear equations using Numpy

objective就是求两个一次方程的交点。这两个线性方程是使用 Numpy polyfit 函数导出的。

给定两个时间序列 (xLeftyLeft) 和 (xRightyRight),使用 polyfit如下图:

xLeft = [
    6168, 6169, 6170, 6171, 6172, 6173, 6174, 6175, 6176, 6177,
    6178, 6179, 6180, 6181, 6182, 6183, 6184, 6185, 6186, 6187
]
yLeft = [
    0.98288751, 1.3639959, 1.7550986, 2.1539073, 2.5580614,
    2.9651523, 3.3727503, 3.7784295, 4.1797948, 4.5745049,
    4.9602985, 5.3350167, 5.6966233, 6.0432272, 6.3730989,
    6.6846867, 6.9766307, 7.2477727, 7.4971657, 7.7240791
]

xRight = [
    6210, 6211, 6212, 6213, 6214, 6215, 6216, 6217, 6218, 6219,
    6220, 6221, 6222, 6223, 6224, 6225, 6226, 6227, 6228, 6229,
    6230, 6231, 6232, 6233, 6234, 6235, 6236, 6237, 6238, 6239,
    6240, 6241, 6242, 6243, 6244, 6245, 6246, 6247, 6248, 6249,
    6250, 6251, 6252, 6253, 6254, 6255, 6256, 6257, 6258, 6259,
    6260, 6261, 6262, 6263, 6264, 6265, 6266, 6267, 6268, 6269,
    6270, 6271, 6272, 6273, 6274, 6275, 6276, 6277, 6278, 6279,
    6280, 6281, 6282, 6283, 6284, 6285, 6286, 6287, 6288]
yRight=[
    7.8625913, 7.7713094, 7.6833806, 7.5997391, 7.5211883,
    7.4483986, 7.3819046, 7.3221073, 7.2692747, 7.223547,
    7.1849418, 7.1533613, 7.1286001, 7.1103559,  7.0982385,
    7.0917811, 7.0904517, 7.0936642, 7.100791, 7.1111741,
    7.124136, 7.1389918, 7.1550579, 7.1716633, 7.1881566,
    7.2039142, 7.218349, 7.2309117, 7.2410989, 7.248455,
    7.2525721, 7.2530937, 7.249711, 7.2421637, 7.2302341,
    7.213747, 7.1925621, 7.1665707, 7.1356878, 7.0998487,
    7.0590014, 7.0131001, 6.9621005, 6.9059525, 6.8445964,
    6.7779589, 6.7059474, 6.6284504, 6.5453324, 6.4564347,
    6.3615761, 6.2605534, 6.1531439, 6.0391097, 5.9182019,
    5.7901659, 5.6547484, 5.5117044, 5.360805, 5.2018456,
    5.034656, 4.8591075, 4.6751242, 4.4826899, 4.281858,
    4.0727611, 3.8556159, 3.6307325, 3.3985188, 3.1594861,
    2.9142516, 2.6635408, 2.4081881, 2.1491354, 1.8874279,
    1.6242117,1.3607255,1.0982931,0.83831298
]

left_line = np.polyfit(xleft, yleft, 1)
right_line = np.polyfit(xRight, yRight, 1)

在这种情况下,polyfit 分别输出 y = mx + b 的系数 mb

然后两个线性方程的交集可以计算如下:

x0 = -(left_line[1] - right_line[1]) / (left_line[0] - right_line[0])
y0 = x0 * left_line[0] + left_line[1]

但是,我想知道是否存在 Numpy 内置方法来计算最后两个步骤?

您可以将其建模为线性系统并使用简单的线性代数:

def get_intersection(m1,b1,m2,b2):
    A = np.array([[-m1, 1], [-m2, 1]])
    b = np.array([[b1], [b2]])
    # you have to solve linear System AX = b where X = [x y]'
    X = np.linalg.pinv(A) @ b
    x, y = np.round(np.squeeze(X), 4)
    return x, y # returns point of intersection (x,y) with 4 decimal precision

m1,b1,m2,b2 = left_line(0), left_line(1), right_line(0), right_line(1)
print(get_intersection(m1,b1,m2,b2))

例如,对于行 y - x = 1 和 y + x = 1,我们期望交点为 (0,1):

m1,b1,m2,b2 = 1, 1, -1, 1
print(get_intersection(m1,b1,m2,b2))

输出:(0.0, 1.0) 符合预期。

不完全是内置方法,但您可以简化问题。假设我有行给定 y = m1 * x + b1y = m2 * x + b2。您可以轻松找到差值的方程式,这也是一条线:

y = (m1 - m2) * x + (b1 - b2)

注意这条线将在两条原始线的交点处有一个根,如果它们相交的话。您可以使用 numpy.polynomial.Polynomial class 来执行这些操作:

>>> (np.polynomial.Polynomial(left_line[::-1]) - np.polynomial.Polynomial(right_line[::-1])).roots()
array([6192.0710885])

请注意,我必须交换系数的顺序,因为 Polynomial 期望从最小到最大,而 np.polyfit returns the opposite. In fact, np.polyfit is not recommended. Instead, you can get Polynomial obejcts directly using np.polynomial.Polynomial.fit class 方法。您的代码将如下所示:

left_line = np.polynomial.Polynomial.fit(xLeft, yLeft, 1, domain=[-1, 1])
right_line = np.polynomial.Polynomial.fit(xRight, yRight, 1, domain=[-1, 1])
x0 = (left_line - right_line).roots()
y0 = left_line(x0)

domain 映射到 window [-1, 1]。如果您未指定 domain,则将使用 x 值的峰峰值。您不希望这样,因为它会导致输入值的映射。相反,我们明确指定域 [-1, 1] 映射到相同的 window。另一种方法是使用默认域并设置例如window=[xLeft.min(), xLeft.max()]。这种方法的问题是它会为多项式创建不同的域,从而阻止操作 left_line - right_line.

有关详细信息,请参阅 https://numpy.org/doc/stable/reference/routines.polynomials.classes.html