Get Closest Point to a Line

get closest point to a line

Here's Ruby disguised as Pseudo-Code, assuming Point objects each have a x and y field.

def GetClosestPoint(A, B, P)

  a_to_p = [P.x - A.x, P.y - A.y]     # Storing vector A->P
  a_to_b = [B.x - A.x, B.y - A.y]     # Storing vector A->B

  atb2 = a_to_b[0]**2 + a_to_b[1]**2  # **2 means "squared"
                                      #   Basically finding the squared magnitude
                                      #   of a_to_b

  atp_dot_atb = a_to_p[0]*a_to_b[0] + a_to_p[1]*a_to_b[1]
                                      # The dot product of a_to_p and a_to_b

  t = atp_dot_atb / atb2              # The normalized "distance" from a to
                                      #   your closest point

  return Point.new( :x => A.x + a_to_b[0]*t,
                    :y => A.y + a_to_b[1]*t )
                                      # Add the distance to A, moving
                                      #   towards B

end

Alternatively:

From Line-Line Intersection, at Wikipedia. First, find Q, which is a second point that is to be had from taking a step from P in the "right direction". This gives us four points.

def getClosestPointFromLine(A, B, P)

  a_to_b = [B.x - A.x, B.y - A.y]   # Finding the vector from A to B
                                        This step can be combined with the next
  perpendicular = [ -a_to_b[1], a_to_b[0] ]
                                    # The vector perpendicular to a_to_b;
                                        This step can also be combined with the next

  Q = Point.new(:x => P.x + perpendicular[0], :y => P.y + perpendicular[1])
                                    # Finding Q, the point "in the right direction"
                                    # If you want a mess, you can also combine this
                                    # with the next step.

  return Point.new (:x => ((A.x*B.y - A.y*B.x)*(P.x - Q.x) - (A.x-B.x)*(P.x*Q.y - P.y*Q.x)) / ((A.x - B.x)*(P.y-Q.y) - (A.y - B.y)*(P.y-Q.y)),
                    :y => ((A.x*B.y - A.y*B.x)*(P.y - Q.y) - (A.y-B.y)*(P.x*Q.y - P.y*Q.x)) / ((A.x - B.x)*(P.y-Q.y) - (A.y - B.y)*(P.y-Q.y)) )

end

Caching, Skipping steps, etc. is possible, for performance reasons.

How to find closest point on line?

Infinite length:

If you have line with infinite length with start and direction, calculate the dot product of the line direction then multiply it by the direction and add the starting point to it.

public Vector2 FindNearestPointOnLine(Vector2 origin, Vector2 direction, Vector2 point)
{
    direction.Normalize();
    Vector2 lhs = point - origin;

    float dotP = Vector2.Dot(lhs, direction);
    return origin + direction * dotP;
}

Finite length:

If you have line with finite length with start to end positions, get the heading the perform a projection from the starting point to the. Also, use Mathf.Clamp to clap it just in case the line is off.

public Vector2 FindNearestPointOnLine(Vector2 origin, Vector2 end, Vector2 point)
{
    //Get heading
    Vector2 heading = (end - origin);
    float magnitudeMax = heading.magnitude;
    heading.Normalize();

    //Do projection from the point but clamp it
    Vector2 lhs = point - origin;
    float dotP = Vector2.Dot(lhs, heading);
    dotP = Mathf.Clamp(dotP, 0f, magnitudeMax);
    return origin + heading * dotP;
}

How to find the closest point on a line segment to an arbitrary point?

The general answer is to project the point onto the line.
One way to see it is to transform the point into the reference frame defined by your segment (p1 is the new origin (0, 0), p2 the new (1, 0)).
Then, you get rid of the new y coordinate (that's where the actual projection occurs) and transform the new point (x, 0) back into the original frame.

Concretely, you'll have to find the transformations.
The second one, from the new space into the original space is easy to write (just draw it on paper, you'll see):

x = (x2 - x1) * nx + (y2 - y1) * ny + x1
y = (y1 - y2) * nx + (x2 - x1) * ny + y1

But you can invert these equations to find (nx, ny) that correspond to a point (x, y).

When you do that, and assuming neither of us has made any mistakes nor typo, you should get something like:

dx = x2 - x1
dy = y2 - y1
d2 = dx*dx + dy*dy
nx = ((x3-x1)*dx + (y3-y1)*dy) / d2
return (dx*nx + x1, dy*nx + y1)

edit: If you actually have to find the closest point on the segment instead of the line, it is easy to find because if the projection falls within the segment, you have 0 <= nx <= 1 (it's a necessary and sufficient condition).
Before the return statement, you can just force nx to stay in this interval:

nx = min(1, max(0, nx))

reedit:
The statement above is equivalent to:

if nx<0:
    nx = 0
if nx>1:
    nx = 1

This way, the projection of the point on the line (which can be outside the segment) gets pushed back inside the segment (defined by 0 <= nx <= 1) at the closest point.

how to find the nearest LINESTRING to a POINT?

• your sample geometry is invalid for line strings, have modified
• it's simple to achieve with sjoin_nearest()

import geopandas as gpd
import shapely.wkt
import shapely.geometry

line_string = ["LINESTRING (-1.15.12 9.9, -1.15.13 9.93)", "LINESTRING (-2.15.12 8.9, -2.15.13 8.93)"]
# fix invalid wkt string...
line_string = ["LINESTRING (-1.15 9.9, -1.15 9.93)", "LINESTRING (-2.15 8.9, -2.15 8.93)"]
point = "POINT (5.41 3.9)"

gdf_p = gpd.GeoDataFrame(geometry=[shapely.wkt.loads(point)])
gdf_l = gpd.GeoDataFrame(geometry=pd.Series(line_string).apply(shapely.wkt.loads))

df_n = gpd.sjoin_nearest(gdf_p, gdf_l).merge(gdf_l, left_on="index_right", right_index=True)

df_n["distance"] = df_n.apply(lambda r: r["geometry_x"].distance(r["geometry_y"]), axis=1)

df_n

import numpy as np

def isectSphere(p0, p1, cloud):
        """
        >>> isectSphere([1, 0], [3, 0], [[0, -4], [2, 3]])
        1
        >>> isectSphere([1, 0, 0], [3, 0, 0], [[0, -4, 0], [2, 3, 0]])
        1
        """
        p0 = np.asarray(p0)
        p1 = np.asarray(p1)
        cloud = np.asarray(cloud)
        product = np.cross(cloud - p0, p1 - p0)
        if product.ndim == 2:
            distances = np.linalg.norm(product, axis=1)
        else:
            distances = np.abs(product)
        return distances.argmin()