language agnostic - How to calculate distance from a point to a line segment, on a sphere?

Question

I have a line segment (great circle part) on earth. The line segment is defined by the coordinates of its ends. Obviously, two points define two line segments, so assume I am interested in the shorter one.

I am given a third point, and I am looking for the (shortest) distance between the line and the point.

All the coordinates are given in longitudelatitude (WGS 84).

How do I calculate the distance?

A solution in any reasonable programming language will do.

Answer 1

Here's my own solution, based on the idea in ask Dr. Math. I'd be happy to see your feedback.

Disclaimer first. This solution is correct for spheres. Earth isn't a sphere, and the coordinates system (WGS 84) doesn't assume it's a sphere. So this is just an approximation, and I can't really estimate is error. Also, for very small distances, it's probably also possible to get good approximation by assuming everything is a just a coplanar. Again I don't know how "small" the distances have to be.

Now to business. I will call the ends of the lines A,B and the third point C. Basically, the algorithm is to:

1. convert the coordinates first to Cartesian coordinates (with the origin at the center of earth) - e.g. here.

2. Calculate T, the point on the line AB that is nearest to C, using the following 3 vector products:

   G = A x B

   F = C x G

   T = G x F

3. Normalize T and multiply by the radius of earth.

4. Convert T back to longitudelatitude.
5. Calculate the distance between T and C - e.g. here.

These steps are enough if you are looking for the distance between C and the great circle defined by A and B. If like me you are interested in the distance between C and the shorter line segment, you need to take the extra step of verifying that T is indeed on this segment. If it isn't, then necessarily the nearest point is one of the ends A or B - the easiest way is to check which one.

In general terms, the idea behind the three vector products is the following. The first one (G) gives us the the plane of the great circle of A and B (so the plane containing A,B and the origin). The second (F) gives us the great circle the goes through C and is perpendicular to the G. Then T is the intersection of the great circles defined by F and G, brought to the correct length by normalization and multiplication by R.

Here's some partial Java code for doing it.

Finding the nearest point on the great circle. The inputs and output are length-2 arrays. The intermediate arrays are of length 3.

double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
    double[] a_ = toCartsian(a);
    double[] b_ = toCartsian(b);
    double[] c_ = toCartsian(c);

    double[] G = vectorProduct(a_, b_);
    double[] F = vectorProduct(c_, G);
    double[] t = vectorProduct(G, F);
    normalize(t);
    multiplyByScalar(t, R_EARTH);
    return fromCartsian(t);
}

Finding the nearest point on the segment:

double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
   double[] t= nearestPointGreatCircle(a,b,c);
   if (onSegment(a,b,t))
     return t;
   return (distance(a,c) < distance(b,c)) ? a : c;
} 

This is a simple method of testing if point T, which we know is on the same great circle as A and B, is on the shorter segment of this great circle. However there are more efficient methods to do it:

   boolean onSegment (double[] a, double[] b, double[] t)
   {
     // should be   return distance(a,t)+distance(b,t)==distance(a,b), 
     // but due to rounding errors, we use: 
     return Math.abs(distance(a,b)-distance(a,t)-distance(b,t)) < PRECISION;
   }