ceph/src/boost/libs/geometry/test/strategies/andoyer.cpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // Unit Test
   3
   4 // Copyright (c) 2007-2016 Barend Gehrels, Amsterdam, the Netherlands.
   5 // Copyright (c) 2008-2016 Bruno Lalande, Paris, France.
   6 // Copyright (c) 2009-2016 Mateusz Loskot, London, UK.
   7
   8 // This file was modified by Oracle on 2014-2017.
   9 // Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
  12 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  13
  14 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  15 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  16
  17 // Use, modification and distribution is subject to the Boost Software License,
  18 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  19 // http://www.boost.org/LICENSE_1_0.txt)
  20
  21
  22 #include <geometry_test_common.hpp>
  23
  24 #include <boost/concept_check.hpp>
  25
  26 #include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
  27 #include <boost/geometry/strategies/geographic/side_andoyer.hpp>
  28
  29 #include <boost/geometry/core/srs.hpp>
  30 #include <boost/geometry/strategies/strategies.hpp>
  31 #include <boost/geometry/algorithms/assign.hpp>
  32 #include <boost/geometry/geometries/point.hpp>
  33 #include <test_common/test_point.hpp>
  34
  35 #ifdef HAVE_TTMATH
  36 #  include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
  37 #endif
  38
  39
  40
  41 double make_deg(double deg, double min, double sec)
  42 {
  43     return deg + min / 60.0 + sec / 3600.0;
  44 }
  45
  46 double to_rad(double deg)
  47 {
  48     return bg::math::pi<double>() * deg / 180.0;
  49 }
  50
  51 double to_deg(double rad)
  52 {
  53     return 180.0 * rad / bg::math::pi<double>();
  54 }
  55
  56 double normlized_deg(double deg)
  57 {
  58     if (deg > 180)
  59         return deg - 360;
  60     else if (deg < -180)
  61         return deg + 360;
  62     else
  63         return deg;
  64 }
  65
  66
  67 template <typename P1, typename P2>
  68 void test_distance(double lon1, double lat1, double lon2, double lat2, double expected_km)
  69 {
  70     // Set radius type, but for integer coordinates we want to have floating point radius type
  71     typedef typename bg::promote_floating_point
  72         <
  73             typename bg::coordinate_type<P1>::type
  74         >::type rtype;
  75
  76     typedef bg::srs::spheroid<rtype> stype;
  77
  78     typedef bg::strategy::distance::andoyer<stype> andoyer_type;
  79     typedef bg::formula::andoyer_inverse<rtype, true, false> andoyer_inverse_type;
  80
  81     BOOST_CONCEPT_ASSERT
  82         (
  83             (bg::concepts::PointDistanceStrategy<andoyer_type, P1, P2>)
  84         );
  85
  86     andoyer_type andoyer;
  87     typedef typename bg::strategy::distance
  88         ::services::return_type<andoyer_type, P1, P2>::type return_type;
  89
  90     P1 p1;
  91     P2 p2;
  92
  93     bg::assign_values(p1, lon1, lat1);
  94     bg::assign_values(p2, lon2, lat2);
  95
  96     return_type d_strategy = andoyer.apply(p1, p2);
  97     return_type d_function = bg::distance(p1, p2, andoyer);
  98
  99     double diff = bg::math::longitude_distance_signed<bg::degree>(lon1, lon2);
 100     return_type d_formula;
 101
 102     // if the points lay on a meridian, distance strategy calls the special formula
 103     // for meridian distance that returns different result than andoyer formula
 104     // for nearly antipodal points
 105     if (bg::math::equals(diff, 0.0)
 106        || bg::math::equals(bg::math::abs(diff), 180.0))
 107     {
 108         d_formula = d_strategy;
 109     }
 110     else
 111     {
 112         d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1),
 113                                                 to_rad(lon2), to_rad(lat2),
 114                                                 stype()).distance;
 115     }
 116
 117     BOOST_CHECK_CLOSE(d_strategy / 1000.0, expected_km, 0.001);
 118     BOOST_CHECK_CLOSE(d_function / 1000.0, expected_km, 0.001);
 119     BOOST_CHECK_CLOSE(d_formula / 1000.0, expected_km, 0.001);
 120 }
 121
 122 template <typename PS, typename P>
 123 void test_azimuth(double lon1, double lat1,
 124                   double lon2, double lat2,
 125                   double expected_azimuth_deg)
 126 {
 127     // Set radius type, but for integer coordinates we want to have floating point radius type
 128     typedef typename bg::promote_floating_point
 129         <
 130             typename bg::coordinate_type<PS>::type
 131         >::type rtype;
 132
 133     typedef bg::srs::spheroid<rtype> stype;
 134     typedef bg::formula::andoyer_inverse<rtype, false, true> andoyer_inverse_type;
 135
 136     rtype a_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1), to_rad(lon2), to_rad(lat2), stype()).azimuth;
 137
 138     rtype azimuth_deg = to_deg(a_formula);
 139
 140     if (bg::math::equals(azimuth_deg, -180.0))
 141         azimuth_deg = 180.0;
 142     if (bg::math::equals(expected_azimuth_deg, -180.0))
 143         expected_azimuth_deg = 180.0;
 144
 145     if (bg::math::equals(expected_azimuth_deg, 0.0))
 146     {
 147         BOOST_CHECK(bg::math::equals(azimuth_deg, expected_azimuth_deg));
 148     }
 149     else
 150     {
 151         BOOST_CHECK_CLOSE(azimuth_deg, expected_azimuth_deg, 0.001);
 152     }
 153 }
 154
 155 template <typename P1, typename P2>
 156 void test_distazi(double lon1, double lat1, double lon2, double lat2,
 157                   double expected_km, double expected_azimuth_deg)
 158 {
 159     test_distance<P1, P2>(lon1, lat1, lon2, lat2, expected_km);
 160     test_azimuth<P1, P2>(lon1, lat1, lon2, lat2, expected_azimuth_deg);
 161 }
 162
 163 // requires SW->NE
 164 template <typename P1, typename P2>
 165 void test_distazi_symm(double lon1, double lat1, double lon2, double lat2,
 166                        double expected_km, double expected_azimuth_deg,
 167                        bool is_antipodal = false)
 168 {
 169     double d180 = is_antipodal ? 0 : 180;
 170     test_distazi<P1, P2>(lon1, lat1, lon2, lat2, expected_km, expected_azimuth_deg);
 171     test_distazi<P1, P2>(-lon1, lat1, -lon2, lat2, expected_km, -expected_azimuth_deg);
 172     test_distazi<P1, P2>(lon1, -lat1, lon2, -lat2, expected_km, d180 - expected_azimuth_deg);
 173     test_distazi<P1, P2>(-lon1, -lat1, -lon2, -lat2, expected_km, -d180 + expected_azimuth_deg);
 174 }
 175
 176 template <typename P1, typename P2>
 177 void test_distazi_symmNS(double lon1, double lat1, double lon2, double lat2,
 178                          double expected_km, double expected_azimuth_deg)
 179 {
 180     test_distazi<P1, P2>(lon1, lat1, lon2, lat2, expected_km, expected_azimuth_deg);
 181     test_distazi<P1, P2>(lon1, -lat1, lon2, -lat2, expected_km, 180 - expected_azimuth_deg);
 182 }
 183
 184
 185 template <typename PS, typename P>
 186 void test_side(double lon1, double lat1,
 187                double lon2, double lat2,
 188                double lon, double lat,
 189                int expected_side)
 190 {
 191     // Set radius type, but for integer coordinates we want to have floating point radius type
 192     typedef typename bg::promote_floating_point
 193         <
 194             typename bg::coordinate_type<PS>::type
 195         >::type rtype;
 196
 197     typedef bg::srs::spheroid<rtype> stype;
 198
 199     typedef bg::strategy::side::andoyer<stype> strategy_type;
 200
 201     strategy_type strategy;
 202
 203     PS p1, p2;
 204     P p;
 205
 206     bg::assign_values(p1, lon1, lat1);
 207     bg::assign_values(p2, lon2, lat2);
 208     bg::assign_values(p, lon, lat);
 209
 210     int side = strategy.apply(p1, p2, p);
 211
 212     BOOST_CHECK_EQUAL(side, expected_side);
 213 }
 214
 215 template <typename P1, typename P2>
 216 void test_all()
 217 {
 218     // polar
 219     test_distazi<P1, P2>(0, 90, 1, 80,
 220                          1116.814237, 179);
 221
 222     // no point difference
 223     test_distazi<P1, P2>(4, 52, 4, 52,
 224                          0.0, 0.0);
 225
 226     // normal cases
 227     test_distazi<P1, P2>(4, 52, 3, 40,
 228                          1336.039890, -176.3086);
 229     test_distazi<P1, P2>(3, 52, 4, 40,
 230                          1336.039890, 176.3086);
 231     test_distazi<P1, P2>(make_deg(17, 19, 43.28),
 232                          make_deg(40, 30, 31.151),
 233                          18, 40,
 234                          80.323245,
 235                          make_deg(134, 27, 50.05));
 236
 237     // antipodal
 238     // ok? in those cases shorter path would pass through a pole
 239     // but 90 or -90 would be consistent with distance?
 240     test_distazi<P1, P2>(0, 0,  180, 0, 20003.9, 0.0);
 241     test_distazi<P1, P2>(0, 0, -180, 0, 20003.9, 0.0);
 242     test_distazi<P1, P2>(-90, 0, 90, 0, 20003.9, 0.0);
 243     test_distazi<P1, P2>(90, 0, -90, 0, 20003.9, 0.0);
 244
 245     // 0, 45, 90 ...
 246     for (int i = 0 ; i < 360 ; i += 45)
 247     {
 248         // 0 45 90 ...
 249         double l = normlized_deg(i);
 250         // -1 44 89 ...
 251         double l1 = normlized_deg(i - 1);
 252         // 1 46 91 ...
 253         double l2 = normlized_deg(i + 1);
 254
 255         // near equator
 256         test_distazi_symm<P1, P2>(l1, -1, l2, 1, 313.7956, 45.1964);
 257
 258         // near poles
 259         test_distazi_symmNS<P1, P2>(l, -89.5, l, 89.5, 19892.2, 0.0);
 260         test_distazi_symmNS<P1, P2>(l, -89.6, l, 89.6, 19914.6, 0.0);
 261         test_distazi_symmNS<P1, P2>(l, -89.7, l, 89.7, 19936.9, 0.0);
 262         test_distazi_symmNS<P1, P2>(l, -89.8, l, 89.8, 19959.2, 0.0);
 263         test_distazi_symmNS<P1, P2>(l, -89.9, l, 89.9, 19981.6, 0.0);
 264         test_distazi_symmNS<P1, P2>(l, -89.99, l, 89.99, 20001.7, 0.0);
 265         test_distazi_symmNS<P1, P2>(l, -89.999, l, 89.999, 20003.7, 0.0);
 266         // antipodal
 267         test_distazi_symmNS<P1, P2>(l, -90, l, 90, 20003.9, 0.0);
 268
 269         test_distazi_symm<P1, P2>(normlized_deg(l-10.0), -10.0, normlized_deg(l+135), 45, 14892.1, 34.1802);
 270         test_distazi_symm<P1, P2>(normlized_deg(l-30.0), -30.0, normlized_deg(l+135), 45, 17890.7, 33.7002);
 271         test_distazi_symm<P1, P2>(normlized_deg(l-40.0), -40.0, normlized_deg(l+135), 45, 19319.7, 33.4801);
 272         test_distazi_symm<P1, P2>(normlized_deg(l-41.0), -41.0, normlized_deg(l+135), 45, 19459.1, 33.2408);
 273         test_distazi_symm<P1, P2>(normlized_deg(l-42.0), -42.0, normlized_deg(l+135), 45, 19597.8, 32.7844);
 274         test_distazi_symm<P1, P2>(normlized_deg(l-43.0), -43.0, normlized_deg(l+135), 45, 19735.8, 31.7784);
 275         test_distazi_symm<P1, P2>(normlized_deg(l-44.0), -44.0, normlized_deg(l+135), 45, 19873.0, 28.5588);
 276         test_distazi_symm<P1, P2>(normlized_deg(l-44.1), -44.1, normlized_deg(l+135), 45, 19886.7, 27.8304);
 277         test_distazi_symm<P1, P2>(normlized_deg(l-44.2), -44.2, normlized_deg(l+135), 45, 19900.4, 26.9173);
 278         test_distazi_symm<P1, P2>(normlized_deg(l-44.3), -44.3, normlized_deg(l+135), 45, 19914.1, 25.7401);
 279         test_distazi_symm<P1, P2>(normlized_deg(l-44.4), -44.4, normlized_deg(l+135), 45, 19927.7, 24.1668);
 280         test_distazi_symm<P1, P2>(normlized_deg(l-44.5), -44.5, normlized_deg(l+135), 45, 19941.4, 21.9599);
 281         test_distazi_symm<P1, P2>(normlized_deg(l-44.6), -44.6, normlized_deg(l+135), 45, 19955.0, 18.6438);
 282         test_distazi_symm<P1, P2>(normlized_deg(l-44.7), -44.7, normlized_deg(l+135), 45, 19968.6, 13.1096);
 283         test_distazi_symm<P1, P2>(normlized_deg(l-44.8), -44.8, normlized_deg(l+135), 45, 19982.3, 2.0300);
 284         // nearly antipodal
 285         test_distazi_symm<P1, P2>(normlized_deg(l-44.9), -44.9, normlized_deg(l+135), 45, 19995.9, 0.0);
 286         test_distazi_symm<P1, P2>(normlized_deg(l-44.95), -44.95, normlized_deg(l+135), 45, 20002.7, 0.0);
 287         test_distazi_symm<P1, P2>(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0);
 288         test_distazi_symm<P1, P2>(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0);
 289         // antipodal
 290         test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20003.92, 0.0, true);
 291     }
 292
 293     /* SQL Server gives:
 294         1116.82586908528, 0, 1336.02721932545
 295
 296        with:
 297 SELECT 0.001 * geography::STGeomFromText('POINT(0 90)', 4326).STDistance(geography::STGeomFromText('POINT(1 80)', 4326))
 298 union SELECT 0.001 * geography::STGeomFromText('POINT(4 52)', 4326).STDistance(geography::STGeomFromText('POINT(4 52)', 4326))
 299 union SELECT 0.001 * geography::STGeomFromText('POINT(4 52)', 4326).STDistance(geography::STGeomFromText('POINT(3 40)', 4326))
 300      */
 301
 302
 303     test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0);
 304     test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0);
 305     test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0);
 306     test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0);
 307
 308     test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left
 309     test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right
 310
 311     test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right
 312     test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left
 313     test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right
 314     test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left
 315 }
 316
 317 template <typename P>
 318 void test_all()
 319 {
 320     test_all<P, P>();
 321 }
 322
 323 int test_main(int, char* [])
 324 {
 325     //test_all<float[2]>();
 326     //test_all<double[2]>();
 327     //test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();
 328     //test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
 329     test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
 330
 331 #if defined(HAVE_TTMATH)
 332     test_all<bg::model::point<ttmath::Big<1,4>, 2, bg::cs::geographic<bg::degree> > >();
 333     test_all<bg::model::point<ttmath_big, 2, bg::cs::geographic<bg::degree> > >();
 334 #endif
 335
 336     return 0;
 337 }