Improve short-distance accuracy of Andoyer inverse

[mjacobse]

Use haversine formula instead of law of cosines to avoid numerical loss of precision for close points as suggested in #1217.

This has the side-effect of returning non-zero azimuths for smaller angles than before, which results in quite inaccurate non-zero reduced length and geodesic scale values too. To fix this, division by the cosine of the latitudes is changed to be done implicitly within atan2, which together with use of the haversine formula improves accuracy for all result values.

I planned to add tests for very close points like {8.81, 53.08}, {8.81000001, 53.08} to inverse_cases.hpp, since Andoyer returns a distance of about 9cm for such examples before this change, which has a very large relative error. Unfortunately the other formulas (especially Thomas) have some issues with examples like that too, so I could not easily get such tests to pass. Suggestions would be welcome.

The following plots (generated using a simple test executable and a Python script, see mjacobse@d2930ae) show the errors of the different quantities compared to GeographicLib for many randomly generated pairs of points depending on the order of magnitude of the distance between them (top row absolute error, bottom row relative error. red before these changes, green after):

distance:

azimuth:

reverse_azimuth:

reduced_length:

geodesic_scale:

[tinko92]

Thanks for this PR. I find your test method and results convincing regarding the numeric improvements.

The CI failures are unrelated to the code changes, I believe (saw them in other runs as well). Running the test suite manually on an x86_64 (GCC14) machine with this PR, I get new test failures, e.g. one in test/algorithms/area/area_sph_geo.cpp and two failures in test/algorithms/area/area.cpp . Can you look into those? You can run the full test suite with b2 test from the geometry main directory, if you have a setup as discussed at https://github.com/boostorg/geometry/wiki/Contribution-Tutorial

Regarding performance, with a naive benchmark based on your testing code (calling only the inverse function in a loop, separated from the random generation of the test cases), I find a performance cost of ~14%. To me this seems acceptable for the improved precision but I don't know the typical use cases, maybe @vissarion can comment on this (I remember you did work in the past on the tradeoffs of different geographic distance formulas).

[tinko92]

(Minor) These names and comments may be less obvious now with the cos_d identifier no longer being defined in this file (meaning can still be inferred from the other _d identifiers though and I don't have proposals for other names).

[mjacobse]

Maybe adding a comment like

// follows from definition of haversine and trigonometric power reduction formula:
// hav(d) = sin^2(d/2) = (1 - cos(d))/2

could help?

[barendgehrels]

Also here the * could use spaces for readability

[mjacobse]

Running the test suite manually on an x86_64 (GCC14) machine with this PR, I get new test failures, e.g. one in test/algorithms/area/area_sph_geo.cpp and two failures in test/algorithms/area/area.cpp . Can you look into those? You can run the full test suite with b2 test from the geometry main directory, if you have a setup as discussed at https://github.com/boostorg/geometry/wiki/Contribution-Tutorial

Sure, sorry for not checking this before.

The failure in area_sph_geo.cpp only seems to be due to the area using andoyer being more accurate now (closer to thomas, vincenty, karney and GeographicLib), so I updated the reference value for andoyer.

The failures in area.cpp seem to be caused by a difference in behavior for out-of-range latitudes. The test uses (lon=45°, lat=95°) and (lon=45°, lat=-95°) as inputs. Using the what I believe should be the equivalent in-bounds coordinates (lon=-135°, lat=85°) and (lon=-135°, lat=-85°) instead makes develop and this branch behave the same way on the test, although they both fail and I am not quite sure why, I would have to dig deeper into the area computation.

Does the library guarantee some behavior for such out-of-range values or are such calls out of contract? Even if it is out of contract, should I try to restore the same behavior for out-of-range coordinates with the haversine formula to before to not break users who rely on this behavior?

[barendgehrels]

hi @mjacobse , could you maybe add this testcase?

// Short-distance accuracy test for andoyer_inverse.
// Reference for the regression: https://github.com/boostorg/geometry/issues/1217
// Fix proposed in PR #1461.
//
// On develop the issue_1217 case fails: andoyer returns ~9 cm for a true
// distance of ~0.67 mm. After the PR, andoyer agrees with vincenty to
// well under a millimetre across the whole range tested here.

#include <geometry_test_common.hpp>

#include <boost/geometry/formulas/andoyer_inverse.hpp>
#include <boost/geometry/formulas/vincenty_inverse.hpp>
#include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/util/math.hpp>

namespace
{

// Maximum allowed disagreement between andoyer and vincenty for the cases
// below. 1 mm comfortably catches the ~9 cm regression from issue #1217 and
// the cancellation-to-zero failures, while leaving headroom for the
// first-order Andoyer-Lambert approximation error at the metre scale.
double const tolerance_m = 1.0e-3;

struct short_case
{
    char const* name;
    double lon1, lat1, lon2, lat2;
};

void test_short(short_case const& c, double abs_tol_meters = tolerance_m)
{
    double const d2r = bg::math::d2r<double>();
    bg::srs::spheroid<double> const spheroid; // WGS84 by default

    using vinc_t = bg::formula::vincenty_inverse<double, true, false, false, false, false>;
    using ando_t = bg::formula::andoyer_inverse <double, true, false, false, false, false>;

    double const lon1r = c.lon1 * d2r;
    double const lat1r = c.lat1 * d2r;
    double const lon2r = c.lon2 * d2r;
    double const lat2r = c.lat2 * d2r;

    double const ref = vinc_t::apply(lon1r, lat1r, lon2r, lat2r, spheroid).distance;
    double const got = ando_t::apply(lon1r, lat1r, lon2r, lat2r, spheroid).distance;

    double const diff = bg::math::abs(got - ref);

    BOOST_CHECK_MESSAGE(diff <= abs_tol_meters,
        c.name
            << ": vincenty=" << ref << " m"
            << ", andoyer="  << got << " m"
            << ", diff="     << diff << " m"
            << ", tol="      << abs_tol_meters << " m");
}

} // namespace

int test_main(int, char*[])
{
    // Marquee case from issue #1217: ~0.67 mm true distance at mid-latitude.
    // Develop returns ~9 cm here (the regression the PR fixes).
    test_short({"issue_1217",      8.81, 53.08, 8.81000001, 53.08      });

    // East-west steps at the equator, sweeping sub-mm to ~10 m.
    test_short({"sub_mm_equator",  0.0,   0.0,  1.0e-8,      0.0       });
    test_short({"cm_equator",      0.0,   0.0,  1.0e-7,      0.0       });
    test_short({"m_equator",       0.0,   0.0,  1.0e-5,      0.0       });
    test_short({"10m_equator",     0.0,   0.0,  1.0e-4,      0.0       });

    // North-south steps along a meridian.
    test_short({"sub_mm_meridian", 0.0,  45.0,  0.0,         45.0 + 1.0e-8 });
    test_short({"m_meridian",      0.0,  45.0,  0.0,         45.0 + 1.0e-5 });

    // High-latitude oblique step (cos(lat) is small, so longitude differences shrink).
    test_short({"oblique_70N",     10.0, 70.0, 10.0 + 1.0e-7, 70.0 + 1.0e-7 });

    return 0;
}

in a new file andoyer_inverse_short.cpp (for "short distances")? It currently fails (on develop), but I expect that most, or maybe all, of the test cases improve with your PR.

It should then also be added to ...../test/formulas/CMakeLists.txt

It compares with Vincenty, that has a higher precision (we cannot compare with geographicLib without adding hardcoded expectations)

The cases themselves are suggested by Claude, feel free to adapt it - or adapt the tolerance

Current behaviour:

Running 1 test case...
error: in "test_main_caller( argc_ argv )": issue_1217: vincenty=0.00066931725477152083 m, andoyer=0.095245320237387923 m, diff=0.094576002982616397 m, tol=0.001 m
error: in "test_main_caller( argc_ argv )": sub_mm_equator: vincenty=0.0011094625761734762 m, andoyer=0 m, diff=0.0011094625761734762 m, tol=0.001 m
error: in "test_main_caller( argc_ argv )": cm_equator: vincenty=0.011131823941310088 m, andoyer=0 m, diff=0.011131823941310088 m, tol=0.001 m
error: in "test_main_caller( argc_ argv )": sub_mm_meridian: vincenty=0.0011113174432166772 m, andoyer=0 m, diff=0.0011113174432166772 m, tol=0.001 m
error: in "test_main_caller( argc_ argv )": m_meridian: vincenty=1.111317774736599 m, andoyer=1.1146369167819898 m, diff=0.0033191420453908549 m, tol=0.001 m
error: in "test_main_caller( argc_ argv )": oblique_70N: vincenty=0.011791161304736061 m, andoyer=0 m, diff=0.011791161304736061 m, tol=0.001 m

I didn't test on your PR