Python：点在最接近第三点的线上

Question

我在两个 XY 点 (p1 和 p2) 和第三个 XY 点 (p3) 之间有一条线/向量，它位于该线之外。根据this post，我知道如何获得该点到线的距离。但我实际上正在寻找的是该线上的一个点（p4），它与第三点（p3）的最小距离（d）。我找到了this post，但我觉得这不是正确的解决方案。也许 Numpy 或 Python 中包含一些东西？

根据@allo，我尝试了以下方法。您可以将我的代码下载为Python file 或Jupyter Notebook（均为Python3）。

points = [[1, 1], [3, 1], [2.5, 2], [2.5, 1]]
import matplotlib.pyplot as plt
%matplotlib inline

fig, ax = plt.subplots()
fig.set_size_inches(6,6)

x, y = zip(*points[:2])
l1, = ax.plot(x,y, color='blue')
scatter1 = ax.scatter(x=x,y=y, color='blue', marker='x', s=80, alpha=1.0)

x, y = zip(*points[2:])
l2, = ax.plot(x,y, color='red')
scatter2 = ax.scatter(x=x,y=y, color='red', marker='x', s=80, alpha=1.0)

p1 = Vector2D(*points[0])
p2 = Vector2D(*points[1])
p3 = Vector2D(*points[2])

p1p2 = p2.sub_vector(p1)
p1p3 = p3.sub_vector(p1)

angle_p1p2_p1p3 = p1p2.get_angle_radians(p1p3)
length_p1p3 = p1p3.get_length()
length_p1p2 = p1p2.get_length()

p4 = p1.add_vector(p1p2.multiply(p1p3.get_length()/p1p2.get_length()).multiply(math.cos(p1p2.get_angle_radians(p1p3))))

#p4 = p1 + p1p2 * length(p1p3)/length(p1p2)*cos(angle(p1p2, p1p3))

p4 = p1.add_vector(p1p2.multiply(length_p1p3/length_p1p2*math.cos(angle_p1p2_p1p3)))
p4

这导致 p4 = (1.8062257748298551, 1.0) 但显然应该是 (2.5, 1.0)。

Answer 1

解析几何

让我们从指定的线开始，我们根据上面的两个点来定义线(x1, y1) 和(x2, y2)。

使用dx = x2-x1 和dy = y2-y1，我们可以将直线上的每个点正式写为(x12, y12) = (x1, y1) + a*(dx, dy)，其中a 是一个实数。

使用类似的符号，通过(x3, y3) 并垂直于指定点的线上的点是(x34, y34) = (x3, y3) + b*(-dy, +dx)。

要找到交叉点，我们必须施加 (x12, y12) = (x34, y34) 或 (x1, y1) + a*(dx, dy) = (x3, y3) + b*(-dy, +dx).

分别编写x 和y 的方程式

y1 + a dy - y3 - b dx = 0
x1 + a dx + b dy - x3 = 0

它是a和b中的一个线性系统，其解是

a = (dy y3 - dy y1 + dx x3 - dx x1) / (dy^2 + dx^2)
b = (dy x3 - dy x1 - dx y3 + dx y1) / (dy^2 + dx^2)

离(x3, y3)最近的点的坐标位于线上 是(x1+a*dx, y1+a*dy)——你只需要计算系数a。

从数值上讲，线性系统的行列式是dx**2+dy**2，因此只有当两个初始点到第三个点的距离非常接近时，才会出现问题。

Python 实现

我们使用一个 2-uple 浮点数来表示一个 2D 点，我们定义一个函数，其参数是 3 个 2-uple 表示定义线的点 (p1, p2) 和确定的点 (p3) p4 在该行上的位置。

In [16]: def p4(p1, p2, p3):
    ...:     x1, y1 = p1
    ...:     x2, y2 = p2
    ...:     x3, y3 = p3
    ...:     dx, dy = x2-x1, y2-y1
    ...:     det = dx*dx + dy*dy
    ...:     a = (dy*(y3-y1)+dx*(x3-x1))/det
    ...:     return x1+a*dx, y1+a*dy

为了测试实现，我使用了 OP 使用的三个点 展示他们对这个问题的看法：

In [17]: p4((1.0, 1.0), (3.0, 1.0), (2.5, 2))
Out[17]: (2.5, 1.0)

看来p4(...)的结果和OP的预期不谋而合。

Answer 2

Shapely 的distance() 函数返回最小距离：

>>> from shapely.geometry import LineString as shLs
>>> from shapely.geometry import Point as shPt
>>> l = shLs([ (1,1), (3,1)])
>>> p = shPt(2,2)
>>> dist = p.distance(l)
1.0
>>> l.interpolate(dist).wkt
'POINT (2 1)'

Answer 3

你想做的是vector projection。

边p1p3 旋转到边p1p2 上，您需要找到段p1p4 的正确长度。然后你可以使用p1+FACTOR*p1p2 / length(p1p2)。所需因子由p1p2 和p1p3 之间的夹角余弦给出。然后你得到

p4 = p1 + p1p2 * length(p1p3)/length(p1p2)*cos(angle(p1p2, p1p3))

这里以两种极端情况为例：

• 如果p1p3 与p1p2 正交，则余弦为0，因此p4 位于p1 上。
• 当p1p3 位于p1p2 上时，余弦为1，所以p1p2 只是通过length(p1p3)/length(p1p2) 缩放得到p1p4。

您还可以将余弦替换为点积dot(p1p2 / length(p1p2), p1p3 / length(p1p3)。

您可以在in the wikibook about linear algebra找到更多细节和精美插图。

这里有一个从您的 python 代码派生的完整示例。我在这里使用 numpy 而不是 Vector2D。

points = [[1, 1], [3, 1], [2.5, 2], [2.5, 1]]
import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np

fig, ax = plt.subplots()
fig.set_size_inches(6,6)

x, y = zip(*points[:2])
l1, = ax.plot(x,y, color='blue')
scatter1 = ax.scatter(x=x,y=y, color='blue', marker='x', s=80, alpha=1.0)

x, y = zip(*points[2:])
l2, = ax.plot(x,y, color='red')
scatter2 = ax.scatter(x=x,y=y, color='red', marker='x', s=80, alpha=1.0)

p1 = np.array(points[0])
p2 = np.array(points[1])
p3 = np.array(points[2])

p1p2 = p2 - p1
p1p3 = p3 - p1

p4 = p1 + p1p2 / np.linalg.norm(p1p2) * np.linalg.norm(p1p3) * ((p1p2/np.linalg.norm(p1p2)).T * (p1p3/np.linalg.norm(p1p3)))

p1, p2, p3, p4, p1p2, p1p3

我们可以像这样缩短p4 行，使用标量积的线性：

p4 = p1 + p1p2 / np.linalg.norm(p1p2) * ((p1p2/np.linalg.norm(p1p2)).T * (p1p3))
p4 = p1 + p1p2 / np.linalg.norm(p1p2)**2 * (p1p2.T * (p1p3))