geometry/util/normalize_spheroidal_coordinates.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2015-2025, Oracle and/or its affiliates.
0006
0007 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0008 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010 // Contributed and/or modified by Adeel Ahmad, as part of Google Summer of Code 2018 program
0011
0012 // Licensed under the Boost Software License version 1.0.
0013 // http://www.boost.org/users/license.html
0014
0015 #ifndef BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP
0016 #define BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP
0017
0018 #include <boost/geometry/core/assert.hpp>
0019 #include <boost/geometry/core/cs.hpp>
0020 #include <boost/geometry/util/math.hpp>
0021
0022
0023 namespace boost { namespace geometry
0024 {
0025
0026 namespace math
0027 {
0028
0029 #ifndef DOXYGEN_NO_DETAIL
0030 namespace detail
0031 {
0032
0033 // CoordinateType, radian, true
0034 template <typename CoordinateType, typename Units, bool IsEquatorial = true>
0035 struct constants_on_spheroid
0036 {
0037     static inline CoordinateType period()
0038     {
0039         return math::two_pi<CoordinateType>();
0040     }
0041
0042     static inline CoordinateType half_period()
0043     {
0044         return math::pi<CoordinateType>();
0045     }
0046
0047     static inline CoordinateType quarter_period()
0048     {
0049         static CoordinateType const
0050             pi_half = math::pi<CoordinateType>() / CoordinateType(2);
0051         return pi_half;
0052     }
0053
0054     static inline CoordinateType min_longitude()
0055     {
0056         static CoordinateType const minus_pi = -math::pi<CoordinateType>();
0057         return minus_pi;
0058     }
0059
0060     static inline CoordinateType max_longitude()
0061     {
0062         return math::pi<CoordinateType>();
0063     }
0064
0065     static inline CoordinateType min_latitude()
0066     {
0067         static CoordinateType const minus_half_pi
0068             = -math::half_pi<CoordinateType>();
0069         return minus_half_pi;
0070     }
0071
0072     static inline CoordinateType max_latitude()
0073     {
0074         return math::half_pi<CoordinateType>();
0075     }
0076 };
0077
0078 template <typename CoordinateType>
0079 struct constants_on_spheroid<CoordinateType, radian, false>
0080     : constants_on_spheroid<CoordinateType, radian, true>
0081 {
0082     static inline CoordinateType min_latitude()
0083     {
0084         return CoordinateType(0);
0085     }
0086
0087     static inline CoordinateType max_latitude()
0088     {
0089         return math::pi<CoordinateType>();
0090     }
0091 };
0092
0093 template <typename CoordinateType>
0094 struct constants_on_spheroid<CoordinateType, degree, true>
0095 {
0096     static inline CoordinateType period()
0097     {
0098         return CoordinateType(360.0);
0099     }
0100
0101     static inline CoordinateType half_period()
0102     {
0103         return CoordinateType(180.0);
0104     }
0105
0106     static inline CoordinateType quarter_period()
0107     {
0108         return CoordinateType(90.0);
0109     }
0110
0111     static inline CoordinateType min_longitude()
0112     {
0113         return CoordinateType(-180.0);
0114     }
0115
0116     static inline CoordinateType max_longitude()
0117     {
0118         return CoordinateType(180.0);
0119     }
0120
0121     static inline CoordinateType min_latitude()
0122     {
0123         return CoordinateType(-90.0);
0124     }
0125
0126     static inline CoordinateType max_latitude()
0127     {
0128         return CoordinateType(90.0);
0129     }
0130 };
0131
0132 template <typename CoordinateType>
0133 struct constants_on_spheroid<CoordinateType, degree, false>
0134     : constants_on_spheroid<CoordinateType, degree, true>
0135 {
0136     static inline CoordinateType min_latitude()
0137     {
0138         return CoordinateType(0);
0139     }
0140
0141     static inline CoordinateType max_latitude()
0142     {
0143         return CoordinateType(180.0);
0144     }
0145 };
0146
0147
0148 } // namespace detail
0149 #endif // DOXYGEN_NO_DETAIL
0150
0151
0152 template <typename Units, typename CoordinateType>
0153 inline CoordinateType latitude_convert_ep(CoordinateType const& lat)
0154 {
0155     using constants = math::detail::constants_on_spheroid<CoordinateType, Units>;
0156
0157     return constants::quarter_period() - lat;
0158 }
0159
0160
0161 template <typename Units, bool IsEquatorial, typename T>
0162 static bool is_latitude_pole(T const& lat)
0163 {
0164     using constants = math::detail::constants_on_spheroid<T, Units>;
0165
0166     return math::equals(math::abs(IsEquatorial
0167                                     ? lat
0168                                     : math::latitude_convert_ep<Units>(lat)),
0169                         constants::quarter_period());
0170 }
0171
0172
0173 template <typename Units, typename T>
0174 static bool is_longitude_antimeridian(T const& lon)
0175 {
0176     using constants = math::detail::constants_on_spheroid<T, Units>;
0177
0178     return math::equals(math::abs(lon), constants::half_period());
0179 }
0180
0181
0182 #ifndef DOXYGEN_NO_DETAIL
0183 namespace detail
0184 {
0185
0186
0187 template <typename Units, bool IsEquatorial>
0188 struct latitude_convert_if_polar
0189 {
0190     template <typename T>
0191     static inline void apply(T & /*lat*/) {}
0192 };
0193
0194 template <typename Units>
0195 struct latitude_convert_if_polar<Units, false>
0196 {
0197     template <typename T>
0198     static inline void apply(T & lat)
0199     {
0200         lat = latitude_convert_ep<Units>(lat);
0201     }
0202 };
0203
0204
0205 template <typename Units, typename CoordinateType, bool IsEquatorial = true>
0206 class normalize_spheroidal_coordinates
0207 {
0208     using constants = constants_on_spheroid<CoordinateType, Units>;
0209
0210 protected:
0211     static inline CoordinateType normalize_up(CoordinateType const& value)
0212     {
0213         return
0214             math::mod(value + constants::half_period(), constants::period())
0215             - constants::half_period();
0216     }
0217
0218     static inline CoordinateType normalize_down(CoordinateType const& value)
0219     {
0220         return
0221             math::mod(value - constants::half_period(), constants::period())
0222             + constants::half_period();
0223     }
0224
0225 public:
0226     static inline void apply(CoordinateType& longitude, bool exact = true)
0227     {
0228         // normalize longitude
0229         CoordinateType const epsilon = std::numeric_limits<float>::epsilon();
0230         static constexpr bool is_integer = std::numeric_limits<CoordinateType>::is_integer;
0231
0232         if (exact || is_integer ? math::equals(math::abs(longitude), constants::half_period())
0233             : math::abs(math::abs(longitude) - constants::half_period()) <= epsilon)
0234         {
0235             longitude = constants::half_period();
0236         }
0237         else if (longitude > constants::half_period())
0238         {
0239             longitude = normalize_up(longitude);
0240             if (exact || is_integer ? math::equals(longitude, -constants::half_period())
0241                 : math::abs(longitude + constants::half_period()) <= epsilon)
0242             {
0243                 longitude = constants::half_period();
0244             }
0245         }
0246         else if (longitude < -constants::half_period())
0247         {
0248             longitude = normalize_down(longitude);
0249         }
0250     }
0251
0252     static inline void apply(CoordinateType& longitude,
0253                              CoordinateType& latitude,
0254                              bool normalize_poles = true,
0255                              bool exact = true)
0256     {
0257         latitude_convert_if_polar<Units, IsEquatorial>::apply(latitude);
0258
0259 #ifdef BOOST_GEOMETRY_NORMALIZE_LATITUDE
0260         // normalize latitude
0261         if (math::larger(latitude, constants::half_period()))
0262         {
0263             latitude = normalize_up(latitude);
0264         }
0265         else if (math::smaller(latitude, -constants::half_period()))
0266         {
0267             latitude = normalize_down(latitude);
0268         }
0269
0270         // fix latitude range
0271         if (latitude < constants::min_latitude())
0272         {
0273             latitude = -constants::half_period() - latitude;
0274             longitude -= constants::half_period();
0275         }
0276         else if (latitude > constants::max_latitude())
0277         {
0278             latitude = constants::half_period() - latitude;
0279             longitude -= constants::half_period();
0280         }
0281 #endif // BOOST_GEOMETRY_NORMALIZE_LATITUDE
0282
0283         // normalize longitude
0284         apply(longitude, exact);
0285
0286         // finally normalize poles
0287         if (normalize_poles)
0288         {
0289             if (math::equals(math::abs(latitude), constants::max_latitude()))
0290             {
0291                 // for the north and south pole we set the longitude to 0
0292                 // (works for both radians and degrees)
0293                 longitude = CoordinateType(0);
0294             }
0295         }
0296
0297         latitude_convert_if_polar<Units, IsEquatorial>::apply(latitude);
0298
0299 #ifdef BOOST_GEOMETRY_NORMALIZE_LATITUDE
0300         BOOST_GEOMETRY_ASSERT(! math::larger(constants::min_latitude(), latitude));
0301         BOOST_GEOMETRY_ASSERT(! math::larger(latitude, constants::max_latitude()));
0302 #endif // BOOST_GEOMETRY_NORMALIZE_LATITUDE
0303
0304         BOOST_GEOMETRY_ASSERT(! math::larger_or_equals(constants::min_longitude(), longitude));
0305         BOOST_GEOMETRY_ASSERT(! math::larger(longitude, constants::max_longitude()));
0306     }
0307 };
0308
0309
0310 template <typename Units, typename CoordinateType>
0311 inline void normalize_angle_loop(CoordinateType& angle)
0312 {
0313     using constants = constants_on_spheroid<CoordinateType, Units>;
0314     CoordinateType const pi = constants::half_period();
0315     CoordinateType const two_pi = constants::period();
0316     while (angle > pi)
0317         angle -= two_pi;
0318     while (angle <= -pi)
0319         angle += two_pi;
0320 }
0321
0322 template <typename Units, typename CoordinateType>
0323 inline void normalize_angle_cond(CoordinateType& angle)
0324 {
0325     using constants = constants_on_spheroid<CoordinateType, Units>;
0326     CoordinateType const pi = constants::half_period();
0327     CoordinateType const two_pi = constants::period();
0328     if (angle > pi)
0329         angle -= two_pi;
0330     else if (angle <= -pi)
0331         angle += two_pi;
0332 }
0333
0334
0335 } // namespace detail
0336 #endif // DOXYGEN_NO_DETAIL
0337
0338
0339 /*!
0340 \brief Short utility to normalize the coordinates on a spheroid
0341 \tparam Units The units of the coordindate system in the spheroid
0342 \tparam CoordinateType The type of the coordinates
0343 \param longitude Longitude
0344 \param latitude Latitude
0345 \ingroup utility
0346 */
0347 template <typename Units, typename CoordinateType>
0348 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
0349                                              CoordinateType& latitude,
0350                                              bool exact = true)
0351 {
0352     detail::normalize_spheroidal_coordinates
0353         <
0354             Units, CoordinateType
0355         >::apply(longitude, latitude, true, exact);
0356 }
0357
0358 template <typename Units, bool IsEquatorial, typename CoordinateType>
0359 inline void normalize_spheroidal_coordinates(CoordinateType& longitude,
0360                                              CoordinateType& latitude,
0361                                              bool exact = true)
0362 {
0363     detail::normalize_spheroidal_coordinates
0364         <
0365             Units, CoordinateType, IsEquatorial
0366         >::apply(longitude, latitude, true, exact);
0367 }
0368
0369 /*!
0370 \brief Short utility to normalize the longitude on a spheroid.
0371        Note that in general both coordinates should be normalized at once.
0372        This utility is suitable e.g. for normalization of the difference of longitudes.
0373 \tparam Units The units of the coordindate system in the spheroid
0374 \tparam CoordinateType The type of the coordinates
0375 \param longitude Longitude
0376 \ingroup utility
0377 */
0378 template <typename Units, typename CoordinateType>
0379 inline void normalize_longitude(CoordinateType& longitude, bool exact = true)
0380 {
0381     detail::normalize_spheroidal_coordinates
0382         <
0383             Units, CoordinateType
0384         >::apply(longitude, exact);
0385 }
0386
0387 /*!
0388 \brief Short utility to normalize the azimuth on a spheroid
0389        in the range (-180, 180].
0390 \tparam Units The units of the coordindate system in the spheroid
0391 \tparam CoordinateType The type of the coordinates
0392 \param angle Angle
0393 \ingroup utility
0394 */
0395 template <typename Units, typename CoordinateType>
0396 inline void normalize_azimuth(CoordinateType& angle)
0397 {
0398     math::normalize_longitude<Units, CoordinateType>(angle, true);
0399 }
0400
0401 /*!
0402 \brief Normalize the given values.
0403 \tparam ValueType The type of the values
0404 \param x Value x
0405 \param y Value y
0406 TODO: adl1995 - Merge this function with
0407 formulas/vertex_longitude.hpp
0408 */
0409 template<typename ValueType>
0410 inline void normalize_unit_vector(ValueType& x, ValueType& y)
0411 {
0412     ValueType h = boost::math::hypot(x, y);
0413
0414     BOOST_GEOMETRY_ASSERT(h > 0);
0415
0416     x /= h; y /= h;
0417 }
0418
0419 /*!
0420 \brief Short utility to calculate difference between two longitudes
0421        normalized in range (-180, 180].
0422 \tparam Units The units of the coordindate system in the spheroid
0423 \tparam CoordinateType The type of the coordinates
0424 \param longitude1 Longitude 1
0425 \param longitude2 Longitude 2
0426 \ingroup utility
0427 */
0428 template <typename Units, typename CoordinateType>
0429 inline CoordinateType longitude_distance_signed(CoordinateType const& longitude1,
0430                                                 CoordinateType const& longitude2)
0431 {
0432     CoordinateType diff = longitude2 - longitude1;
0433     math::normalize_longitude<Units, CoordinateType>(diff, true);
0434     return diff;
0435 }
0436
0437
0438 /*!
0439 \brief Short utility to calculate difference between two longitudes
0440        normalized in range [0, 360).
0441 \tparam Units The units of the coordindate system in the spheroid
0442 \tparam CoordinateType The type of the coordinates
0443 \param longitude1 Longitude 1
0444 \param longitude2 Longitude 2
0445 \ingroup utility
0446 */
0447 template <typename Units, typename CoordinateType>
0448 inline CoordinateType longitude_distance_unsigned(CoordinateType const& longitude1,
0449                                                   CoordinateType const& longitude2)
0450 {
0451     using constants = math::detail::constants_on_spheroid<CoordinateType, Units>;
0452
0453     CoordinateType const c0 = 0;
0454     CoordinateType diff = longitude_distance_signed<Units>(longitude1, longitude2);
0455     if (diff < c0) // (-180, 180] -> [0, 360)
0456     {
0457         diff += constants::period();
0458     }
0459     return diff;
0460 }
0461
0462 /*!
0463 \brief The abs difference between longitudes in range [0, 180].
0464 \tparam Units The units of the coordindate system in the spheroid
0465 \tparam CoordinateType The type of the coordinates
0466 \param longitude1 Longitude 1
0467 \param longitude2 Longitude 2
0468 \ingroup utility
0469 */
0470 template <typename Units, typename CoordinateType>
0471 inline CoordinateType longitude_difference(CoordinateType const& longitude1,
0472                                            CoordinateType const& longitude2)
0473 {
0474     return math::abs(math::longitude_distance_signed<Units>(longitude1, longitude2));
0475 }
0476
0477 template <typename Units, typename CoordinateType>
0478 inline CoordinateType longitude_interval_distance_signed(CoordinateType const& longitude_a1,
0479                                                          CoordinateType const& longitude_a2,
0480                                                          CoordinateType const& longitude_b)
0481 {
0482     CoordinateType const c0 = 0;
0483     CoordinateType dist_a12 = longitude_distance_signed<Units>(longitude_a1, longitude_a2);
0484     CoordinateType dist_a1b = longitude_distance_signed<Units>(longitude_a1, longitude_b);
0485     if (dist_a12 < c0)
0486     {
0487         dist_a12 = -dist_a12;
0488         dist_a1b = -dist_a1b;
0489     }
0490
0491     return dist_a1b < c0 ? dist_a1b
0492          : dist_a1b > dist_a12 ? dist_a1b - dist_a12
0493          : c0;
0494 }
0495
0496
0497 } // namespace math
0498
0499
0500 }} // namespace boost::geometry
0501
0502 #endif // BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP