Distance Point to Line Segment

Understanding the dot product will enable one to come up with an algorithm to compute the distance between a point and line in two or three dimensions without much fuss. This is a math problem but discussion and focus is on writing a function to compute distance.

Dot Product Intuition

The dot product comes from the law of cosines. Here’s the formula:

\begin{equation}
c^2 = a^2 + b^2 – 2ab * cos\;\gamma
\label{eq1}
\end{equation}

This is just an equation that relates the cosine of an angle within a triangle to its various side lengths a, b and c. The Wikipedia page (link above) does a nice job of explaining this. Equation \eqref{eq1} can be rewritten as:

\begin{equation}
c^2 – a^2 – b^2 = -2ab * cos\;\gamma
\label{eq2}
\end{equation}

The right hand side equation \eqref{eq2} is interesting! Lets say that instead of writing the equation with side lengths a, b and c, it is written with two vectors: u and v. The third side can be represented as u – v. Re-writing equation \eqref{eq2} in vector notation yields:

\begin{equation}
|u\;-\;v|^2 – |u|^2 – |v|^2 = -2|u||v| * cos\;\gamma
\label{eq3}
\end{equation}

Which can be expressed in scalar form as:

\begin{equation}
(u_x\;-\;v_x)^2 + (u_y\;-\;v_y)^2 + (u_z\;-\;v_z)^2\;- \\
(u_{x}^2\;+\;u_{y}^2\;+\;u_{z}^2) – (v_{x}^2\;+\;v_{y}^2\;+\;v_{z}^2)\;= \\
-2|u||v| * cos\;\gamma
\label{eq4}
\end{equation}

By crossing out some redundant terms, and getting rid of the -2 on each side of the equation, this ugly equation can be turned into a much more approachable version:

\begin{equation}
u_x v_x + u_y v_y + u_w v_w = |u||v| * cos\;\gamma
\label{eq5}
\end{equation}

Equation \eqref{eq5} is the equation for the dot product. If both u and v are unit vectors then the equation will simplify to:

\begin{equation}
dot(\;\hat{u},\;\hat{v}\;) = cos\;\gamma
\label{eq6}
\end{equation}

If u and v are not unit vectors equation \eqref{eq5} says that the dot product between both vectors is equal to cos( γ ) that has been scaled by the lengths of u and v. This is a nice thing to know! For example: the squared length of a vector is just itself dotted with itself.

If u is a unit vector and v is not, then dot( u, v ) will return the distance in which v travels in the u direction. This is useful for understanding the plane equation in three dimensions (or any other dimension):

\begin{equation}
ax\;+\;by\;+\;cz\;-\;d\;=\;0
\end{equation}

The normal of a plane would be the vector: { a, b, c }. If this normal is a unit vector, then d represents the distance to the plane from the origin. If the normal is not a unit vector then d is scaled by the length of the normal.

To compute the distance of a point to this plane any point can be substituted into the plane equation, assuming the normal of the plane equation is of unit length. This operation is computing the distance along the normal a given point travels. The subtraction by d can be viewed as “translating the plane to the origin” in order to convert the distance along the normal, to a distance to the plane.

Writing the Function: Distance Point to Line

The simplest function for computing the distance to a plane (or line in 2D) would be to place a point into the plane equation. This means that we’ll have to either compute the plane equation in 2D if all we have are two points to represent the plane, and in 3D find a new tactic altogether since planes in 3D are not lines.

In my own experience I’ve found it most common to have a line in the form of two points in order to represent the parametric equation of a line. Two points can come from a triangle, a mesh edge, or two pieces of world geometry.

To setup the problem lets outline the function to be created as so:

The two parameters a and b are used to define the line segment itself. The direction of the line would be the vector b – a.

After a brief visit to the Wikipedia page for this exact problem I quickly wrote down my own derivation of the formula they have on their page. Take a look at this image I drew:

2014-07-23 00.16.37

The problem of finding the distance of a point to a line makes use of finding the vector that points from p to the closest point on the line ab. From the above picture: a simple way to calculate this vector would be to subtract away the portion of a – p that travels along the vector ab.

The part of a – p that travels along ab can be calculated by projecting a – p onto ab. This projection is described in the previous section about the dot product intuition.

Given the vector d the distance from p to ab is just sqrt( dot( d, d ) ). The sqrt operation can be omit entirely to compute a distance squared. Our function may now look like:

This function is quite nice because it will never return a negative number. There is a popular version of this function that performs a division operation. Given a very small line segment as input for ab it is entirely possible to have the following function return a negative number:

It’s very misleading to have a function called “square distance” or “distance” to return a negative number. Passing in the result of this function to a sqrt function call can result in NaNs and be really nasty to deal with.

Barycentric Coordinates – Segments

A full discussion of barycentric coordinates is way out of scope here. However, they can be used to compute distance from a point to line segment. The segment portion of the code just clamps a point projected into the line within the bounds of a and b.

Assuming readers are a little more comfortable with the dot product than I was when I first started programming, the following function should make sense:

This function can be adapted pretty easily to compute the closest point on the line segment to p instead of returning a scalar. The idea is to use the vector from p to the closest position on ab to project p onto the segment ab.

The above function works by computing barycentric coordinates of p relative to ab. The coordinates are scaled by the length of ab so the second if statement must be adapted slightly. If the direction ab were normalized then the second if statement would be a comparison with the value 1, which should make sense for barycentric coordinates.

Sample Program

Here’s a sample program you can try out yourself:

The output is: “Distance squared: 2.117647”.

TwitterRedditFacebookShare

13 thoughts on “Distance Point to Line Segment

  1. Pingback: Inertia Tensor – Capsule | Randy Gaul's Game Programming Blog

  2. satheesh

    Hi,

    I am working in fluent(CFD software) and performed a tutorial that involves moving mesh(grid points) there is a c code to move the mesh and it was written some 15 years ago and i do not have the explanation of code and i am trying to understand it and i am struck at this place.

    My problem is purely mathematical (i am poor in maths), this part of the code computes a value t1 and t2.Then when some criteria is satisfied t1 is used to update the mesh they dont use t2.

    Now i cannot figure out geometrically what t1 and t2 means and totally struck.

    I think t1 and t2 parameters of barycentric coordinates(not sure) and really searching over the net hoplessly, could you please help me graphically understand what these t1 and t2 mean.

    (Xa,Ya)-coordinate of first node

    (Xb,Yb)-coordinate of second node

    (Xc,Yc)-origin about which the mesh rotates(0,0)

    (xd,yd)-The node selected to get the coordinate updated

    i am attaching the part of the code,

    static void find_t(real xa, real ya, real xb, real yb,real xc, real yc, real xd, real yd,
    real *t1, real *t2, int * flag)

    {
    if((xd==xc)&&(xb!=xa))

    {
    if(yd==yc)

    {
    Message0(“\nc and d are the same points-aborting!!!\n”);

    exit(0);

    }
    (*flag) = 2;

    (*t1) = (xc-xa)/(xb-xa);

    (*t2) = ((yb-ya)*(*t1)-(yc-ya))/(yd-yc);

    }

    else if((xd!=xc)&&(xb==xa))

    {

    if(yb==ya)

    {
    Message0(“\na and b are the same points-aborting!!!\n”);

    exit(0);

    }
    (*flag) = 2;

    (*t2) = (xc-xa)/(xc-xd);

    (*t1) = ((yd-yc)*(*t2)+(yc-ya))/(yb-ya);

    }

    else if((xd==xc)&&(xb==xa))

    {
    (*flag) = 0;
    return;

    }

    else

    {
    real k1, k2;
    k1 = (yd-yc)/(xd-xc);

    k2 = (yb-ya)/(xb-xa);
    if(fabs(k2)>=fabs(k1))

    {

    if(k1==k2)

    {

    (*flag) = 0;

    return;

    }
    (*flag) = 1;

    (*t1) = (-k1*(xc-xa)+(yc-ya))/((yb-ya)-k1*(xb-xa));

    (*t2) = ((xb-xa)*(*t1)-(xc-xa))/(xd-xc);

    }

    else

    {
    (*flag) = 1;

    (*t2) = (-k2*(xc-xa)+(yc-ya))/(-(yd-yc)+k2*(xd-xc));

    (*t1) = ((xd-xc)*(*t2)+(xc-xa))/(xb-xa);

    }
    }
    return;

    }

     

     

    then this function returns the value to this part of program,

    if(((flag==1)||(flag==2))&&(t1<=1+cal_small_number(j))&&(t1>=-cal_small_number(j))&&(t2>=0))

    {

    (*x) = (1-t1)*pt1[0] + t1*pt2[0];

    (*y) = (1-t1)*pt1[1] + t1*pt2[1];
    return;

    }

    Reply
    1. Randy Gaul Post author

      Sorry, I don’t have time to look into anyone’s code, especially math related code. Best of luck to you! Feel free to email me specific questions if you have them.

      Reply
    1. Randy Gaul Post author

      A lot of time distance routines return a squared distance in order to avoid the sqrt operation. So, like you pointed out, if the real distance is needed a normalization would occur.

      Reply
      1. Aaron Krajeski

        If I run the code in the Sample Program the output is “Distance squared: 2” not “Distance squared: 2.117647.” I believe this is because n is meant to be a unit vector, according to the wikipedia article, whereas b – a is of length sqrt(17). Unless I’m missing something incredibly obvious, which is also possible.

        Reply
        1. Aaron Krajeski

          Wait, I must be doing something wrong, I just ran your code as is and go the correct answer. I’ve been running a python version this whole time…

          Regardless, I am confused about the fact that n is not of unit length.

          Reply
          1. Aaron Krajeski

            Gah! Ints v. floats was my problem, of course. Sorry to have wasted your time. Thank you so much for the excellent post. Again, still confused about why n isn’t normalized, but I’ll solve that on my own time.

          2. Randy Gaul Post author

            Oh I’m sorry. I answered a different question than the one you were asking! Yes, like you pointed out in the sample code n is not normalized. If you take a close look at the definition of the dot product it is: dot( a, b ) = cos(theta) * |a| * |b|. So, if a or b are not normalized the dot product is scaled by their lengths. This is fine since the demo code does: divide by dot( n, n ). This division here is used instead of a sqrt operation. Sometimes this can be preferred (choosing a division over a square root). In a sense this division is accounting for the length of n. Alternatively the division can be removed if the length of n is 1.

            The choice between the two is mostly a matter of preference, and can sometimes be a matter of efficiency. If we normalize n we must compute a dot product and one sqrt. If we divide, we must compute a dot product and one division. Is division cheaper than sqrt? Probably! The point is by carefully crafting the code some operations can be avoided or changed without compromising the result. To do this we have to understand the dot product definition very well.

          3. Aaron Krajeski

            Gotcha, it was the fact that the magnitude of n gets squared in the numerator that I missed, so you can simply divide by dot(n,n). Thanks for your prompt replies and super blog post!

          4. Randy Gaul Post author

            dot( d, d ) is equal to: xx + yy, where [x, y] is d. We know that the distance between two points A and B is sqrt( dot( B – A, B – A ) ). B – A is a subtraction of two points, and a subtraction of two points results in a vector, say, vector d. Therefore we substitute A – B with d and get: sqrt( dot( d, d ) ). Here is where often times the sqrt can be ignored, for example:

            Is equivalent to:

            Does that make sense? Hope that helps!

            Edit: OOPS. This seems not to be what you were looking for, please see my other answer :)

  3. Fernando

    First of all thanks for the great article.

    Second, can you please explain a little further how to get the closest point?

     

    You mentioned that “This function can be adapted pretty easily to compute the closest point on the line segment to p instead of returning a scalar. The idea is to use the vector from p to the closest position on ab to project p onto the segment ab.”

     

    How would that be?

     

    Thanks

    Reply
    1. Randy Gaul Post author

      Think of the segment AB as lining up along the x axis, and A is on the left. Given P, first see if P is to the left, or to the right of A. This can be done by testing the sign of dot( A – P, P ). If P is to the left of A, the closest point is A. Otherwise see if P is to the right of B by checking the sign of dot( B – P, P ). If P is to the left of B, project P onto the plane of AB (or in other words, remove the y value), otherwise the closest point is B.

      To handle the case of the AB not lining up nicely on the x axis you will need to project P onto an arbitrary plane. You will also need to figure out how to do the sign checks correctly.

      Reply

Leave a Reply

Your email address will not be published. Required fields are marked *