Distance between point and line segment

[arah]

I'm having one of those brain fart days. :( (Lack of sleep on a Friday and all that.)

Could anyone please help by sending in a nice and simple algorithm for determining the shortest distance between a point (px,py) and a line segment(x1,y1)-(x2,y2)? Thanks. :)

[Philip Nicoletti]

Here is a function that calculates both the distance to the line segment and the distance to the line (assuming infinite extent).
I have a sample OpenGL graphics project I can attach which tests the function out if you want (as well as a couple of other small "mathematical" tasks such as this).


void DistanceFromLine(double cx, double cy, double ax, double ay ,
double bx, double by, double &distanceSegment,
double &distanceLine)
{

//
// find the distance from the point (cx,cy) to the line
// determined by the points (ax,ay) and (bx,by)
//
// distanceSegment = distance from the point to the line segment
// distanceLine = distance from the point to the line (assuming
// infinite extent in both directions
//

/*

Subject 1.02: How do I find the distance from a point to a line?


Let the point be C (Cx,Cy) and the line be AB (Ax,Ay) to (Bx,By).
Let P be the point of perpendicular projection of C on AB. The parameter
r, which indicates P's position along AB, is computed by the dot product
of AC and AB divided by the square of the length of AB:

(1) AC dot AB
r = ---------
||AB||^2

r has the following meaning:

r=0 P = A
r=1 P = B
r<0 P is on the backward extension of AB
r>1 P is on the forward extension of AB
0<r<1 P is interior to AB

The length of a line segment in d dimensions, AB is computed by:

L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 + ... + (Bd-Ad)^2)

so in 2D:

L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )

and the dot product of two vectors in d dimensions, U dot V is computed:

D = (Ux * Vx) + (Uy * Vy) + ... + (Ud * Vd)

so in 2D:

D = (Ux * Vx) + (Uy * Vy)

So (1) expands to:

(Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
r = -------------------------------
L^2

The point P can then be found:

Px = Ax + r(Bx-Ax)
Py = Ay + r(By-Ay)

And the distance from A to P = r*L.

Use another parameter s to indicate the location along PC, with the
following meaning:
s<0 C is left of AB
s>0 C is right of AB
s=0 C is on AB

Compute s as follows:

(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
s = -----------------------------
L^2


Then the distance from C to P = |s|*L.

*/


double r_numerator = (cx-ax)*(bx-ax) + (cy-ay)*(by-ay);
double r_denomenator = (bx-ax)*(bx-ax) + (by-ay)*(by-ay);
double r = r_numerator / r_denomenator;
//
double px = ax + r*(bx-ax);
double py = ay + r*(by-ay);
//
double s = ((ay-cy)*(bx-ax)-(ax-cx)*(by-ay) ) / r_denomenator;

distanceLine = fabs(s)*sqrt(r_denomenator);

//
// (xx,yy) is the point on the lineSegment closest to (cx,cy)
//
double xx = px;
double yy = py;

if ( (r >= 0) && (r <= 1) )
{
distanceSegment = distanceLine;
}
else
{

double dist1 = (cx-ax)*(cx-ax) + (cy-ay)*(cy-ay);
double dist2 = (cx-bx)*(cx-bx) + (cy-by)*(cy-by);
if (dist1 < dist2)
{
xx = ax;
yy = ay;
distanceSegment = sqrt(dist1);
}
else
{
xx = bx;
yy = by;
distanceSegment = sqrt(dist2);
}


}

return;
}

[Anthony Mai]

Obviously Philip Nicoletti's code is a bad example, since it used "if" and "else", No boolean logic should be involved in a simple calculation like determing the distance from a point to a line segment.

Suppose you have points A(xa, ya), B(xb, yb) and C(xc,yc). The distance between point C and line segment AB equals the area of parallelgram ABCC' divided by the length of AB.

So it can be written as simple as:
distance = |AB X AC| / sqrt(AB * AB)
Here X mean cross product of vectors, and * mean dot product of vectors. This applied in both 2 dimentional and three dimentioanl space.

In 2-D it becomes:
sqrt(((yb-ya)*(xc-xa)+(xb-xa)*(yc-ya))^2/((xb-xa)^2 + (yb-ya)^2))

[Philip Nicoletti]

Hi Anthony. Thanks (I think) for your comments on my code.
I am not really interested in getting into a discussion
about that, but I will mention something about that later.

I am more interested in the solution algorithm. Maybe I am
missing something, but I do not think your solution is correct.


The distance between point C and line segment AB equals the
area of parallelgram ABCC' divided by the length of AB.


If I am missing something, please let me know as I am interested
in these types of small mathematical problems. Or maybe I am
mis-interpreting your algorithm.

Some tests. Following is what my code calculates :

A(0,0) , B(1,0) , C(0.707,0.707) : distanceToLine = distanceToSegment = 0.707

A(0,0) , B(1,0) , C(-0.707,0.707) : distanceToLine = 0.707 , distanceToSegment = 1.0

A(0,0) , B(100,100) , C(200,-106) : distanceToLine = distanceToSegment = 216.37

======

As for my use of if statements in the solution.


notice in calculating the distance from the point to the
line not line segment, I did not use boolean
logic. I did code it in multiple lines instead of one. I did
this on purpose, so that the code would follow the given theory
more closely. It would be an easy exercise for someone interested
in the code to change this if desired.

Maybe there is an simple equation for finding the distance from
the point to the line segment, but I am not aware of it
(as I mentioned, I do not think your solution is correct). The
algorithm I used for the case of the line yields info on
whether the closest point is within the line segment or not. If
it is not within the line segment, then the closest point is one
of the end points of the segment. My code was designed to make that
clear.

the code does more than just calculated the two distances, it also
calculates the closest point to the line segment (although I
did forget to pass those as arguments).

This was one of the first functions I wrote when learning C++ (I
had Fortran code which did the same calculations). Although you did
not mention it, there are a couple of things I would do differently
now :

1) I would check the the end points of the line segment are not the same
(that is, the user really supplied a point only instead of a line segment).
In which case I would simply return the distance from C to A. I do not
really want to have a code that divides by zero !

2) the prototype and typical call for my function is :


void DistanceFromLine(double cx, double cy, double ax, double ay ,
double bx, double by, double &distanceSegment,
double &distanceLine);
//
//
//
DistanceFromLine(cx,cy,ax,ay,bx,by,distanceSegment,distanceLine);


The problem : from the typical call, there is no way to know that
distanceSegment and distanceLine are modified. So intead, maybe I
would do this instead :


void DistanceFromLine(double cx, double cy, double ax, double ay ,
double bx, double by, double *distanceSegment,
double *distanceLine);
//
//
//
DistanceFromLine(cx,cy,ax,ay,bx,by,&distanceSegment,&distanceLine);

[arah]

Thanks for the help! It proved three things.
1) That I was remembering how to calculate this correctly. While not the prettiest of code, my original routines produced the same results.
2) That my stuff wasn't working not because of the point to segment distance calculation, but because my line segment placement was thrown off-kilter by my not taking the window-space to OpenGL-space aspect ratios into consideration. (Stupid of me, I know. It was just a looooong day.)
3) That this forum is still full of cool and helpful people. :)

Again, thanks. Without your help I'd still be trying to figure out how to fix the wrong thing.

[arctanb]

Obviously Philip Nicoletti's code is a bad example, since it used "if" and "else", No boolean logic should be involved in a simple calculation like determing the distance from a point to a line segment.

Suppose you have points A(xa, ya), B(xb, yb) and C(xc,yc). The distance between point C and line segment AB equals the area of parallelgram ABCC' divided by the length of AB.

So it can be written as simple as:
distance = |AB X AC| / sqrt(AB * AB)
Here X mean cross product of vectors, and * mean dot product of vectors. This applied in both 2 dimentional and three dimentioanl space.

In 2-D it becomes:
sqrt(((yb-ya)*(xc-xa)+(xb-xa)*(yc-ya))^2/((xb-xa)^2 + (yb-ya)^2))


I believe your answer is incorrect - the area of the parallelogram is the det of the 2x2 matrix formed by the points: it should therefore be:

(yb-ya)*(xc-xa)-(xb-xa)*(yc-ya)

[nuzzle]

Obviously Philip Nicoletti's code is a bad example, since it used "if" and "else", No boolean logic should be involved in a simple calculation like determing the distance from a point to a line segment.


Finding the shortest distance to a line segment requires you to consider three distances; One to the line and two to the endpoints. To select the shortest of these requires logical statements because selection always does.

[JohnW@Wessex]

There is this equation,

http://mathworld.wolfram.com/images/equations/Point-LineDistance2-Dimensional/NumberedEquation8.gif


http://mathworld.wolfram.com/images/eps-gif/PointLineDistance2DVec_1000.gif


from this website. (http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html)

[JohnW@Wessex]

This may help too.

http://www.worsleyschool.net/science/files/linepoint/distance.html

[nuzzle]

This may help too.


I hope you know the difference between a line and a line segment?

[Richard.J]

can anyone explain why you are digging out a seven year old thread?

[arctanb]

can anyone explain why you are digging out a seven year old thread?

My bad - I was researching how to find the distance between a line segment and a point. I found this on google, implemented the ingenious suggestion of using areas and discovered it to be flawed so registered to post a correction in the spirit of making the web a better place...

[nuzzle]

My bad - I was researching how to find the distance between a line segment and a point. I found this on google, implemented the ingenious suggestion of using areas and discovered it to be flawed so registered to post a correction in the spirit of making the web a better place...

It's a good initiative.

Now if people only new the difference between a line and a line segment. :)

[arctanb]

Now if people only new the difference between a line and a line segment. :)

... and indeed a ray ;)

OK I think it's time for me to stop unwittingly bumping this topic

[nuzzle]

... and indeed a ray ;)

OK I think it's time for me to stop unwittingly bumping this topic

Well, the problem has an optimal algoritmic solution I think (see for example 3D Game Engine Design by Eberly).

But who knows, maybe someone comes up with something better. This happens too often to risk anything so please keep bumping the thread. :)

[blackhole12]

This thing shows up on google if you search for this, and I really hate getting 15 search results that give you "Here's an algorithm that's got your solution in it somewhere! If you can find it! Maybe!" The most important problem here is that while this equation here ALMOST works:

http://mathworld.wolfram.com/images/...dEquation8.gif

It actually treats your segment as a line, which is useless. I found the answer buried in some obscure c code that does a variation on this technique: http://local.wasp.uwa.edu.au/~pbourke/geometry/pointline/source.c - this contains a check that tells us if the point is within the line segment or not. Due to the similarity of the approaches i was able to integrate this into the previous equation for a function that, after almost 6 freaking hours of research, actually does exactly what we want it to do, under all circumstances.


float DistanceToPointSquared(float X, float Y)
{
float vx = x1-X, vy = y1-Y, ux = x2-x1, uy = y2-y1;
float length = ux*ux+uy*uy;

float det = (-vx*ux)+(-vy*uy); //if this is < 0 or > length then its outside the line segment
if(det<0 || det>length)
{
ux=x2-X;
uy=y2-Y;
return std::min<float>(vx*vx+vy*vy, ux*ux+uy*uy);
}

det = ux*vy-uy*vx;
return (det*det)/length;
}

[basiliskus]

keke, Necropost, but.

The code described above fails on the case point: 0,0. Line point0: 1,1, Line point1: 2,2.

The code:
public float DistanceSquared(float X, float Y)
{
float vx = v0.x - X, vy = v0.y - Y, ux = v1.x - v0.x, uy = v1.y - v0.y;
float length = ux * ux + uy * uy;

float det = (-vx * ux) + (-vy * uy); //if this is < 0 or > length then its outside the line segment
if (det < 0)
return (v0.x - X) * (v0.x - X) + (v0.y - Y) * (v0.y - Y);
if (det > length)
return (v1.x - X) * (v1.x - X) + (v1.y - Y) * (v1.y - Y);

det = ux * vy - uy * vx;
return (det * det) / length;
}

takes this into consideration. I think ( run a couple of tests, when point is outside the segment from both ends, above and below).

Also, v0 is first vertex of the segment. v1 is the second one.

[JohnW@Wessex]

This page mentions some useful info too.
Look for 'Line-Point Distance', second paragraph.

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=geometry1

[GCDEF]

This thread is over 9 years old.