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-2019, 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 } // namespace detail
 318 #endif // DOXYGEN_NO_DETAIL
 319
 320
 321 /*!
 322 \brief Short utility to normalize the coordinates on a spheroid
 323 \tparam Units The units of the coordindate system in the spheroid
 324 \tparam CoordinateType The type of the coordinates
 325 \param longitude Longitude
 326 \param latitude Latitude
 327 \ingroup utility
 328 */
 329 template <typename Units, typename CoordinateType>
 330 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
 331                                              CoordinateType& latitude)
 332 {
 333     detail::normalize_spheroidal_coordinates
 334         <
 335             Units, CoordinateType
 336         >::apply(longitude, latitude);
 337 }
 338
 339 template <typename Units, bool IsEquatorial, typename CoordinateType>
 340 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
 341                                              CoordinateType& latitude)
 342 {
 343     detail::normalize_spheroidal_coordinates
 344         <
 345             Units, CoordinateType, IsEquatorial
 346         >::apply(longitude, latitude);
 347 }
 348
 349 /*!
 350 \brief Short utility to normalize the longitude on a spheroid.
 351        Note that in general both coordinates should be normalized at once.
 352        This utility is suitable e.g. for normalization of the difference of longitudes.
 353 \tparam Units The units of the coordindate system in the spheroid
 354 \tparam CoordinateType The type of the coordinates
 355 \param longitude Longitude
 356 \ingroup utility
 357 */
 358 template <typename Units, typename CoordinateType>
 359 inline void normalize_longitude(CoordinateType& longitude)
 360 {
 361     detail::normalize_spheroidal_coordinates
 362         <
 363             Units, CoordinateType
 364         >::apply(longitude);
 365 }
 366
 367 /*!
 368 \brief Short utility to normalize the azimuth on a spheroid
 369        in the range (-180, 180].
 370 \tparam Units The units of the coordindate system in the spheroid
 371 \tparam CoordinateType The type of the coordinates
 372 \param angle Angle
 373 \ingroup utility
 374 */
 375 template <typename Units, typename CoordinateType>
 376 inline void normalize_azimuth(CoordinateType& angle)
 377 {
 378     normalize_longitude<Units, CoordinateType>(angle);
 379 }
 380
 381 /*!
 382 \brief Normalize the given values.
 383 \tparam ValueType The type of the values
 384 \param x Value x
 385 \param y Value y
 386 TODO: adl1995 - Merge this function with
 387 formulas/vertex_longitude.hpp
 388 */
 389 template<typename ValueType>
 390 inline void normalize_unit_vector(ValueType& x, ValueType& y)
 391 {
 392     ValueType h = boost::math::hypot(x, y);
 393
 394     BOOST_GEOMETRY_ASSERT(h > 0);
 395
 396     x /= h; y /= h;
 397 }
 398
 399 /*!
 400 \brief Short utility to calculate difference between two longitudes
 401        normalized in range (-180, 180].
 402 \tparam Units The units of the coordindate system in the spheroid
 403 \tparam CoordinateType The type of the coordinates
 404 \param longitude1 Longitude 1
 405 \param longitude2 Longitude 2
 406 \ingroup utility
 407 */
 408 template <typename Units, typename CoordinateType>
 409 inline CoordinateType longitude_distance_signed(CoordinateType const& longitude1,
 410                                                 CoordinateType const& longitude2)
 411 {
 412     CoordinateType diff = longitude2 - longitude1;
 413     math::normalize_longitude<Units, CoordinateType>(diff);
 414     return diff;
 415 }
 416
 417
 418 /*!
 419 \brief Short utility to calculate difference between two longitudes
 420        normalized in range [0, 360).
 421 \tparam Units The units of the coordindate system in the spheroid
 422 \tparam CoordinateType The type of the coordinates
 423 \param longitude1 Longitude 1
 424 \param longitude2 Longitude 2
 425 \ingroup utility
 426 */
 427 template <typename Units, typename CoordinateType>
 428 inline CoordinateType longitude_distance_unsigned(CoordinateType const& longitude1,
 429                                                   CoordinateType const& longitude2)
 430 {
 431     typedef math::detail::constants_on_spheroid
 432         <
 433             CoordinateType, Units
 434         > constants;
 435
 436     CoordinateType const c0 = 0;
 437     CoordinateType diff = longitude_distance_signed<Units>(longitude1, longitude2);
 438     if (diff < c0) // (-180, 180] -> [0, 360)
 439     {
 440         diff += constants::period();
 441     }
 442     return diff;
 443 }
 444
 445 /*!
 446 \brief The abs difference between longitudes in range [0, 180].
 447 \tparam Units The units of the coordindate system in the spheroid
 448 \tparam CoordinateType The type of the coordinates
 449 \param longitude1 Longitude 1
 450 \param longitude2 Longitude 2
 451 \ingroup utility
 452 */
 453 template <typename Units, typename CoordinateType>
 454 inline CoordinateType longitude_difference(CoordinateType const& longitude1,
 455                                            CoordinateType const& longitude2)
 456 {
 457     return math::abs(math::longitude_distance_signed<Units>(longitude1, longitude2));
 458 }
 459
 460 template <typename Units, typename CoordinateType>
 461 inline CoordinateType longitude_interval_distance_signed(CoordinateType const& longitude_a1,
 462                                                          CoordinateType const& longitude_a2,
 463                                                          CoordinateType const& longitude_b)
 464 {
 465     CoordinateType const c0 = 0;
 466     CoordinateType dist_a12 = longitude_distance_signed<Units>(longitude_a1, longitude_a2);
 467     CoordinateType dist_a1b = longitude_distance_signed<Units>(longitude_a1, longitude_b);
 468     if (dist_a12 < c0)
 469     {
 470         dist_a12 = -dist_a12;
 471         dist_a1b = -dist_a1b;
 472     }
 473
 474     return dist_a1b < c0 ? dist_a1b
 475          : dist_a1b > dist_a12 ? dist_a1b - dist_a12
 476          : c0;
 477 }
 478
 479
 480 } // namespace math
 481
 482
 483 }} // namespace boost::geometry
 484
 485 #endif // BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP