ceph/src/boost/boost/geometry/util/normalize_spheroidal_coordinates.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
   4
   5 // Copyright (c) 2015-2020, Oracle and/or its affiliates.
   6
   7 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
   8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   9 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program
  10
  11 // Licensed under the Boost Software License version 1.0.
  12 // http://www.boost.org/users/license.html
  13
  14 #ifndef BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP
  15 #define BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP
  16
  17 #include <boost/geometry/core/assert.hpp>
  18 #include <boost/geometry/core/cs.hpp>
  19 #include <boost/geometry/util/math.hpp>
  20
  21
  22 namespace boost { namespace geometry
  23 {
  24
  25 namespace math
  26 {
  27
  28 #ifndef DOXYGEN_NO_DETAIL
  29 namespace detail
  30 {
  31
  32 // CoordinateType, radian, true
  33 template <typename CoordinateType, typename Units, bool IsEquatorial = true>
  34 struct constants_on_spheroid
  35 {
  36     static inline CoordinateType period()
  37     {
  38         return math::two_pi<CoordinateType>();
  39     }
  40
  41     static inline CoordinateType half_period()
  42     {
  43         return math::pi<CoordinateType>();
  44     }
  45
  46     static inline CoordinateType quarter_period()
  47     {
  48         static CoordinateType const
  49             pi_half = math::pi<CoordinateType>() / CoordinateType(2);
  50         return pi_half;
  51     }
  52
  53     static inline CoordinateType min_longitude()
  54     {
  55         static CoordinateType const minus_pi = -math::pi<CoordinateType>();
  56         return minus_pi;
  57     }
  58
  59     static inline CoordinateType max_longitude()
  60     {
  61         return math::pi<CoordinateType>();
  62     }
  63
  64     static inline CoordinateType min_latitude()
  65     {
  66         static CoordinateType const minus_half_pi
  67             = -math::half_pi<CoordinateType>();
  68         return minus_half_pi;
  69     }
  70
  71     static inline CoordinateType max_latitude()
  72     {
  73         return math::half_pi<CoordinateType>();
  74     }
  75 };
  76
  77 template <typename CoordinateType>
  78 struct constants_on_spheroid<CoordinateType, radian, false>
  79     : constants_on_spheroid<CoordinateType, radian, true>
  80 {
  81     static inline CoordinateType min_latitude()
  82     {
  83         return CoordinateType(0);
  84     }
  85
  86     static inline CoordinateType max_latitude()
  87     {
  88         return math::pi<CoordinateType>();
  89     }
  90 };
  91
  92 template <typename CoordinateType>
  93 struct constants_on_spheroid<CoordinateType, degree, true>
  94 {
  95     static inline CoordinateType period()
  96     {
  97         return CoordinateType(360.0);
  98     }
  99
 100     static inline CoordinateType half_period()
 101     {
 102         return CoordinateType(180.0);
 103     }
 104
 105     static inline CoordinateType quarter_period()
 106     {
 107         return CoordinateType(90.0);
 108     }
 109
 110     static inline CoordinateType min_longitude()
 111     {
 112         return CoordinateType(-180.0);
 113     }
 114
 115     static inline CoordinateType max_longitude()
 116     {
 117         return CoordinateType(180.0);
 118     }
 119
 120     static inline CoordinateType min_latitude()
 121     {
 122         return CoordinateType(-90.0);
 123     }
 124
 125     static inline CoordinateType max_latitude()
 126     {
 127         return CoordinateType(90.0);
 128     }
 129 };
 130
 131 template <typename CoordinateType>
 132 struct constants_on_spheroid<CoordinateType, degree, false>
 133     : constants_on_spheroid<CoordinateType, degree, true>
 134 {
 135     static inline CoordinateType min_latitude()
 136     {
 137         return CoordinateType(0);
 138     }
 139
 140     static inline CoordinateType max_latitude()
 141     {
 142         return CoordinateType(180.0);
 143     }
 144 };
 145
 146
 147 } // namespace detail
 148 #endif // DOXYGEN_NO_DETAIL
 149
 150
 151 template <typename Units, typename CoordinateType>
 152 inline CoordinateType latitude_convert_ep(CoordinateType const& lat)
 153 {
 154     typedef math::detail::constants_on_spheroid
 155             <
 156                 CoordinateType,
 157                 Units
 158             > constants;
 159
 160     return constants::quarter_period() - lat;
 161 }
 162
 163
 164 template <typename Units, bool IsEquatorial, typename T>
 165 static bool is_latitude_pole(T const& lat)
 166 {
 167     typedef math::detail::constants_on_spheroid
 168         <
 169             T,
 170             Units
 171         > constants;
 172
 173     return math::equals(math::abs(IsEquatorial
 174                                     ? lat
 175                                     : math::latitude_convert_ep<Units>(lat)),
 176                         constants::quarter_period());
 177
 178 }
 179
 180
 181 template <typename Units, typename T>
 182 static bool is_longitude_antimeridian(T const& lon)
 183 {
 184     typedef math::detail::constants_on_spheroid
 185         <
 186             T,
 187             Units
 188         > constants;
 189
 190     return math::equals(math::abs(lon), constants::half_period());
 191
 192 }
 193
 194
 195 #ifndef DOXYGEN_NO_DETAIL
 196 namespace detail
 197 {
 198
 199
 200 template <typename Units, bool IsEquatorial>
 201 struct latitude_convert_if_polar
 202 {
 203     template <typename T>
 204     static inline void apply(T & /*lat*/) {}
 205 };
 206
 207 template <typename Units>
 208 struct latitude_convert_if_polar<Units, false>
 209 {
 210     template <typename T>
 211     static inline void apply(T & lat)
 212     {
 213         lat = latitude_convert_ep<Units>(lat);
 214     }
 215 };
 216
 217
 218 template <typename Units, typename CoordinateType, bool IsEquatorial = true>
 219 class normalize_spheroidal_coordinates
 220 {
 221     typedef constants_on_spheroid<CoordinateType, Units> constants;
 222
 223 protected:
 224     static inline CoordinateType normalize_up(CoordinateType const& value)
 225     {
 226         return
 227             math::mod(value + constants::half_period(), constants::period())
 228             - constants::half_period();
 229     }
 230
 231     static inline CoordinateType normalize_down(CoordinateType const& value)
 232     {
 233         return
 234             math::mod(value - constants::half_period(), constants::period())
 235             + constants::half_period();
 236     }
 237
 238 public:
 239     static inline void apply(CoordinateType& longitude)
 240     {
 241         // normalize longitude
 242         if (math::equals(math::abs(longitude), constants::half_period()))
 243         {
 244             longitude = constants::half_period();
 245         }
 246         else if (longitude > constants::half_period())
 247         {
 248             longitude = normalize_up(longitude);
 249             if (math::equals(longitude, -constants::half_period()))
 250             {
 251                 longitude = constants::half_period();
 252             }
 253         }
 254         else if (longitude < -constants::half_period())
 255         {
 256             longitude = normalize_down(longitude);
 257         }
 258     }
 259
 260     static inline void apply(CoordinateType& longitude,
 261                              CoordinateType& latitude,
 262                              bool normalize_poles = true)
 263     {
 264         latitude_convert_if_polar<Units, IsEquatorial>::apply(latitude);
 265
 266 #ifdef BOOST_GEOMETRY_NORMALIZE_LATITUDE
 267         // normalize latitude
 268         if (math::larger(latitude, constants::half_period()))
 269         {
 270             latitude = normalize_up(latitude);
 271         }
 272         else if (math::smaller(latitude, -constants::half_period()))
 273         {
 274             latitude = normalize_down(latitude);
 275         }
 276
 277         // fix latitude range
 278         if (latitude < constants::min_latitude())
 279         {
 280             latitude = -constants::half_period() - latitude;
 281             longitude -= constants::half_period();
 282         }
 283         else if (latitude > constants::max_latitude())
 284         {
 285             latitude = constants::half_period() - latitude;
 286             longitude -= constants::half_period();
 287         }
 288 #endif // BOOST_GEOMETRY_NORMALIZE_LATITUDE
 289
 290         // normalize longitude
 291         apply(longitude);
 292
 293         // finally normalize poles
 294         if (normalize_poles)
 295         {
 296             if (math::equals(math::abs(latitude), constants::max_latitude()))
 297             {
 298                 // for the north and south pole we set the longitude to 0
 299                 // (works for both radians and degrees)
 300                 longitude = CoordinateType(0);
 301             }
 302         }
 303
 304         latitude_convert_if_polar<Units, IsEquatorial>::apply(latitude);
 305
 306 #ifdef BOOST_GEOMETRY_NORMALIZE_LATITUDE
 307         BOOST_GEOMETRY_ASSERT(! math::larger(constants::min_latitude(), latitude));
 308         BOOST_GEOMETRY_ASSERT(! math::larger(latitude, constants::max_latitude()));
 309 #endif // BOOST_GEOMETRY_NORMALIZE_LATITUDE
 310
 311         BOOST_GEOMETRY_ASSERT(math::smaller(constants::min_longitude(), longitude));
 312         BOOST_GEOMETRY_ASSERT(! math::larger(longitude, constants::max_longitude()));
 313     }
 314 };
 315
 316
 317 template <typename Units, typename CoordinateType>
 318 inline void normalize_angle_loop(CoordinateType& angle)
 319 {
 320     typedef constants_on_spheroid<CoordinateType, Units> constants;
 321     CoordinateType const pi = constants::half_period();
 322     CoordinateType const two_pi = constants::period();
 323     while (angle > pi)
 324         angle -= two_pi;
 325     while (angle <= -pi)
 326         angle += two_pi;
 327 }
 328
 329 template <typename Units, typename CoordinateType>
 330 inline void normalize_angle_cond(CoordinateType& angle)
 331 {
 332     typedef constants_on_spheroid<CoordinateType, Units> constants;
 333     CoordinateType const pi = constants::half_period();
 334     CoordinateType const two_pi = constants::period();
 335     if (angle > pi)
 336         angle -= two_pi;
 337     else if (angle <= -pi)
 338         angle += two_pi;
 339 }
 340
 341
 342 } // namespace detail
 343 #endif // DOXYGEN_NO_DETAIL
 344
 345
 346 /*!
 347 \brief Short utility to normalize the coordinates on a spheroid
 348 \tparam Units The units of the coordindate system in the spheroid
 349 \tparam CoordinateType The type of the coordinates
 350 \param longitude Longitude
 351 \param latitude Latitude
 352 \ingroup utility
 353 */
 354 template <typename Units, typename CoordinateType>
 355 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
 356                                              CoordinateType& latitude)
 357 {
 358     detail::normalize_spheroidal_coordinates
 359         <
 360             Units, CoordinateType
 361         >::apply(longitude, latitude);
 362 }
 363
 364 template <typename Units, bool IsEquatorial, typename CoordinateType>
 365 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
 366                                              CoordinateType& latitude)
 367 {
 368     detail::normalize_spheroidal_coordinates
 369         <
 370             Units, CoordinateType, IsEquatorial
 371         >::apply(longitude, latitude);
 372 }
 373
 374 /*!
 375 \brief Short utility to normalize the longitude on a spheroid.
 376        Note that in general both coordinates should be normalized at once.
 377        This utility is suitable e.g. for normalization of the difference of longitudes.
 378 \tparam Units The units of the coordindate system in the spheroid
 379 \tparam CoordinateType The type of the coordinates
 380 \param longitude Longitude
 381 \ingroup utility
 382 */
 383 template <typename Units, typename CoordinateType>
 384 inline void normalize_longitude(CoordinateType& longitude)
 385 {
 386     detail::normalize_spheroidal_coordinates
 387         <
 388             Units, CoordinateType
 389         >::apply(longitude);
 390 }
 391
 392 /*!
 393 \brief Short utility to normalize the azimuth on a spheroid
 394        in the range (-180, 180].
 395 \tparam Units The units of the coordindate system in the spheroid
 396 \tparam CoordinateType The type of the coordinates
 397 \param angle Angle
 398 \ingroup utility
 399 */
 400 template <typename Units, typename CoordinateType>
 401 inline void normalize_azimuth(CoordinateType& angle)
 402 {
 403     normalize_longitude<Units, CoordinateType>(angle);
 404 }
 405
 406 /*!
 407 \brief Normalize the given values.
 408 \tparam ValueType The type of the values
 409 \param x Value x
 410 \param y Value y
 411 TODO: adl1995 - Merge this function with
 412 formulas/vertex_longitude.hpp
 413 */
 414 template<typename ValueType>
 415 inline void normalize_unit_vector(ValueType& x, ValueType& y)
 416 {
 417     ValueType h = boost::math::hypot(x, y);
 418
 419     BOOST_GEOMETRY_ASSERT(h > 0);
 420
 421     x /= h; y /= h;
 422 }
 423
 424 /*!
 425 \brief Short utility to calculate difference between two longitudes
 426        normalized in range (-180, 180].
 427 \tparam Units The units of the coordindate system in the spheroid
 428 \tparam CoordinateType The type of the coordinates
 429 \param longitude1 Longitude 1
 430 \param longitude2 Longitude 2
 431 \ingroup utility
 432 */
 433 template <typename Units, typename CoordinateType>
 434 inline CoordinateType longitude_distance_signed(CoordinateType const& longitude1,
 435                                                 CoordinateType const& longitude2)
 436 {
 437     CoordinateType diff = longitude2 - longitude1;
 438     math::normalize_longitude<Units, CoordinateType>(diff);
 439     return diff;
 440 }
 441
 442
 443 /*!
 444 \brief Short utility to calculate difference between two longitudes
 445        normalized in range [0, 360).
 446 \tparam Units The units of the coordindate system in the spheroid
 447 \tparam CoordinateType The type of the coordinates
 448 \param longitude1 Longitude 1
 449 \param longitude2 Longitude 2
 450 \ingroup utility
 451 */
 452 template <typename Units, typename CoordinateType>
 453 inline CoordinateType longitude_distance_unsigned(CoordinateType const& longitude1,
 454                                                   CoordinateType const& longitude2)
 455 {
 456     typedef math::detail::constants_on_spheroid
 457         <
 458             CoordinateType, Units
 459         > constants;
 460
 461     CoordinateType const c0 = 0;
 462     CoordinateType diff = longitude_distance_signed<Units>(longitude1, longitude2);
 463     if (diff < c0) // (-180, 180] -> [0, 360)
 464     {
 465         diff += constants::period();
 466     }
 467     return diff;
 468 }
 469
 470 /*!
 471 \brief The abs difference between longitudes in range [0, 180].
 472 \tparam Units The units of the coordindate system in the spheroid
 473 \tparam CoordinateType The type of the coordinates
 474 \param longitude1 Longitude 1
 475 \param longitude2 Longitude 2
 476 \ingroup utility
 477 */
 478 template <typename Units, typename CoordinateType>
 479 inline CoordinateType longitude_difference(CoordinateType const& longitude1,
 480                                            CoordinateType const& longitude2)
 481 {
 482     return math::abs(math::longitude_distance_signed<Units>(longitude1, longitude2));
 483 }
 484
 485 template <typename Units, typename CoordinateType>
 486 inline CoordinateType longitude_interval_distance_signed(CoordinateType const& longitude_a1,
 487                                                          CoordinateType const& longitude_a2,
 488                                                          CoordinateType const& longitude_b)
 489 {
 490     CoordinateType const c0 = 0;
 491     CoordinateType dist_a12 = longitude_distance_signed<Units>(longitude_a1, longitude_a2);
 492     CoordinateType dist_a1b = longitude_distance_signed<Units>(longitude_a1, longitude_b);
 493     if (dist_a12 < c0)
 494     {
 495         dist_a12 = -dist_a12;
 496         dist_a1b = -dist_a1b;
 497     }
 498
 499     return dist_a1b < c0 ? dist_a1b
 500          : dist_a1b > dist_a12 ? dist_a1b - dist_a12
 501          : c0;
 502 }
 503
 504
 505 } // namespace math
 506
 507
 508 }} // namespace boost::geometry
 509
 510 #endif // BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP