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.