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:

1 2 3 |
float DistancePtLine( Vec2 a, Vec2 b, Vec2 p ) { } |

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:

The problem of finding the distance of a point to a line makes use of finding the vector *d *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:

1 2 3 4 5 6 7 8 |
float DistancePtLine( Vec2 a, Vec2 b, Vec2 p ) { Vec2 n = b - a; Vec2 pa = a - p; Vec2 c = n * (Dot( pa, n ) / Dot( n, n )); Vec2 d = pa - c; return sqrt( Dot( d, d ) ); } |

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:

1 2 3 4 5 6 7 8 |
float SqDistancePtLine( Vec2 a, Vec2 b, Vec2 p ) { Vec2 ab = b - a, ap = p - a, bp = p - b; float e = Dot( ap, ab ); return Dot( ap, ap ) - e * e / Dot( ab, ab ); } |

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 *NaN*s 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:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
float SqDistancePtSegment( Vec2 a, Vec2 b, Vec2 p ) { Vec2 n = b - a; Vec2 pa = a - p; float c = Dot( n, pa ); // Closest point is a if ( c > 0.0f ) return Dot( pa, pa ); Vec2 bp = p - b; // Closest point is b if ( Dot( n, bp ) > 0.0f ) return Dot( bp, bp ); // Closest point is between a and b Vec2 e = pa - n * (c / Dot( n, n )); return Dot( e, e ); } |

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:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 |
#include <cstdio> struct Vec2 { float x; float y; const Vec2 operator-( const Vec2& a ) { Vec2 v; v.x = x - a.x; v.y = y - a.y; return v; } const Vec2 operator*( float a ) const { Vec2 v; v.x = x * a; v.y = y * a; return v; } }; float Dot( const Vec2& a, const Vec2& b ) { return a.x * b.x + a.y * b.y; } int main( ) { Vec2 a, b, p; a.x = 1.0f; a.y = 1.0f; b.x = 5.0f; b.y = 2.0f; p.x = 3.0f; p.y = 3.0f; Vec2 n = b - a; Vec2 pa = a - p; Vec2 c = n * (Dot( n, pa ) / Dot( n, n )); Vec2 d = pa - c; float d2 = Dot( d, d ); printf( "Distance squared: %f\n", d2 ); } |

The output is: “Distance squared: 2.117647”.

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

satheeshHi,

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;

}

Randy GaulPost authorSorry, 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.

Aaron KrajeskiI think you forgot to normalize the vector

nin the code above.main()returns a square magnitude of 2.Randy GaulPost authorA 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.

Aaron KrajeskiIf I run the code in the

Sample Programthe output is “Distance squared: 2” not “Distance squared: 2.117647.” I believe this is becausenis meant to be a unit vector, according to the wikipedia article, whereasb – ais of lengthsqrt(17). Unless I’m missing something incredibly obvious, which is also possible.Aaron KrajeskiWait, 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

amconfused about the fact thatnis not of unit length.Aaron KrajeskiGah! 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

nisn’t normalized, but I’ll solve that on my own time.Randy GaulPost authorOh 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.

Aaron KrajeskiGotcha, it was the fact that the magnitude of

ngets squared in the numerator that I missed, so you can simply divide bydot(Thanks for your prompt replies and super blog post!n,n).Randy GaulPost authordot( 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 :)

FernandoFirst 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

pinstead of returning a scalar. The idea is to use the vector frompto the closest position onabto projectponto the segmentab.”How would that be?

Thanks

Randy GaulPost authorThink 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.