Python: point on a line closest to third point

Question

I have a line/vector between two XY points (p1 and p2) and a third XY point (p3) that is outside the line. According to this post I know how to get the distance of that point to the line. But what I'm actually looking for is a point (p4) on that line that is in a minimum distance (d) to the third point (p3). I found this post, but I feel it's not the correct solution. Maybe there's something included in Numpy or Python?

According to @allo I tried the following. You can download my code as Python file or Jupyter Notebook (both Python3).

points = [[1, 1], [3, 1], [2.5, 2], [2.5, 1]]
import matplotlib.pyplot as plt
%matplotlib inline

fig, ax = plt.subplots()
fig.set_size_inches(6,6)

x, y = zip(*points[:2])
l1, = ax.plot(x,y, color='blue')
scatter1 = ax.scatter(x=x,y=y, color='blue', marker='x', s=80, alpha=1.0)

x, y = zip(*points[2:])
l2, = ax.plot(x,y, color='red')
scatter2 = ax.scatter(x=x,y=y, color='red', marker='x', s=80, alpha=1.0)

p1 = Vector2D(*points[0])
p2 = Vector2D(*points[1])
p3 = Vector2D(*points[2])

p1p2 = p2.sub_vector(p1)
p1p3 = p3.sub_vector(p1)

angle_p1p2_p1p3 = p1p2.get_angle_radians(p1p3)
length_p1p3 = p1p3.get_length()
length_p1p2 = p1p2.get_length()

p4 = p1.add_vector(p1p2.multiply(p1p3.get_length()/p1p2.get_length()).multiply(math.cos(p1p2.get_angle_radians(p1p3))))

#p4 = p1 + p1p2 * length(p1p3)/length(p1p2)*cos(angle(p1p2, p1p3))

p4 = p1.add_vector(p1p2.multiply(length_p1p3/length_p1p2*math.cos(angle_p1p2_p1p3)))
p4

Which results in p4 = (1.8062257748298551, 1.0) but should obviously be (2.5, 1.0).

Answer 1

Analytical Geometry

Let's start with the assigned line, we define the line in terms of two points on it (x1, y1) and (x2, y2).

With dx = x2-x1 and dy = y2-y1 we can formally write every point on the line as (x12, y12) = (x1, y1) + a*(dx, dy) where a is a real number.

Using an analogous notation a point on the line passing in (x3, y3) and perpendicular to the assigned one is (x34, y34) = (x3, y3) + b*(-dy, +dx).

To find the intersection we have to impose (x12, y12) = (x34, y34) or (x1, y1) + a*(dx, dy) = (x3, y3) + b*(-dy, +dx).

Writing separately the equations for x and y

y1 + a dy - y3 - b dx = 0   | +dy  -dx |   | a |   | y3-y1 |
                          → |          | × |   | = |       |
x1 + a dx + b dy - x3 = 0   | +dx  +dy |   | b |   | x3-x1 |

it is a linear system in a and b whose solution is

| a |       1       | +dy  +dx |   | y3-y1 |    1 | dy·δy+dx·δx |
|   | = --------- × |          | × |       | =  - |             |
| b |   dx² + dy²   | -dx  +dy |   | x3-x1 |    Δ | dy·δx-dx·δy |

dx=x2-x1, δx=x3-x1, dy=y2-t1, δy=y3-y1, Δ=dx²+dy².

The coordinates of the closest point to (x3, y3) lying on the line are (x1+a*dx, y1+a*dy) — you need to compute only the coefficient a = (dy·δy+dx·δx)/Δ.

Numerically speaking, the determinant of the linear system is Δ = dx**2+dy**2 so you have problems only when the two initial points are extremely close to each other with respect to their distance w/r to the third point.

Python Implementation

We use a 2-uple of floats to represent a 2D point and we define a function whose arguments are 3 2-uples representing the points that define the line (p1, p2) and the point (p3) that determines the position of p4 on said line.

In [16]: def p4(p1, p2, p3):
    ...:     x1, y1 = p1
    ...:     x2, y2 = p2
    ...:     x3, y3 = p3
    ...:     dx, dy = x2-x1, y2-y1
    ...:     det = dx*dx + dy*dy
    ...:     a = (dy*(y3-y1)+dx*(x3-x1))/det
    ...:     return x1+a*dx, y1+a*dy

To test the implementation I am using the three points used by the OP to demonstrate their issues with this problem:

In [17]: p4((1.0, 1.0), (3.0, 1.0), (2.5, 2))
Out[17]: (2.5, 1.0)

It seems that the result of p4(...) coincides with the OP expectation. Also, when p1, p2, p3 are collinear, the function return p3 as expected from geometrical considerations.

---

A Matplotlib Example

Note that your plot will show a right angle where a right angle is expected if and only if you set the aspect ratio of the Axes to 1, i.e., the scaling of the axes is the same.

import matplotlib.pyplot as plt

def p(p1, p2, p3):
    (x1, y1), (x2, y2), (x3, y3) = p1, p2, p3
    dx, dy = x2-x1, y2-y1
    det = dx*dx + dy*dy
    a = (dy*(y3-y1)+dx*(x3-x1))/det
    return x1+a*dx, y1+a*dy

p1, p2, p3 = (2, 4), (7, 3), (1, 1)
p4 = p(p1, p2, p3)

fig, ax = plt.subplots()

# if we want to see right angles where they should be,
# the aspect ratio y/x must be equal to one,
ax.set_aspect(1)

plt.plot(*zip(p1, p2, p4, p3), marker='*')

Answer 2

What you want to do is a vector projection.

The edge p1p3 is rotated onto the edge p1p2 and you need to find the correct length of the segment p1p4. Then you can just use p1+FACTOR*p1p2 / length(p1p2). The needed factor is given by the cosine of the angle between p1p2 and p1p3. Then you get

p4 = p1 + p1p2 * length(p1p3)/length(p1p2)*cos(angle(p1p2, p1p3))

Here the two edge cases as example:

• The cosinus is 0 if p1p3 is orthogonal to p1p2, so p4 lies on p1.
• The cosinus is 1 when p1p3 lies on p1p2, so p1p2 is just scaled by length(p1p3)/length(p1p2) to get p1p4.

You further can replace the cosinus by the dot product dot(p1p2 / length(p1p2), p1p3 / length(p1p3).

You can find more details and nice illustrations in the wikibook about linear algebra.

Here an full example derived from your python code. I used numpy instead of Vector2D here.

points = [[1, 1], [3, 1], [2.5, 2], [2.5, 1]]
import matplotlib.pyplot as plt
%matplotlib inline

import numpy as np

fig, ax = plt.subplots()
fig.set_size_inches(6,6)

x, y = zip(*points[:2])
l1, = ax.plot(x,y, color='blue')
scatter1 = ax.scatter(x=x,y=y, color='blue', marker='x', s=80, alpha=1.0)

x, y = zip(*points[2:])
l2, = ax.plot(x,y, color='red')
scatter2 = ax.scatter(x=x,y=y, color='red', marker='x', s=80, alpha=1.0)

p1 = np.array(points[0])
p2 = np.array(points[1])
p3 = np.array(points[2])

p1p2 = p2 - p1
p1p3 = p3 - p1

p4 = p1 + p1p2 / np.linalg.norm(p1p2) * np.linalg.norm(p1p3) * ((p1p2/np.linalg.norm(p1p2)).T * (p1p3/np.linalg.norm(p1p3)))

p1, p2, p3, p4, p1p2, p1p3

We can shorten the p4 line a bit like this, using the linearity of the scalar product:

p4 = p1 + p1p2 / np.linalg.norm(p1p2) * ((p1p2/np.linalg.norm(p1p2)).T * (p1p3))
p4 = p1 + p1p2 / np.linalg.norm(p1p2)**2 * (p1p2.T * (p1p3))

Answer 3

Shapely's distance() function returns the minimum distance:

>>> from shapely.geometry import LineString as shLs
>>> from shapely.geometry import Point as shPt
>>> l = shLs([ (1,1), (3,1)])
>>> p = shPt(2,2)
>>> dist = p.distance(l)
1.0
>>> l.interpolate(dist).wkt
'POINT (2 1)'