Shortest distance between a point and a line segment

Question

I need a basic function to find the shortest distance between a point and a line segment. Feel free to write the solution in any language you want; I can translate it into what I'm using (Javascript).

EDIT: My line segment is defined by two endpoints. So my line segment AB is defined by the two points A (x1,y1) and B (x2,y2). I'm trying to find the distance between this line segment and a point C (x3,y3). My geometry skills are rusty, so the examples I've seen are confusing, I'm sorry to admit.

Solution 1

Eli, the code you've settled on is incorrect. A point near the line on which the segment lies but far off one end of the segment would be incorrectly judged near the segment. Update: The incorrect answer mentioned is no longer the accepted one.

Here's some correct code, in C++. It presumes a class 2D-vector class vec2 {float x,y;}, essentially, with operators to add, subract, scale, etc, and a distance and dot product function (i.e. x1 x2 + y1 y2).

float minimum_distance(vec2 v, vec2 w, vec2 p) {
  // Return minimum distance between line segment vw and point p
  const float l2 = length_squared(v, w);  // i.e. |w-v|^2 -  avoid a sqrt
  if (l2 == 0.0) return distance(p, v);   // v == w case
  // Consider the line extending the segment, parameterized as v + t (w - v).
  // We find projection of point p onto the line. 
  // It falls where t = [(p-v) . (w-v)] / |w-v|^2
  // We clamp t from [0,1] to handle points outside the segment vw.
  const float t = max(0, min(1, dot(p - v, w - v) / l2));
  const vec2 projection = v + t * (w - v);  // Projection falls on the segment
  return distance(p, projection);
}

EDIT: I needed a Javascript implementation, so here it is, with no dependencies (or comments, but it's a direct port of the above). Points are represented as objects with x and y attributes.

function sqr(x) { return x * x }
function dist2(v, w) { return sqr(v.x - w.x) + sqr(v.y - w.y) }
function distToSegmentSquared(p, v, w) {
  var l2 = dist2(v, w);
  if (l2 == 0) return dist2(p, v);
  var t = ((p.x - v.x) * (w.x - v.x) + (p.y - v.y) * (w.y - v.y)) / l2;
  t = Math.max(0, Math.min(1, t));
  return dist2(p, { x: v.x + t * (w.x - v.x),
                    y: v.y + t * (w.y - v.y) });
}
function distToSegment(p, v, w) { return Math.sqrt(distToSegmentSquared(p, v, w)); }

EDIT 2: I needed a Java version, but more important, I needed it in 3d instead of 2d.

float dist_to_segment_squared(float px, float py, float pz, float lx1, float ly1, float lz1, float lx2, float ly2, float lz2) {
  float line_dist = dist_sq(lx1, ly1, lz1, lx2, ly2, lz2);
  if (line_dist == 0) return dist_sq(px, py, pz, lx1, ly1, lz1);
  float t = ((px - lx1) * (lx2 - lx1) + (py - ly1) * (ly2 - ly1) + (pz - lz1) * (lz2 - lz1)) / line_dist;
  t = constrain(t, 0, 1);
  return dist_sq(px, py, pz, lx1 + t * (lx2 - lx1), ly1 + t * (ly2 - ly1), lz1 + t * (lz2 - lz1));
}

Here, in the function parameters, <px,py,pz> is the point in question and the line segment has the endpoints <lx1,ly1,lz1> and <lx2,ly2,lz2>. The function dist_sq (which is assumed to exist) finds the square of the distance between two points.

Solution 2

Here is the simplest complete code in Javascript.

x, y is your target point and x1, y1 to x2, y2 is your line segment.

UPDATED: fix for 0 length line problem from comments.

function pDistance(x, y, x1, y1, x2, y2) {

  var A = x - x1;
  var B = y - y1;
  var C = x2 - x1;
  var D = y2 - y1;

  var dot = A * C + B * D;
  var len_sq = C * C + D * D;
  var param = -1;
  if (len_sq != 0) //in case of 0 length line
      param = dot / len_sq;

  var xx, yy;

  if (param < 0) {
    xx = x1;
    yy = y1;
  }
  else if (param > 1) {
    xx = x2;
    yy = y2;
  }
  else {
    xx = x1 + param * C;
    yy = y1 + param * D;
  }

  var dx = x - xx;
  var dy = y - yy;
  return Math.sqrt(dx * dx + dy * dy);
}

Solution 3

This is an implementation made for FINITE LINE SEGMENTS, not infinite lines like most other functions here seem to be (that's why I made this).

Implementation of theory by Paul Bourke.

Python:

def dist(x1, y1, x2, y2, x3, y3): # x3,y3 is the point
    px = x2-x1
    py = y2-y1

    norm = px*px + py*py

    u =  ((x3 - x1) * px + (y3 - y1) * py) / float(norm)

    if u > 1:
        u = 1
    elif u < 0:
        u = 0

    x = x1 + u * px
    y = y1 + u * py

    dx = x - x3
    dy = y - y3

    # Note: If the actual distance does not matter,
    # if you only want to compare what this function
    # returns to other results of this function, you
    # can just return the squared distance instead
    # (i.e. remove the sqrt) to gain a little performance

    dist = (dx*dx + dy*dy)**.5

    return dist

AS3:

public static function segmentDistToPoint(segA:Point, segB:Point, p:Point):Number
{
    var p2:Point = new Point(segB.x - segA.x, segB.y - segA.y);
    var something:Number = p2.x*p2.x + p2.y*p2.y;
    var u:Number = ((p.x - segA.x) * p2.x + (p.y - segA.y) * p2.y) / something;

    if (u > 1)
        u = 1;
    else if (u < 0)
        u = 0;

    var x:Number = segA.x + u * p2.x;
    var y:Number = segA.y + u * p2.y;

    var dx:Number = x - p.x;
    var dy:Number = y - p.y;

    var dist:Number = Math.sqrt(dx*dx + dy*dy);

    return dist;
}

Java

private double shortestDistance(float x1,float y1,float x2,float y2,float x3,float y3)
    {
        float px=x2-x1;
        float py=y2-y1;
        float temp=(px*px)+(py*py);
        float u=((x3 - x1) * px + (y3 - y1) * py) / (temp);
        if(u>1){
            u=1;
        }
        else if(u<0){
            u=0;
        }
        float x = x1 + u * px;
        float y = y1 + u * py;

        float dx = x - x3;
        float dy = y - y3;
        double dist = Math.sqrt(dx*dx + dy*dy);
        return dist;

    }

Solution 4

In my own question thread how to calculate shortest 2D distance between a point and a line segment in all cases in C, C# / .NET 2.0 or Java? I was asked to put a C# answer here when I find one: so here it is, modified from http://www.topcoder.com/tc?d1=tutorials&d2=geometry1&module=Static :

//Compute the dot product AB . BC
private double DotProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] BC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    BC[0] = pointC[0] - pointB[0];
    BC[1] = pointC[1] - pointB[1];
    double dot = AB[0] * BC[0] + AB[1] * BC[1];

    return dot;
}

//Compute the cross product AB x AC
private double CrossProduct(double[] pointA, double[] pointB, double[] pointC)
{
    double[] AB = new double[2];
    double[] AC = new double[2];
    AB[0] = pointB[0] - pointA[0];
    AB[1] = pointB[1] - pointA[1];
    AC[0] = pointC[0] - pointA[0];
    AC[1] = pointC[1] - pointA[1];
    double cross = AB[0] * AC[1] - AB[1] * AC[0];

    return cross;
}

//Compute the distance from A to B
double Distance(double[] pointA, double[] pointB)
{
    double d1 = pointA[0] - pointB[0];
    double d2 = pointA[1] - pointB[1];

    return Math.Sqrt(d1 * d1 + d2 * d2);
}

//Compute the distance from AB to C
//if isSegment is true, AB is a segment, not a line.
double LineToPointDistance2D(double[] pointA, double[] pointB, double[] pointC, 
    bool isSegment)
{
    double dist = CrossProduct(pointA, pointB, pointC) / Distance(pointA, pointB);
    if (isSegment)
    {
        double dot1 = DotProduct(pointA, pointB, pointC);
        if (dot1 > 0) 
            return Distance(pointB, pointC);

        double dot2 = DotProduct(pointB, pointA, pointC);
        if (dot2 > 0) 
            return Distance(pointA, pointC);
    }
    return Math.Abs(dist);
} 

I'm @SO not to answer but ask questions so I hope I don't get million down votes for some reasons but constructing critic. I just wanted (and was encouraged) to share somebody else's ideas since the solutions in this thread are either with some exotic language (Fortran, Mathematica) or tagged as faulty by somebody. The only useful one (by Grumdrig) for me is written with C++ and nobody tagged it faulty. But it's missing the methods (dot etc.) that are called.

Solution 5

In F#, the distance from the point c to the line segment between a and b is given by:

let pointToLineSegmentDistance (a: Vector, b: Vector) (c: Vector) =
  let d = b - a
  let s = d.Length
  let lambda = (c - a) * d / s
  let p = (lambda |> max 0.0 |> min s) * d / s
  (a + p - c).Length

The vector d points from a to b along the line segment. The dot product of d/s with c-a gives the parameter of the point of closest approach between the infinite line and the point c. The min and max function are used to clamp this parameter to the range 0..s so that the point lies between a and b. Finally, the length of a+p-c is the distance from c to the closest point on the line segment.

Example use:

pointToLineSegmentDistance (Vector(0.0, 0.0), Vector(1.0, 0.0)) (Vector(-1.0, 1.0))

View more solutions

Comments

[Eli Courtwright]

Thanks - this is exactly the kind of code I was looking for. I've posted my own answer below, since I managed to put something together that works in current-era-browser-Javascript, but I've marked your answer as accepted because it's simple, well-written, easy-to-understand, and much appreciated.

[Eli Courtwright]

Hey, thanks! This would have been really helpful when I was coding my solution, and it'll be a nice reference for the next time I need to do something like this. I'm marking this accepted because your Python solution is short, explains its steps, and is a complete working example that I can copy/paste for testing.

[nit]

isn't this computing the distance of a point to a line instead of the segment?

[Benji Mizrahi]

How about the projection length? What would its algorithm be?

[Benji Mizrahi]

I think I found my own answer: return fabs(ac[0]*t[0]+ac[1]*t[1])

[Grumdrig]

This code has a bug. A point near the line on which the segment lies, but far off one end of the segment, would be incorrectly judged to be near the segment.

[Grumdrig]

This is indeed the distance to the line the segment is on, not to the segment.

[Eli Courtwright]

Interesting, I'll look into this the next time I'm working on this codebase to confirm your assertion. Thanks for the tip.

[quano]

Isn't this missing the dot-method? In any case, it is easy to calculate: vec1.x * vec2.x + vec1.y * vec2.y

[quano]

This doesn't seem to work. If you've got a segment of (0,0) and (5,0), and try against point (7,0), it will return 0, which isn't true. The distance should be 2.

[phkahler]

He's failed to consider the case where the projection of the point onto the segment is outside the interval from A to B. That might be what the questioner wanted, but not what he asked.

[Sambatyon]

This is not what was originally asked.

[Frederik]

Sorry, but I tried this and it still gives me the results as if the line was extending into infinity. I've found Grumdig's answer to work, though.

[quano]

In that case you're using it wrong or meaning something else with non-infinite. See an example of this code here: boomie.se/upload/Drawdebug.swf

[LinuxRsa]

As others have said, this algorithm is for a line, not a line segment. The quano comment explains the difference.

[Trinidad]

Does not answer the question. This only works for lines (the ones that extend infinitely in space) not line segments (which have a finite length).

[Kromster]

Looks like a mistake in code or something, I get the same result as Frederik/

[quano]

Strange. Well, in the swf example I use exactly this code (it's copy pasted), and if the behavior in that swf is what you're looking for, this is the code.

[Rudolf Meijering]

Thanks, this Matlab code indeed calculates the shortest distance to the line SEGMENT and not the distance to the infinite line on which the segment lies.

[M Katz]

I've added a fleshed-out version of this as a separate answer.

[miguelSantirso]

The choice of variable names is far from good (p2, something, u, ...)

[awolf]

Thanks @Grumdrig, your javascript solution was spot on and a huge time saver. I ported your solution to Objective-C and added it below.

[Vladimir Obrizan]

I've tried the Python version of the function and found that it shows incorrect results if the parameters are integers. distAnother(0, 0, 4, 0, 2, 2) gives 2.8284271247461903 (incorrect). distAnother(0., 0., 4., 0., 2., 2.) gives 2.0 (correct). Please be aware of this. I think the code can be improved to have float conversion somewhere.

[RenniePet]

Thanks for posting this. But it looks like there's an obvious optimization possible in the last method: Don't compute dist until after it's determined that it's needed.

[RenniePet]

Thanks for posting this. Very well structured and commented and formatted - almost made me forget how much I dislike C++. I've used this to make a corresponding C# version, which I've now posted here.

[RenniePet]

"above function" is an ambiguous reference. (Irritates me because sometimes this answer is shown beneath my answer.)

[gongzhitaao]

@JaredMcAteer shouldn't the float equality be determined by another function instead of directly ==?

[Grumdrig]

We're really just trying to avoid a divide by zero there.

[user1815201]

Of all the code I've seen to solve this problem, I like this one the best. It is very clear and easy to read. The math behind it though, is a little bit mystical. What does the dot-product divided by the length squared really represent, for example?

[Totti]

You haven't defined mx1, b or m and your solution is not for line segments.

[wmac]

Excuse me for asking. If the distance is 0 does it mean the point has cut the line segment?

[Grumdrig]

Yes, the point should be somewhere on the line segment (possibly at an endpoint) in that case.

[Gregir]

I get 'nan' returned from this line. Any idea why? (Thanks for typing this up in Obj-C, by the way!) return dist2(p, CGPointMake(v.x + t * (w.x - v.x), v.y + t * (w.y - v.y)))

[Gregir]

This appears to work well for me. Thanks for converting.

[Logan Pickup]

The dot product divided by length squared gives you the projection distance from (x1, y1). This is the fraction of the line that the point (x,y) is closest to. Notice the final else clause where (xx, yy) is calculated - this the projection of the point (x,y) onto the segment (x1,y1)-(x2,y2).

[Senseful]

sqrtf() is squaring x, not getting its square root

[awolf]

@Senseful Not sure what you mean. sqrtf is square root. developer.apple.com/library/mac/documentation/Darwin/Referen‌​ce/…

[Senseful]

@awolf: Take a look at the first line of code above. It defines the method sqrtf(x) = x*x.

[awolf]

@Senseful thanks, it was misnamed rather than performing the wrong operation.

[Blair Holloway]

I think the last line is incorrect, and should read: (a + p - c).Length

[HostedMetrics.com]

The check for line segments of length 0 is too far down in the code. 'len_sq' will be zero and the code will divide by 0 before it gets to the safety check.

[user1132959]

This is will crash if the line points are the same. Read comment from @Heliodor.

[Anne van Rossum]

Kudos for actually writing a formula down. However, this is the formula to calculate the distance to a line, not to a line segment. Picking (x0,y0)=(-10,0), (x1,y1)=(0,0), and (x2,y2)=(10,0) gives a distance of 0, while it should have been 10.

[Joshua]

Updated to fix the 0 length line problem.

[Nicolas Miari]

@Grumdrig (and to the people who ported into other languages): This is a really handy, short piece of code. Have you considered making this answer into a Gist? (gist.github.com)

[Steve Johnson]

Works like a charm!! Saved me countless hours. Thanks so much!!

[Nat]

It's worth noticing, that (xx, yy) is location of closest point. I've edited a bit your code, so it return both the point and distance, refactored names so they describe what is what and provided example at: stackoverflow.com/a/28028023/849616.

[mikkoma]

That still does not fully fix the issue. One way to correct the function would be to redefine lambda and p as let lambda = (c - a) * d / (s * s) and let p = a + (lambda |> max 0.0 |> min 1.0) * d, respectively. After that the function returns correct distance e.g. for the case where a = (0,1), b = (1,0) and c = (1,1).

[xaxxon]

Can someone give a more descriptive definition of t and projection please? specifically how they relate to whether the projection lands on the line segment or not?

[Grumdrig]

The projection of point p onto a line is the point on the line closest to p. (And a perpendicular to the line at the projection will pass through p.) The number t is how far along the line segment from v to w that the projection falls. So if t is 0 the projection falls right on v; if it's 1, it's on w; if it's 0.5, for example, then it's halfway between. If t is less than 0 or greater than 1 it falls on the line past one end or the other of the segment. In that case the distance to the segment will be the distance to the nearer end.

[ViniBiso]

The distance returned in this function is in what? Meters, Kilometers?

[Joshua]

Meters. It is returned in meters.

[mengg]

I think the second condition if (param > 1) should be if(param>len_sq). Any one?

[Metal450]

The comment on DotProduct says it's computing AB.AC, but it's computing AB.BC.

[ChrisJJ]

In the JS version, distance() is redundant, no?

[1j01]

Yeah, in both versions it's unused. Feel free to remove it.

[WDUK]

Tried this code, doesn't seem to work quite correctly. Seems to get the wrong distance some times.

[clankill3r]

How can I do this in 3d?

[FirstOne]

JS code converted to PHP: http://stackoverflow.com/a/43571942/4577762

[urschrei]

@grumdrig It looks like this function returns a positive distance if the point lies on the line segment and the numbers are small; i.e. if p = (0.0001, 0.0), v = (0.0, 0.0), w = (0.004, 0.0), it returns 1.3552527156068805e-20 (see jsfiddle: jsfiddle.net/qg0zerkc). Is this a floating-point stability issue?

[Grumdrig]

@urschrei That would be floating point imprecision - if you're testing for whether a point is on the line you'd want to check that the distance is less than some small value adequate to your needs, rather than equality with zero.

[yolo sora]

nice idea, but height will be distance to line, in some cases distance to segment = distance to segment's start or end.

[Admin]

What does it mean when the param < 0 or > 1?

[Admin]

Can someone who understands each step add comments or provide an explanation? I know we're calculating the projection, I'm not sure why we're diving again by the norm and why we're checking of it's < 0 or > 1 and then finally what is dx, dy? Thanks

[Admin]

Let's say the line segment was 2 lat/lon points and the point was another lat/lon point, how could I get the distance returned in miles?

[Grumdrig]

@RogiSolorzano That's a 3D problem (earth isn't flat) so if you want accuracy this isn't the solution for you. If you want a very rough estimate for points not too far apart there's very very roughly 60 miles per degree.

[Grumdrig]

Oops - didn't notice someone had supplied a 3D version. @RogiSolorzano, you'll need to convert the lat,long coordinates into x,y,z coordinates in 3-space first.

[Admin]

Thank you @Grumdrig

[Sean]

@nevermind, let's call our point p0 and the points that define the line as p1 and p2. Then you get the vectors A = p0 - p1 and B = p2 - p1. Param is the scalar value that when multiplied to B gives you the point on the line closest to p0. If param <= 0, the closest point is p1. If param >= 1, the closest point is p1. If it's between 0 and 1, it's somewhere between p1 and p2 so we interpolate. XX and YY is then the closest point on the line segment, dx/dy is the vector from p0 to that point, and finally we return the length that vector.

[Alcamtar]

Thank you! Needed to find the intersection in Java, and this translated perfectly!

[Martin Asenov]

so you are saying the last snippet is in Java, huh ?

[Morris Franken]

The implementation can even be sped up by avoiding to compute the projection. Since the distance from p to line v,w is equal to sqrt(|p-v|^2 - (|w-v|*t)^2). Thus the last 2 lines of the c++ implementation can be replaced by return sqrt(length_squared(p, v) - l2 * (t * t));

[Humoyun Ahmad]

brilliant.org/wiki/…

[Humoyun Ahmad]

brilliant.org/wiki/…

[Nolo]

For anyone who needs an explanation for why this code works, please see this wikipedia permalink for vector projection and refer to the image at the beginning of the article for reference.

[Soonts]

It works but it’s inefficient. The code will do branching anyway, for min/max. When min or max clamps t, there’s no need to calculate the rest of the formula, you already know you need the distance between p and v or w.

[Eric Duminil]

@Joshua it doesn't necessarily return meters. It returns the same unit that is used for the input coordinates.

[Clonkex]

@ViniBiso In case it wasn't obvious from the other comments that jokingly replied "metres", the distance is returned in whatever you put in. So if you feed in metres, it will return metres. If you feed in miles, it will return miles. If you feed in shklobmorphs, it will return shklobmorphs. Same as virtually all geometric maths.

[SteakOverflow]

The cross product by definition returns a vector but returns a scalar here.

[Kamran Bashir]

You saved my day. Can you please tell me which equation you have used here? I had tried the shortest distance from line to a point but it didn't work for line segments.

[Aaron Franke]

What about the special case where one of the points is (0, 0)?

[user7924113]

anyone has 3D implementation in python?

[fede s.]

All these are calculating t with something equivalent to dot(p - v, w - v) / l2), which is t = u.v / || u ||². Shouldn't they be using the square root of l2? Meaning, the segment's length, IOW the magnitude of the vector w-v, and not its square. Am I missing something?

[Admin]

How can I make it so that it can also work with lines that extend forever and lines that have one end that extends forever?

[Savannah Madison]

I think the sqrt can be avoided since we apply this for all points and use the return value

[JasonG]

I found something odd, where two lines will have equal distance calculations for a point in fairly complex space: > pDistance(-79.73851776123047, 43.65549087524414, -79.7388727891, 43.6549957985, -79.7388654821, 43.6549894151); 0.0006092173427681732 > pDistance(-79.73851776123047, 43.65549087524414, -79.7396021293, 43.6545140287, -79.7388727891, 43.6549957985); 0.0006092173427681732 I have to try with another approach but I don't expect the same distance to appear as often as it does.

[Goku]

@Grumdrig Having trouble guessing what dist_sq(lx1, ly1, lz1, lx2, ly2, lz2) does. Does it find the distance between the two points (sum of sqares and then sqare root)? or does it not take a square root? I tried using the code and was getting strange results when I didnt take a square root in definitio of the dist_sq() function

[Grumdrig]

@Goku It finds the square of the distance between two points. (Sum of squares of component differences.) I've added some explanation to the end of the answer.