计算点到轴的距离

Question

我有一个在三维空间中的点数组：

P = np.random.random((10, 3))

现在我想找出这些点到某个特定轴的距离，以及在那个轴上的距离。

Ax_support = array([3, 2, 1])
Ax_direction = array([1, 2, 3])

我找到了一种解决方案，首先计算每个点到方向向量的垂直向量……不过我觉得这个方法有点复杂，对于这样一个常见的问题，应该已经有现成的numpy或scipy的函数可以用（就像在scipy.spatial.distance中找到点之间的距离一样）。

Answer 1

我一直觉得这个功能很重要，但市面上没有，于是我在我维护的工具库中添加了一个叫 haggis.math.segment_distance 的功能，库的名字是 haggis。

不过有一点需要注意，这个函数是通过两个点来定义一条线，而不是通过一个支点和方向。你可以输入任意数量的点和线，只要它们的维度能够匹配。通常的情况是很多点投影到一条线上（就像你那样），或者一个点投影到多条线上。

distances = haggis.math.segment_distance(
        P, Ax_support, Ax_support + Ax_direction,
        axis=1, segment=False)

这里有一个可以复现的例子：

np.random.seed(0xBEEF)
P = np.random.random((10,3))
Ax_support = np.array([3, 2, 1])
Ax_direction = np.array([1, 2, 3])
d, t = haggis.math.segment_distance(
        P, Ax_support, Ax_support + Ax_direction,
        axis=1, return_t=True, segment=False)

return_t 会返回沿着这条线的法线点的位置，表示为从支点开始的线段长度的比例（也就是说，Ax_support + t * Ax_direction 是投影点的位置）。

>>> d
array([2.08730838, 2.73314321, 2.1075711 , 2.5672012 , 1.96132443,
       2.53325436, 2.15278454, 2.77763701, 2.50545181, 2.75187883])

>>> t
array([-0.47585462, -0.60843258, -0.46755277, -0.4273361 , -0.53393468,
       -0.58564737, -0.38732655, -0.53212317, -0.54956548, -0.41748691])

这样你就可以绘制出以下的图：

fig, ax = plt.subplots(subplot_kw={'projection': '3d'})

ax.plot(*Ax_support, 'ko')
ax.plot(*Ax_support + Ax_direction, 'ko')

start = min(t.min(), 0)
end = max(t.max(), 1)
margin = 0.1 * (end - start)
start -= margin
end += margin

ax.plot(*Ax_support[:, None] + [start, end] * Ax_direction[:, None], 'k--')
for s, te in zip(P, t):
    pt = Ax_support + te * Ax_direction
    ax.plot(*np.stack((s, pt), axis=-1), 'k-')
    ax.plot(*pt, 'ko', markerfacecolor='w')
ax.plot(*P.T, 'ko', markerfacecolor='r')
plt.show()

Answer 2

我会很惊讶在numpy/scipy的标准操作中看到这样的操作。你想要的其实是投影到你那条线上的距离。首先要做的是减去Ax_support：

P_centered = P - Ax_support

通过0点并且方向为Ax_direction的那条线，离每个P_centered的点最近的点可以用下面的公式表示：

P_projected = P_centered.dot(np.linalg.pinv(
             Ax_direction[np.newaxis, :])).dot(Ax_direction[np.newaxis, :])

所以你要找的公式就是：

distances = np.sqrt(((P_centered - P_projected) ** 2).sum(axis=1))

没错，这正是你提到的，以一种向量化的方式来处理，所以对于合理数量的数据点来说，这应该会很快。

注意：如果有人知道有内置的函数可以做到这一点，我会非常感兴趣！