8

If i have three points P1, P2, P3 with their coordinates(x,y)

P1(x,y) and P3(x,y) are coordinate of line(start, end) and P3 is a point need to be projected.

how can i find the coordinate of point r(x,y) which is projection of P3 over P1 and P2 enter image description here

Mohamad A Sallal
  • 614
  • 6
  • 12
  • Does this answer your question? [Projection of a point on line defined by 2 points](https://stackoverflow.com/questions/15232356/projection-of-a-point-on-line-defined-by-2-points) – Lesiak Apr 21 '20 at 10:55
  • No that is different and it use Microsoft.Maps.Location API, I am looking to solve it in Python and Numpy – Mohamad A Sallal Apr 21 '20 at 10:57
  • It is trivial to port it to python - use your representation of point – Lesiak Apr 21 '20 at 10:58
  • I am checking here because i need it to be vectorized. i have thousands of points need to be projected. and i did not mentioned my points are 3d (x,y,z) for simplicity – Mohamad A Sallal Apr 21 '20 at 11:08
  • In your question you say P1(x,y) and P3(x,y) , which clearly defines @D points, not you say you need 3D solution. Make up your mind :) – Lesiak Apr 21 '20 at 11:25
  • Are you trying to project your point to a line or a segment (has clearly defined start/end)? – ec2604 Apr 21 '20 at 12:36

2 Answers2

11

This solution extends to points with any geometric dimensions (2D, 3D, 4D, ...). It assumes all points are one dimensional numpy arrays (or two dimensional with one dimension shape 1). I am not sure if you require the projection to fall onto line segment or the extension of segment so I include both. You can pick whichever fits your question the best:

#distance between p1 and p2
l2 = np.sum((p1-p2)**2)
if l2 == 0:
  print('p1 and p2 are the same points')

#The line extending the segment is parameterized as p1 + t (p2 - p1).
#The projection falls where t = [(p3-p1) . (p2-p1)] / |p2-p1|^2

#if you need the point to project on line extention connecting p1 and p2
t = np.sum((p3 - p1) * (p2 - p1)) / l2

#if you need to ignore if p3 does not project onto line segment
if t > 1 or t < 0:
  print('p3 does not project onto p1-p2 line segment')

#if you need the point to project on line segment between p1 and p2 or closest point of the line segment
t = max(0, min(1, np.sum((p3 - p1) * (p2 - p1)) / l2))

projection = p1 + t * (p2 - p1)
Ehsan
  • 12,072
  • 2
  • 20
  • 33
  • yes this is what i am looking for, Thank you so much – Mohamad A Sallal Apr 21 '20 at 12:54
  • @MSallal You are welcome. If it solves your problem please accept the answer to help others find it useful. Thank you. – Ehsan Apr 21 '20 at 12:56
  • I am not sure if I understand the question. Using first equation for t in my post results in projection on extended line/segment line which is essentially perpendicular to it. – Ehsan Apr 22 '20 at 04:13
  • thank you for the answer, would it be possible to link the mathematical foundation of this answer? I found several parametric line equations, but have not decided which one is more convenient for implementation – JVGD May 26 '22 at 06:55
5

Numpy code that works for both 2D and 3D (based on https://gamedev.stackexchange.com/questions/72528/how-can-i-project-a-3d-point-onto-a-3d-line):

import numpy as np

def point_on_line(a, b, p):
    ap = p - a
    ab = b - a
    result = a + np.dot(ap, ab) / np.dot(ab, ab) * ab
    return result

A = np.array([2, 0])
B = np.array([4, 4])
P = np.array([1, 3])
projected = point_on_line(A, B, P)
print(projected)

Update

Plot:

A = np.array([ 10.5, 15.6 ])
B = np.array([ 2, 6 ])
P = np.array([ 18.561, -19.451])

projected = point_on_line(A, B, P) 
print(projected)
# [-3.35411076 -0.04699568]

plt.xlim(-20, 20)
plt.ylim(-20, 20)
plt.axis('equal')

x_values = [A[0], B[0]]
y_values = [A[1], B[1]]

plt.plot(B[0], B[1], 'ro')
plt.plot(A[0], A[1], 'ro')
plt.plot(P[0], P[1], 'ro')
plt.plot(x_values, y_values, 'b-')
plt.plot(projected[0], projected[1], 'rx')

Update 2

If you need the point to belong to the segment, you need to make a small amendment

def point_on_line(a, b, p):
    ap = p - a
    ab = b - a
    t = np.dot(ap, ab) / np.dot(ab, ab)
    # if you need the the closest point belonging to the segment
    t = max(0, min(1, t))
    result = a + t * ab
    return result
Lesiak
  • 22,088
  • 2
  • 41
  • 65
  • i tried this before, it is not showing the projected point check this notebook https://projected-msallal.notebooks.azure.com/j/notebooks/ProjectedPoint.ipynb – Mohamad A Sallal Apr 21 '20 at 11:34
  • I took a look at your notebook, you have 2 errors: 1: you call projected = ClosestPointOnLine(A, P, B) while it needs to be projected = ClosestPointOnLine(A, B, P). 2>You plot in an incorrect way – Lesiak Apr 21 '20 at 13:18
  • how can we make it perpendicular on line – Mohamad A Sallal Apr 22 '20 at 02:39
  • See update how to make it belong on the line. Note that accepted answer is basically the same algorithm using more basic functions: `np.sum((p1-p2)**2)` is equal to `np.dot(ab, ab)` etc. – Lesiak Apr 22 '20 at 07:11
  • 1
    lol. this is exactly the same answer :D. It is also good to check for `ab != 0` to avoid throwing division by zero error. – Ehsan Apr 22 '20 at 17:13