0

Inspired by the answer here, I would like to calculate the perpendicular distance (in vector format instead of just magnitude) from a point to a straight line in 3D space.

The above mentioned equation does give the magnitude.

import numpy as np
norm = np.linalg.norm

p1 = np.array([0,0,0])
p2 = np.array([10,0,3])

p3 = np.array([6, -3, 5])

ppDist = np.abs(norm(np.cross(p2-p1, p1-p3)))/norm(p2-p1)
print(ppDist) # 4.28888043816146

ppDist vector in reality would be (0.8807339448,3,-2.935779817) such that its norm() is 4.288880438

Here's a quick visualization - enter image description here

jared
  • 4,165
  • 1
  • 8
  • 31
beta green
  • 113
  • 10

3 Answers3

2

Find projection of P13 vector on P12 and the vector difference of P13 from the projection would be the perpendicular vector.

P12 = p2-p1
P13 = p3-p1

proj_13over12 = np.dot(P13, P12)*P12/norm(P12)**2 # array([6.88073394, 0.        , 2.06422018])
perpendicular = proj_13over12 - P13

print(perpendicular) # [ 0.88073394  3.         -2.93577982]

enter image description here

beta green
  • 113
  • 10
  • 2
    Interesting approach. Note that `np.linalg.norm(P12)**2` is equivalent to `np.dot(P12, P12)`, but the later is faster to compute than the former (due to the implicit square root to compute and internal Numpy overheads). – Jérôme Richard Jul 09 '23 at 00:07
0
python
import numpy as np

def distance_point_to_line(p1, p2, p3):
    v = p2 - p1
    w = p3 - p1
    proj_w_v = np.dot(w, v) / np.dot(v, v) * v
    u = w - proj_w_v
    d = np.linalg.norm(u)
    p = p3 - u
    return d, p
gene
  • 1
  • 3
  • your perpendicular vector returns `array([6.88073394, 0 , 2.06422018]))`, by apparent observation of the input, you could figure out that `y` value of the vector must be `3` not `0` – beta green Jul 08 '23 at 18:17
0

I'm going to give a slightly altered solution to OP's solution. This is because in OP's solution, the meaning of what they call proj_13over12 is missing. What we want from the projection is actually the point that lies on the line defined by p1 and p2 which creates a vector perpendicular to that vector when connecting that point with p3. In the special case here where p1 = (0,0,0), we get that result, but the value does not change when the points are moved. For example, with the given points, the perpendicular point is at [6.88073394, 0., 2.06422018]. But if we translate all the points by 1 in each direction, we should get a perpendicular point that was translated by the same amount, but we do not.

import numpy as np
from numpy.linalg import norm

p1 = np.array([0, 0, 0]) + 1
p2 = np.array([10, 0, 3]) + 1
p3 = np.array([6, -3, 5]) + 1

P12 = p2-p1
P13 = p3-p1

proj_13over12 = np.dot(P13, P12)*P12/norm(P12)**2 
print(proj_13over12) # still [6.88073394, 0., 2.06422018]

I'm going to use some different notation that makes more sense to me. First, the vector that connects p1 to p2 will be called v1 (for vector 1), and the vector connecting p1 to p2 will be called v2. The result of the projection will be called projv2_v1 (this somewhat mirrors the common notation where the vector being projected onto is a subscript of the word "proj"). The perpendicular point will be called pperp.

We compute the projection using the equation v1 (v1 dot v2)/norm(v1)^2, which can also be written as v1 (v1 dot v2)/(v1 dot v1). (Mathematical interlude: The projection works because v1/norm(v1) gives a unit vector in the v1 direction while (v1 dot v2)/norm(v1) gives the magnitude of the v2 vector that is in the v1 direction.) Now, after computing the projection, pperp is found by adding p1 back (we created the vectors by taking the difference with p1, essentially translating everything to the origin, so this step translates back to the actual location in space). And then the distance is found by finding the norm of the vector connecting p1 to pperp. With these changes, the results all have simple physical interpretations.

import numpy as np
from numpy.linalg import norm

p1 = np.array([0,0,0])
p2 = np.array([10,0,3])
p3 = np.array([6, -3, 5])

v1 = p2 - p1
v2 = p3 - p1
projv2_v1 = np.dot(v1,v2)*v1/np.dot(v1, v1)
pperp = projv2_v1 + p1
print(pperp) # [6.88073394, 0., 2.06422018]
dist = norm(pperp - p3)
print(dist)  # 4.28888043816146

Now, if we were to translate all the points as before, we get pperp = [7.88073394, 1., 3.06422018], which makes sense because it would also be translated the same amount.

jared
  • 4,165
  • 1
  • 8
  • 31
  • Thanks for your feedback. Please note that projv2_v1 is a vector (start=p1, end=pperp vector end point) whereas pperp vector itself got start=origin i.e. 0,0,0. [link](https://i.imgur.com/f2RBnmk.png) – beta green Jul 09 '23 at 07:00
  • 1
    @betagreen I agree that `projv2_v1` is a vector, but I'd consider `pperp` a point, which is why I named it like I did. – jared Jul 09 '23 at 13:55