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
x1 + a dx + b dy - x3 = 0

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

a = (dy y3 - dy y1 + dx x3 - dx x1) / (dy^2 + dx^2)
b = (dy x3 - dy x1 - dx y3 + dx y1) / (dy^2 + dx^2)

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.

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.

---

A Matplotlib Example

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 are after  right angles, anything else would be wrong
ax.set_aspect(1)

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

Answer 2

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

Answer 3

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