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

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2 // Unit Test
   3
   4 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
   5 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
   6 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
   7
   8 // This file was modified by Oracle on 2014, 2015, 2016.
   9 // Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12
  13 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  14 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  15
  16 // Use, modification and distribution is subject to the Boost Software License,
  17 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  18 // http://www.boost.org/LICENSE_1_0.txt)
  19
  20
  21 #include <geometry_test_common.hpp>
  22
  23 #include <boost/concept_check.hpp>
  24
  25 #include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
  26 #include <boost/geometry/strategies/geographic/side_vincenty.hpp>
  27 #include <boost/geometry/formulas/vincenty_inverse.hpp>
  28 #include <boost/geometry/formulas/vincenty_direct.hpp>
  29
  30 #include <boost/geometry/core/srs.hpp>
  31 #include <boost/geometry/strategies/strategies.hpp>
  32 #include <boost/geometry/algorithms/assign.hpp>
  33 #include <boost/geometry/geometries/point.hpp>
  34 #include <test_common/test_point.hpp>
  35
  36 #ifdef HAVE_TTMATH
  37 #  include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
  38 #endif
  39
  40 template <typename T>
  41 void normalize_deg(T & deg)
  42 {
  43     while ( deg > T(180) )
  44         deg -= T(360);
  45     while ( deg <= T(-180) )
  46         deg += T(360);
  47 }
  48
  49 template <typename T>
  50 T difference_deg(T const& a1, T const& a2)
  51 {
  52     T d = a1 - a2;
  53     normalize_deg(d);
  54     return d;
  55 }
  56
  57 template <typename T>
  58 void check_deg(std::string const& name, T const& a1, T const& a2, T const& percent, T const& error)
  59 {
  60     T diff = bg::math::abs(difference_deg(a1, a2));
  61
  62     if ( bg::math::equals(a1, T(0)) || bg::math::equals(a2, T(0)) )
  63     {
  64         if ( diff > error )
  65         {
  66             BOOST_ERROR(name << " - the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << error << "}");
  67         }
  68     }
  69     else
  70     {
  71         T greater = (std::max)(bg::math::abs(a1), bg::math::abs(a2));
  72
  73         if ( diff > greater * percent / T(100) )
  74         {
  75             BOOST_ERROR(name << " the difference {" << diff << "} between {" << a1 << "} and {" << a2 << "} exceeds {" << percent << "}%");
  76         }
  77     }
  78 }
  79
  80 double azimuth(double deg, double min, double sec)
  81 {
  82     min = fabs(min);
  83     sec = fabs(sec);
  84
  85     if ( deg < 0 )
  86     {
  87         min = -min;
  88         sec = -sec;
  89     }
  90
  91     return deg + min/60.0 + sec/3600.0;
  92 }
  93
  94 double azimuth(double deg, double min)
  95 {
  96     return azimuth(deg, min, 0.0);
  97 }
  98
  99 template <typename P>
 100 bool non_precise_ct()
 101 {
 102     typedef typename bg::coordinate_type<P>::type ct;
 103     return boost::is_integral<ct>::value || boost::is_float<ct>::value;
 104 }
 105
 106 template <typename P1, typename P2, typename Spheroid>
 107 void test_vincenty(double lon1, double lat1, double lon2, double lat2,
 108                    double expected_distance,
 109                    double expected_azimuth_12,
 110                    double /*expected_azimuth_21*/,
 111                    Spheroid const& spheroid)
 112 {
 113     typedef typename bg::promote_floating_point
 114         <
 115             typename bg::select_calculation_type<P1, P2, void>::type
 116         >::type calc_t;
 117
 118     calc_t tolerance = non_precise_ct<P1>() || non_precise_ct<P2>() ?
 119                        5.0 : 0.001;
 120     calc_t error = non_precise_ct<P1>() || non_precise_ct<P2>() ?
 121                    1e-5 : 1e-12;
 122
 123     // formula
 124     {
 125         double const d2r = bg::math::d2r<double>();
 126         double const r2d = bg::math::r2d<double>();
 127
 128         typedef bg::formula::vincenty_inverse<calc_t, true, true> inverse_formula;
 129         typename inverse_formula::result_type
 130             result_i = inverse_formula::apply(lon1 * d2r,
 131                                              lat1 * d2r,
 132                                              lon2 * d2r,
 133                                              lat2 * d2r,
 134                                              spheroid);
 135         calc_t dist = result_i.distance;
 136         calc_t az12 = result_i.azimuth;
 137         //calc_t az21 = vi.azimuth21();
 138
 139         calc_t az12_deg = az12 * r2d;
 140         //calc_t az21_deg = az21 * r2d;
 141
 142         BOOST_CHECK_CLOSE(dist, calc_t(expected_distance), tolerance);
 143         check_deg("az12_deg", az12_deg, calc_t(expected_azimuth_12), tolerance, error);
 144         //check_deg("az21_deg", az21_deg, calc_t(expected_azimuth_21), tolerance, error);
 145
 146         typedef bg::formula::vincenty_direct<calc_t> direct_formula;
 147         typename direct_formula::result_type
 148             result_d = direct_formula::apply(lon1 * d2r,
 149                                              lat1 * d2r,
 150                                              dist,
 151                                              az12,
 152                                              spheroid);
 153         calc_t direct_lon2 = result_d.lon2;
 154         calc_t direct_lat2 = result_d.lat2;
 155         //calc_t direct_az21 = vd.azimuth21();
 156
 157         calc_t direct_lon2_deg = direct_lon2 * r2d;
 158         calc_t direct_lat2_deg = direct_lat2 * r2d;
 159         //calc_t direct_az21_deg = direct_az21 * r2d;
 160
 161         check_deg("direct_lon2_deg", direct_lon2_deg, calc_t(lon2), tolerance, error);
 162         check_deg("direct_lat2_deg", direct_lat2_deg, calc_t(lat2), tolerance, error);
 163         //check_deg("direct_az21_deg", direct_az21_deg, az21_deg, tolerance, error);
 164     }
 165
 166     // distance strategy
 167     {
 168         typedef bg::strategy::distance::vincenty<Spheroid> vincenty_type;
 169
 170         BOOST_CONCEPT_ASSERT(
 171             (
 172                 bg::concepts::PointDistanceStrategy<vincenty_type, P1, P2>)
 173             );
 174
 175         vincenty_type vincenty(spheroid);
 176         typedef typename bg::strategy::distance::services::return_type<vincenty_type, P1, P2>::type return_type;
 177
 178         P1 p1;
 179         P2 p2;
 180
 181         bg::assign_values(p1, lon1, lat1);
 182         bg::assign_values(p2, lon2, lat2);
 183
 184         BOOST_CHECK_CLOSE(vincenty.apply(p1, p2), return_type(expected_distance), tolerance);
 185         BOOST_CHECK_CLOSE(bg::distance(p1, p2, vincenty), return_type(expected_distance), tolerance);
 186     }
 187 }
 188
 189 template <typename P1, typename P2>
 190 void test_vincenty(double lon1, double lat1, double lon2, double lat2,
 191                    double expected_distance,
 192                    double expected_azimuth_12,
 193                    double expected_azimuth_21)
 194 {
 195     test_vincenty<P1, P2>(lon1, lat1, lon2, lat2,
 196                           expected_distance, expected_azimuth_12, expected_azimuth_21,
 197                           bg::srs::spheroid<double>());
 198 }
 199
 200 template <typename PS, typename P>
 201 void test_side(double lon1, double lat1,
 202                double lon2, double lat2,
 203                double lon, double lat,
 204                int expected_side)
 205 {
 206     // Set radius type, but for integer coordinates we want to have floating point radius type
 207     typedef typename bg::promote_floating_point
 208         <
 209             typename bg::coordinate_type<PS>::type
 210         >::type rtype;
 211
 212     typedef bg::srs::spheroid<rtype> stype;
 213
 214     typedef bg::strategy::side::vincenty<stype> strategy_type;
 215
 216     strategy_type strategy;
 217
 218     PS p1, p2;
 219     P p;
 220
 221     bg::assign_values(p1, lon1, lat1);
 222     bg::assign_values(p2, lon2, lat2);
 223     bg::assign_values(p, lon, lat);
 224
 225     int side = strategy.apply(p1, p2, p);
 226
 227     BOOST_CHECK_EQUAL(side, expected_side);
 228 }
 229
 230 template <typename P1, typename P2>
 231 void test_all()
 232 {
 233     // See:
 234     //  - http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp
 235     //  - http://www.ga.gov.au/geodesy/datums/vincenty_direct.jsp
 236     // Values in the comments below was calculated using the above pages
 237     // in some cases distances may be different, previously used values was left
 238
 239     // use km
 240     double gda_a = 6378.1370;
 241     double gda_f = 1.0 / 298.25722210;
 242     double gda_b = gda_a * ( 1.0 - gda_f );
 243     bg::srs::spheroid<double> gda_spheroid(gda_a, gda_b);
 244
 245     // Test fractional coordinates only for non-integral types
 246     if ( BOOST_GEOMETRY_CONDITION(
 247             ! boost::is_integral<typename bg::coordinate_type<P1>::type>::value
 248          && ! boost::is_integral<typename bg::coordinate_type<P2>::type>::value  ) )
 249     {
 250         // Flinders Peak -> Buninyong
 251         test_vincenty<P1, P2>(azimuth(144,25,29.52440), azimuth(-37,57,3.72030),
 252                               azimuth(143,55,35.38390), azimuth(-37,39,10.15610),
 253                               54.972271, azimuth(306,52,5.37), azimuth(127,10,25.07),
 254                               gda_spheroid);
 255     }
 256
 257     // Lodz -> Trondheim
 258     test_vincenty<P1, P2>(azimuth(19,28), azimuth(51,47),
 259                           azimuth(10,21), azimuth(63,23),
 260                           1399.032724, azimuth(340,54,25.14), azimuth(153,10,0.19),
 261                           gda_spheroid);
 262     // London -> New York
 263     test_vincenty<P1, P2>(azimuth(0,7,39), azimuth(51,30,26),
 264                           azimuth(-74,0,21), azimuth(40,42,46),
 265                           5602.044851, azimuth(288,31,36.82), azimuth(51,10,33.43),
 266                           gda_spheroid);
 267
 268     // Shanghai -> San Francisco
 269     test_vincenty<P1, P2>(azimuth(121,30), azimuth(31,12),
 270                           azimuth(-122,25), azimuth(37,47),
 271                           9899.698550, azimuth(45,12,44.76), azimuth(309,50,20.88),
 272                           gda_spheroid);
 273
 274     test_vincenty<P1, P2>(0, 0, 0, 50, 5540.847042, 0, 180, gda_spheroid); // N
 275     test_vincenty<P1, P2>(0, 0, 0, -50, 5540.847042, 180, 0, gda_spheroid); // S
 276     test_vincenty<P1, P2>(0, 0, 50, 0, 5565.974540, 90, -90, gda_spheroid); // E
 277     test_vincenty<P1, P2>(0, 0, -50, 0, 5565.974540, -90, 90, gda_spheroid); // W
 278
 279     test_vincenty<P1, P2>(0, 0, 50, 50, 7284.879297, azimuth(32,51,55.87), azimuth(237,24,50.12), gda_spheroid); // NE
 280
 281     // The original distance values, azimuths calculated using the web form mentioned above
 282     // Using default spheroid units (meters)
 283     test_vincenty<P1, P2>(0, 89, 1, 80, 1005153.5769, azimuth(178,53,23.85), azimuth(359,53,18.35)); // sub-polar
 284     test_vincenty<P1, P2>(4, 52, 4, 52, 0.0, 0, 0); // no point difference
 285     test_vincenty<P1, P2>(4, 52, 3, 40, 1336039.890, azimuth(183,41,29.08), azimuth(2,58,5.13)); // normal case
 286
 287     test_side<P1, P2>(0, 0, 0, 1, 0, 2, 0);
 288     test_side<P1, P2>(0, 0, 0, 1, 0, -2, 0);
 289     test_side<P1, P2>(10, 0, 10, 1, 10, 2, 0);
 290     test_side<P1, P2>(10, 0, 10, -1, 10, 2, 0);
 291
 292     test_side<P1, P2>(10, 0, 10, 1, 0, 2, 1); // left
 293     test_side<P1, P2>(10, 0, 10, -1, 0, 2, -1); // right
 294
 295     test_side<P1, P2>(-10, -10, 10, 10, 10, 0, -1); // right
 296     test_side<P1, P2>(-10, -10, 10, 10, -10, 0, 1); // left
 297     test_side<P1, P2>(170, -10, -170, 10, -170, 0, -1); // right
 298     test_side<P1, P2>(170, -10, -170, 10, 170, 0, 1); // left
 299 }
 300
 301 template <typename P>
 302 void test_all()
 303 {
 304     test_all<P, P>();
 305 }
 306
 307 int test_main(int, char* [])
 308 {
 309     //test_all<float[2]>();
 310     //test_all<double[2]>();
 311     test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
 312     test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
 313     test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();
 314
 315 #if defined(HAVE_TTMATH)
 316     test_all<bg::model::point<ttmath::Big<1,4>, 2, bg::cs::geographic<bg::degree> > >();
 317     test_all<bg::model::point<ttmath_big, 2, bg::cs::geographic<bg::degree> > >();
 318 #endif
 319
 320
 321     return 0;
 322 }