Point-to-Line and Point-to-Segment Distance, in depth: projection, cross products, clamping and robust code

Setting up with vectors

Call the segment endpoints A and B, the query point P, and define two vectors: d = B - A (the segment direction) and v = P - A (from A to the point). Every point on the infinite line through A and B can be written as A + t*d for some real number t. The segment is the part with t between 0 and 1; a ray from A towards B is the part with t of at least 0.

The closest point on the line is the foot of the perpendicular from P, the point where the vector from it to P is orthogonal to d. Writing that condition as a dot product equal to zero and solving gives the projection parameter:

t = dot(v, d) / dot(d, d)        # where P projects onto the line
closest_on_line = A + t * d

The numerator measures how far v runs along d; dividing by the squared length of d turns it into a fraction of the segment. A value of 0.5 means the foot lies at the midpoint, a negative value means it lies before A, and a value above 1 means it lies past B. That one number is the whole story of point-to-segment distance.

Distance to an infinite line

For the infinite line you do not need the foot point at all. In 2D, the cross product (the determinant d.x*v.y - d.y*v.x) equals the area of the parallelogram spanned by d and v. That area is also the base length times the height, and the height is exactly the distance from P to the line. So the distance is the absolute value of the cross product divided by the length of d.

The sign of the cross product is useful on its own: positive means P lies to the left of the direction from A to B, negative to the right, zero on the line. This is the orientation test that drives convex hull algorithms and segment intersection, and it needs no division or square root, so it can be computed exactly on integer coordinates.

If the line is given in implicit form, a*x + b*y + c = 0, the distance is |a*px + b*py + c| / sqrt(a*a + b*b). Normalise a and b once if you evaluate many points against the same line, and the formula becomes two multiply-adds and an absolute value per point.

In 3D the parallelogram argument still works with the vector cross product: the distance is the length of cross(d, v) divided by the length of d. For a=(0,0,0), b=(2,2,1) and p=(1,3,4), the cross product is (5,-7,4) with length the square root of 90, and the length of d is 3, so the distance to the line is the square root of 10, about 3.162. The projection parameter is 12/9, about 1.333, so the segment's closest point is B and the segment distance is the square root of 11, about 3.317. In higher dimensions use the projection form; it works in any dimension.

Distance to a segment: clamp t

For a segment, clamp t to the interval from 0 to 1 and measure to the resulting point. If t is below 0 the nearest point is A; above 1 it is B; otherwise it is the foot of the perpendicular. Clamping is valid because the squared distance along the line is a parabola in t with its minimum at the unclamped t, so on an interval the minimum is at the nearest end of the interval.

import math

def closest_point_on_segment(p, a, b):
    """Return (closest_point, t) with t clamped to [0, 1]."""
    dx, dy = b[0] - a[0], b[1] - a[1]
    vx, vy = p[0] - a[0], p[1] - a[1]
    dd = dx * dx + dy * dy
    if dd == 0.0:                      # degenerate segment: A == B
        return a, 0.0
    t = (vx * dx + vy * dy) / dd
    t = 0.0 if t < 0.0 else 1.0 if t > 1.0 else t
    return (a[0] + t * dx, a[1] + t * dy), t

def point_segment_distance(p, a, b):
    (cx, cy), _ = closest_point_on_segment(p, a, b)
    return math.hypot(p[0] - cx, p[1] - cy)

def point_line_distance(p, a, b):
    dx, dy = b[0] - a[0], b[1] - a[1]
    length = math.hypot(dx, dy)
    if length == 0.0:                  # no line is defined; fall back to point distance
        return math.hypot(p[0] - a[0], p[1] - a[1])
    return abs(dx * (p[1] - a[1]) - dy * (p[0] - a[0])) / length

When you compare distances, for example to find the nearest of many segments, compare squared distances and take one square root at the end, or none if you only need the winner. For large batches the same arithmetic vectorises cleanly:

import numpy as np

def segment_distances(p, A, B):
    """Distance from point p (2,) to each of n segments A[i]->B[i], arrays of shape (n, 2)."""
    d = B - A
    v = p - A
    dd = np.einsum("ij,ij->i", d, d)
    t = np.divide(np.einsum("ij,ij->i", v, d), dd, out=np.zeros_like(dd), where=dd > 0)
    t = np.clip(t, 0.0, 1.0)
    closest = A + t[:, None] * d
    return np.linalg.norm(p - closest, axis=1)

Worked example: one segment, four points

Take the segment from A=(1,1) to B=(7,4). Then d=(6,3) and the squared length is 45. Four query points:

Work through P1. v=(4,5), so the dot product is 6*4 + 3*5 = 39 and t = 39/45, about 0.867, inside the segment. The cross product is 6*5 - 3*4 = 18, and dividing by the square root of 45 (about 6.708) gives 2.683. As a check, the foot is A + 0.867*d = (6.2, 3.6), and the vector from it to P1 is (-1.2, 2.4), whose length is the square root of 7.2, also 2.683.

P2 is the instructive case. Its cross product is also 18, so its distance to the infinite line is the same 2.683, but t is 1.533, past B. The true distance to the segment is the distance to B: the vector (2, 4) has length the square root of 20, about 4.472. Code that uses the line formula for a segment reports 2.683, 40 percent short of the true distance. P3 projects before A, and P4 lies on the segment itself, with a cross product of exactly zero.

Numerical robustness

Degenerate segments. When A equals B the denominator is zero. The code above returns A as the closest point. Do not replace the exact zero test with a tiny epsilon unless you also handle the nearly-degenerate case, where t can be a huge number before clamping; the clamped result is still correct, but any unclamped use of t is not.

Square roots and overflow. math.hypot avoids intermediate overflow and underflow that sqrt(x*x + y*y) suffers for very large or very small components. Avoid the square root entirely when comparing against a threshold: test dist_sq <= r*r.

Cancellation far from the origin. Coordinates such as projected map positions in the millions of metres leave few significant bits for the differences you actually care about. Subtract a local origin (for example A) before doing anything else, which the vector form does naturally, and keep coordinates in double precision.

Exact decisions on integers. If inputs are integers and you need an exact yes or no, such as "is P within distance r of the segment" in a competitive programming problem or a grid, avoid floating point. Compare the dot product with 0 and with the squared length to classify the case; for the interior case the squared distance is cross*cross / dd, so test cross*cross <= r*r*dd in integers. Python integers do not overflow; in C or Java widen to 64 or 128 bits.

Units. Latitude and longitude are angles, not Cartesian coordinates. For short distances project to a local metric plane first; for long ones use great-circle geometry, where these formulas do not apply.

Where it is used

Polyline simplification. The Ramer-Douglas-Peucker algorithm keeps the point farthest from the segment joining a run's endpoints if that distance exceeds a tolerance, then recurses on both halves. Some implementations measure to the infinite line, which can keep or drop points differently for paths that double back, so check which your library uses.

Hit-testing and snapping. A click selects an edge if its distance to the edge segment is within a few pixels. Map matching snaps a GPS fix to the nearest road segment, using the clamped closest point as the snapped position and t to interpolate attributes such as distance along the road.

Nearest segment among many. Testing every segment is linear per query. For many queries, put segment bounding boxes into a grid or an R-tree, or index segment midpoints in a k-d tree and expand the search radius by half the longest segment, then run the exact test on the candidates.

Building block for harder queries. In 2D, segment-to-segment distance for non-intersecting segments is the minimum of the four point-to-segment distances from each endpoint to the other segment; in 3D that shortcut fails for skew segments, whose closest points can both be interior, so solve for two clamped parameters instead. Polygon distance, swept collision and capsule tests in physics engines reduce to the same primitive. The Graham scan walkthrough shows the orientation test from the previous section at work.

Failure modes

• Using the line formula when the problem means a segment, as P2 shows. Write separate functions with names that say which one they compute.
• Dividing by zero on repeated points, common in GPS traces and user-drawn paths. Deduplicate consecutive points or handle the zero-length case explicitly.
• Taking square roots inside a nearest-neighbour loop, which is both slower and adds rounding to comparisons that could have been exact.
• Mixing degrees with metres, or measuring in screen pixels after a non-uniform zoom.
• Epsilon comparisons with a fixed absolute tolerance that is too small for large coordinates and too large for small ones. Scale tolerances to the data, or use exact arithmetic for decisions.
• Assuming t is meaningful when the segment is degenerate; it is not, so do not interpolate attributes with it.

Trade-offs

The projection form is the most general: it works in any dimension, yields the closest point and t, and costs a few multiplications and one division. The cross-product form is cheaper when you only need the distance to a line, gives the side of the line for free, and has an exact integer version, but it is tied to 2D and 3D. Precomputing the normalised implicit form is fastest when one line is tested against many points.

Floating point is fine for measuring; exact arithmetic is worth it for decisions that must be consistent, such as which side of a line a point falls on in a mesh or a hull. Mixing the two, measuring in floats and then branching on tiny differences, is how geometric code becomes non-deterministic across platforms.

What to do next

1. Decide whether each call site needs a line, a ray or a segment, and name the functions accordingly.
2. Implement the projection form with clamping and an explicit degenerate-segment branch.
3. Compare squared distances in loops, and use hypot for the final value.
4. Test with points before A, past B, on the segment and on a zero-length segment, using the table above as fixtures.
5. Translate coordinates to a local origin, and project latitude and longitude before using these formulas.
6. Use the integer cross-product test for exact side-of-line and within-radius decisions.
7. For many segments, add a spatial index and run the exact test only on candidates.