Find distance between line segments expressed using spherical coordinates

Question

I am trying to find Distance between two linesegments expressed in degrees(latitude, longitude). For example,

LINESTRING1 is ((4.42553, -124.159), (4.42553, -124.159)) and
LINESTRING2 is ((-0.061225, -128.428), (-0.059107, -128.428))

However, the usual way to find this distance is giving me incorrect distance because of the spherical dimensions of the earth. Is there some way by which I may improve the code below such that it returns me the correct distance which is: 558 km

#include <cstdlib>
#include <map>
#include <vector>
#include <algorithm>
#include <string>
#include <future>
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>

using namespace std;

double dot( double* V1, double* V2 )
{
    return V1[0]*V2[0] + V1[1]*V2[1];
}

void lineLineDist2()
        //Compute the square of the minimum distance between two line segments
{

        double seg1[4];
        seg1[0]=4.42553; seg1[1]=-124.159; seg1[2]=4.42553; seg1[3]=-124.159;
        double seg2[4];
        seg2[0]=-0.061225; seg2[1]=-128.428; seg2[2]=-0.059107; seg2[3]=-128.428;


        double u[2],v[2],w[2];

        //Vector   u = S1.P1 - S1.P0;
        u[0] = seg1[2] - seg1[0];
        u[1] = seg1[3] - seg1[1];

        //Vector   v = S2.P1 - S2.P0;
        v[0] = seg2[2] - seg2[0];
        v[1] = seg2[3] - seg2[1];

        //Vector   w = S1.P0 - S2.P0;
        w[0] = seg1[0] - seg2[0];
        w[1] = seg1[1] - seg2[1];

        double    a = dot(u,u);         // always >= 0
        double    b = dot(u,v);
        double    c = dot(v,v);         // always >= 0
        double    d = dot(u,w);
        double    e = dot(v,w);
        double    D = a*c - b*b;        // always >= 0
        double    sc, sN, sD = D;       // sc = sN / sD, default sD = D >= 0
        double    tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0                

        // compute the line parameters of the two closest points
        if (D < 0.00000001) {           // the lines are almost parallel
                sN = 0.0;               // force using point P0 on segment S1
                sD = 1.0;               // to prevent possible division by 0.0 later
                tN = e;
                tD = c;
        }
        else {                          // get the closest points on the infinite lines
                sN = (b*e - c*d);
                tN = (a*e - b*d);

                if (sN < 0.0) {         // sc < 0 => the s=0 edge is visible
                        sN = 0.0;
                        tN = e;
                        tD = c;
                }
                else if (sN > sD) {     // sc > 1  => the s=1 edge is visible
                        sN = sD;
                        tN = e + b;
                        tD = c;
                }
        }

        if (tN < 0.0) {                 // tc < 0 => the t=0 edge is visible

                tN = 0.0;

                // recompute sc for this edge
                if (-d < 0.0)
                        sN = 0.0;
                else if (-d > a)
                        sN = sD;
                else {
                        sN = -d;
                        sD = a;
                }
        }
        else if (tN > tD) {             // tc > 1  => the t=1 edge is visible

                tN = tD;

                // recompute sc for this edge
                if ((-d + b) < 0.0)
                        sN = 0;
                else if ((-d + b) > a)
                    sN = sD;
                else {
                        sN = (-d +  b);
                        sD = a;
                }
        }

        // finally do the division to get sc and tc
        sc = (abs(sN) < 0.00000001 ? 0.0 : sN / sD);
        tc = (abs(tN) < 0.00000001 ? 0.0 : tN / tD);

        // get the difference of the two closest points
        double dp[2];

        //Vector   dP = w + (sc * u) - (tc * v);  // =  S1(sc) - S2(tc)
        dp[0] = w[0] + sc *u[0] - tc*v[0];
        dp[1] = w[1] + sc *u[1] - tc*v[1];                

        cout<<"distance="<<dot(dp,dp)<<"\n";
}
//---------------------------------------------------------------------------
int main()
{
    lineLineDist2();
    return 0; 
}

It's more like a [Code Review](https://codereview.stackexchange.com) type of a question. It is not really related to C++ or C++11 — kazarey, Mar 04 '17 at 20:28
I'm voting to close this question as off-topic because it's a math question. — n. m. could be an AI, Mar 04 '17 at 20:29
See also http://stackoverflow.com/questions/10198985/calculating-the-distance-between-2-latitudes-and-longitudes-that-are-saved-in-a — kennytm, Mar 04 '17 at 20:29
@kennytm They have used haversine formula which works well for finding distance between points, but here I am trying to find distance between line segments — Jannat Arora, Mar 04 '17 at 20:32

score 2 · Answer 1 · answered Mar 05 '17 at 10:37

To find distance between two segments (big circle arcs) one has to check whether they are intersect ("Intersection of two paths" here).

If not, find the smallest cross-track distance from all segment ends to another segment (at the same page)

 dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
 where   δ13 is (angular) distance from start point to third point
 θ13 is (initial) bearing from start point to third point
 θ12 is (initial) bearing from start point to end point
 R is the earth’s radius

Find distance between line segments expressed using spherical coordinates

1 Answers1