libs/geometry/test/formulas/intersection.cpp

   1 // Boost.Geometry
   2 // Unit Test
   3
   4 // Copyright (c) 2016-2019 Oracle and/or its affiliates.
   5
   6 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
   7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   8
   9 // Use, modification and distribution is subject to the Boost Software License,
  10 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  11 // http://www.boost.org/LICENSE_1_0.txt)
  12
  13 #include "test_formula.hpp"
  14 #include "intersection_cases.hpp"
  15
  16 #include <boost/geometry/formulas/andoyer_inverse.hpp>
  17 #include <boost/geometry/formulas/geographic.hpp>
  18 #include <boost/geometry/formulas/gnomonic_intersection.hpp>
  19 #include <boost/geometry/formulas/sjoberg_intersection.hpp>
  20 #include <boost/geometry/formulas/thomas_direct.hpp>
  21 #include <boost/geometry/formulas/thomas_inverse.hpp>
  22 #include <boost/geometry/formulas/vincenty_direct.hpp>
  23 #include <boost/geometry/formulas/vincenty_inverse.hpp>
  24
  25 #include <boost/geometry/strategies/geographic/parameters.hpp>
  26
  27 #include <boost/geometry/srs/spheroid.hpp>
  28
  29 void check_result(expected_result const& result, expected_result const& expected,
  30                   expected_result const& reference, double reference_error,
  31                   bool check_reference_only)
  32 {
  33     //BOOST_CHECK_MESSAGE((false), "(" << result.lon << " " << result.lat << ") vs (" << expected.lon << " " << expected.lat << ")");
  34     check_one(result.lon, expected.lon, reference.lon, reference_error, false, check_reference_only);
  35     check_one(result.lat, expected.lat, reference.lat, reference_error, false, check_reference_only);
  36 }
  37
  38 void test_formulas(expected_results const& results, bool check_reference_only)
  39 {
  40     // reference result
  41     if (results.sjoberg_vincenty.lon == ND)
  42     {
  43         return;
  44     }
  45
  46     double const d2r = bg::math::d2r<double>();
  47     double const r2d = bg::math::r2d<double>();
  48
  49     double lona1r = results.p1.lon * d2r;
  50     double lata1r = results.p1.lat * d2r;
  51     double lona2r = results.p2.lon * d2r;
  52     double lata2r = results.p2.lat * d2r;
  53     double lonb1r = results.q1.lon * d2r;
  54     double latb1r = results.q1.lat * d2r;
  55     double lonb2r = results.q2.lon * d2r;
  56     double latb2r = results.q2.lat * d2r;
  57
  58     expected_result result;
  59
  60     // WGS84
  61     bg::srs::spheroid<double> spheroid(6378137.0, 6356752.3142451793);
  62
  63     if (results.gnomonic_vincenty.lon != ND)
  64     {
  65         bg::formula::gnomonic_intersection<double, bg::formula::vincenty_inverse, bg::formula::vincenty_direct>
  66             ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
  67         result.lon *= r2d;
  68         result.lat *= r2d;
  69         check_result(result, results.gnomonic_vincenty, results.sjoberg_vincenty, 0.00000001, check_reference_only);
  70     }
  71
  72     if (results.gnomonic_thomas.lon != ND)
  73     {
  74         bg::formula::gnomonic_intersection<double, bg::strategy::thomas::inverse, bg::strategy::thomas::direct>
  75                 ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
  76         result.lon *= r2d;
  77         result.lat *= r2d;
  78         check_result(result, results.gnomonic_thomas, results.sjoberg_vincenty, 0.0000001, check_reference_only);
  79     }
  80
  81     if (results.sjoberg_vincenty.lon != ND)
  82     {
  83         bg::formula::sjoberg_intersection<double, bg::formula::vincenty_inverse, 4>
  84             ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
  85         result.lon *= r2d;
  86         result.lat *= r2d;
  87         check_result(result, results.sjoberg_vincenty, results.sjoberg_vincenty, 0.00000001, check_reference_only);
  88     }
  89
  90     if (results.sjoberg_thomas.lon != ND)
  91     {
  92         bg::formula::sjoberg_intersection<double, bg::formula::thomas_inverse, 2>
  93             ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
  94         result.lon *= r2d;
  95         result.lat *= r2d;
  96         check_result(result, results.sjoberg_thomas, results.sjoberg_vincenty, 0.0000001, check_reference_only);
  97     }
  98
  99     if (results.sjoberg_andoyer.lon != ND)
 100     {
 101         bg::formula::sjoberg_intersection<double, bg::formula::andoyer_inverse, 1>
 102             ::apply(lona1r, lata1r, lona2r, lata2r, lonb1r, latb1r, lonb2r, latb2r, result.lon, result.lat, spheroid);
 103         result.lon *= r2d;
 104         result.lat *= r2d;
 105         check_result(result, results.sjoberg_andoyer, results.sjoberg_vincenty, 0.0001, check_reference_only);
 106     }
 107
 108     if (results.great_elliptic.lon != ND)
 109     {
 110         typedef bg::model::point<double, 2, bg::cs::geographic<bg::degree> > point_geo;
 111         typedef bg::model::point<double, 3, bg::cs::cartesian> point_3d;
 112         point_geo a1(results.p1.lon, results.p1.lat);
 113         point_geo a2(results.p2.lon, results.p2.lat);
 114         point_geo b1(results.q1.lon, results.q1.lat);
 115         point_geo b2(results.q2.lon, results.q2.lat);
 116         point_3d a1v = bg::formula::geo_to_cart3d<point_3d>(a1, spheroid);
 117         point_3d a2v = bg::formula::geo_to_cart3d<point_3d>(a2, spheroid);
 118         point_3d b1v = bg::formula::geo_to_cart3d<point_3d>(b1, spheroid);
 119         point_3d b2v = bg::formula::geo_to_cart3d<point_3d>(b2, spheroid);
 120         point_3d resv(0, 0, 0);
 121         point_geo res(0, 0);
 122         bg::formula::great_elliptic_intersection(a1v, a2v, b1v, b2v, resv, spheroid);
 123         res = bg::formula::cart3d_to_geo<point_geo>(resv, spheroid);
 124         result.lon = bg::get<0>(res);
 125         result.lat = bg::get<1>(res);
 126         check_result(result, results.great_elliptic, results.sjoberg_vincenty, 0.01, check_reference_only);
 127     }
 128 }
 129
 130 void test_4_input_combinations(expected_results const& results, bool check_reference_only)
 131 {
 132     test_formulas(results, check_reference_only);
 133
 134 #ifdef BOOST_GEOMETRY_TEST_GEO_INTERSECTION_TEST_SIMILAR
 135     {
 136         expected_results results_alt = results;
 137         std::swap(results_alt.p1, results_alt.p2);
 138         test_formulas(results_alt, true);
 139     }
 140     {
 141         expected_results results_alt = results;
 142         std::swap(results_alt.q1, results_alt.q2);
 143         test_formulas(results_alt, true);
 144     }
 145     {
 146         expected_results results_alt = results;
 147         std::swap(results_alt.p1, results_alt.p2);
 148         std::swap(results_alt.q1, results_alt.q2);
 149         test_formulas(results_alt, true);
 150     }
 151 #endif
 152 }
 153
 154 void test_all(expected_results const& results)
 155 {
 156     test_4_input_combinations(results, false);
 157
 158 #ifdef BOOST_GEOMETRY_TEST_GEO_INTERSECTION_TEST_SIMILAR
 159     expected_results results_alt = results;
 160     results_alt.p1.lat *= -1;
 161     results_alt.p2.lat *= -1;
 162     results_alt.q1.lat *= -1;
 163     results_alt.q2.lat *= -1;
 164     results_alt.gnomonic_vincenty.lat *= -1;
 165     results_alt.gnomonic_thomas.lat *= -1;
 166     results_alt.sjoberg_vincenty.lat *= -1;
 167     results_alt.sjoberg_thomas.lat *= -1;
 168     results_alt.sjoberg_andoyer.lat *= -1;
 169     results_alt.great_elliptic.lat *= -1;
 170     test_4_input_combinations(results_alt, true);
 171 #endif
 172 }
 173
 174 void test_bugs()
 175 {
 176     // https://github.com/boostorg/geometry/issues/612
 177     {
 178         double lon, lat;
 179         bg::formula::sjoberg_intersection<double, bg::formula::andoyer_inverse, 1>
 180             ::apply(-0.0872665, -0.0872665, -0.0872665, 0.0872665,
 181                     0.0, 1.57e-07, -0.392699, 1.57e-07,
 182                     lon, lat, bg::srs::spheroid<double>());
 183         check_one("issue 612", lon, -0.087266500535674751);
 184         check_one("issue 612", lat, 1.5892499139622920e-07);
 185     }
 186 }
 187
 188 int test_main(int, char*[])
 189 {
 190     for (size_t i = 0; i < expected_size; ++i)
 191     {
 192         test_all(expected[i]);
 193     }
 194
 195     test_bugs();
 196
 197     return 0;
 198 }