Anonymous - 2022-09-08

Hello,
Im currently working on a project where it is necessary to find the intersection point between 2 lines with endpoints given in lat long coordinates. I have read some of the other posts in the forum, and seen an algorithm which should solve this - but I can not make it work when I translate it into C# - using the translated version of the GeographicLib (https://github.com/PyxisInt/GeographicLib)

Its the code provided below which is the problem. The first iteration, makes an ok guess, but it does not make a new guess affected by the previous resulting in just making the same calculations over and over. Can someone spot the error?

best regards - and thanks in advance! 

        double lata1 = 42.41593, lona1 = 3.17257, lata2 = 41.96081, lona2 = 4.65972;
        double latb1 = 42.207138, lonb1 = 3.766376, latb2 = 42.2975, lonb2 = 3.89165;

        Geodesic geod = new Geodesic(Constants.WGS84_a, Constants.WGS84_f);
        Gnomonic gn = new Gnomonic(geod);


        double lat0 = (lata1 + lata2 + latb1 + latb2) / 4;
        double lon0 = (lona1 + lona2 + lonb1 + lonb2) / 4;

        for (int i = 0; i < 10; i++)
        {
            double xa1, ya1, xa2, ya2;
            double xb1, yb1, xb2, yb2;

            var a = gn.Forward(lat0, lon0, lata1, lona1);
            var b = gn.Forward(lat0, lon0, lata2, lona2);
            var c = gn.Forward(lat0, lon0, latb1, lonb1);
            var d = gn.Forward(lat0, lon0, latb2, lonb2);

            xa1 = a.PointLatitude;
            ya1 = a.PointLongitude;
            xa2 = b.PointLatitude;
            ya2 = b.PointLongitude;
            xb1 = c.PointLatitude;
            yb1 = c.PointLongitude;
            xb2 = d.PointLatitude;
            yb2 = d.PointLongitude;

            // See Hartley and Zisserman, Multiple View Geometry, Sec. 2.2.1
            Vector3 va1 = new Vector3(xa1, ya1); Vector3 va2 = new Vector3(xa2, ya2);
            Vector3 vb1 = new Vector3(xb1, yb1); Vector3 vb2 = new Vector3(xb2, yb2);
            // la is homogeneous representation of line A1,A2
            // lb is homogeneous representation of line B1,B2
            Vector3 la = va1.Cross(va2);
            Vector3 lb = vb1.Cross(vb2);
            // p0 is homogeneous representation of intersection of la and lb
            Vector3 p0 = la.Cross(lb);
            p0.norm();

            double lat1, lon1;
            var q = gn.Reverse(lat0, lon0, p0._x, p0._y);
            lat1 = q.x;
            lon1 = q.y;
            lat0 = lat1;
            lon0 = lon1;

            Console.WriteLine(lat0.ToString() + ", " + lon0.ToString());
        }
        Console.WriteLine("Final:\n");
        Console.WriteLine(lat0.ToString() + ", " + lon0.ToString());
    }