numpy 巨蟒:距离第三点最近的直线上的点

lsmd5eda  于 2022-11-10  发布在  其他
关注(0)|答案(3)|浏览(156)

我在两个XY点(p1和p2)和线外的第三个XY点(P3)之间有一条直线/向量。根据this post,我知道如何得到该点到直线的距离。但我实际上要找的是直线上的一个点(P4),它离第三个点(P3)有一个最小的距离(D)。我找到了this post,但我觉得这不是正确的解决方案。也许在Numpy或Python中包含了一些东西?

根据@allo的说法,我尝试了以下方法。您可以将我的代码下载为Python fileJupyter 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)**。

bfnvny8b

bfnvny8b1#

解析几何

让我们从指定的直线开始,我们根据其上的两个点(x1, y1)(x2, y2)来定义直线。
对于dx = x2-x1dy = 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)
分别写出xy的方程式

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

它是ab中的线性系统,其解为

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,所以只有当两个初始点相对于它们到第三个点的w/r距离非常接近时,才会有问题。

Python实现

我们使用浮点数的二元组来表示二维点,并定义了一个函数,其参数是三个二元组,表示定义直线的点(p1, p2)和确定p4在该直线上的位置的点(p3)。

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的预期一致。

Matplotlib示例

import matplotlib.pyplot as plt

def p(p1, p2, p3):
    (x1, y1), (x2, y2), (x3, y3) = p1, p2, 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

p1, p2, p3 = (2, 4), (7, 3), (1, 1)
p4 = p(p1, p2, p3)

fig, ax = plt.subplots()

# if we are after  right angles, anything else would be wrong

ax.set_aspect(1)

plt.plot(*zip(p1, p2, p4, p3), marker='*')
z31licg0

z31licg02#

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)'

wsxa1bj1

wsxa1bj13#

您要做的是vector projection
p1p3旋转到边p1p2上,您需要找到线段p1p4的正确长度。然后您可以只使用p1+FACTOR*p1p2 / length(p1p2)。所需的系数由p1p2p1p3之间的夹角的余弦给出。然后你就会得到

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

下面以两种边缘情况为例:

  • 如果p1p3p1p2正交,则余弦为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))

相关问题