Finding the distance of points to an axis

Question

I have an array of points in 3D Cartesian space:

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

Now I'd like to find their distances to a given axis and on that given axis

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

I've found a solution that first finds the vector from each point that is perpendicular to the direction vector... I feel however, that it is really complicated and that for such a standard problem there would be a numpy or scipy routine already out there (as there is to find the distance between points in scipy.spatial.distance)

Answer 1

I would be surprised to see such an operation among the standard operations of numpy/scipy. What you are looking for is the projection distance onto your line. Start by subtracting Ax_support:

P_centered = P - Ax_support

The points on the line through 0 with direction Ax_direction with the shortest distance in the L2 sense to each element of P_centered is given by

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

Thus the formula you are looking for is

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

Yes, this is exactly what you propose, in a vectorized way of doing things, so it should be pretty fast for reasonably many data points.

Note: If anybody knows a built-in function for this, I'd be very interested!

Answer 2

It's been bothering me that something like this does not exist, so I added haggis.math.segment_distance to the library of utilities I maintain, haggis.

One catch is that this function expects a line defined by two points rather than a support and direction. You can supply as many points and lines as you want, as long as the dimensions broadcast. The usual formats are many points projected on one line (as you have), or one point projected on many lines.

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

Here is a reproducible example:

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 returns the location of the normal points along the line as a ratio of the line length from the support (i.e., Ax_support + t * Ax_direction is the location of the projected points).

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

This allows you to make the following plot:

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