ceph/src/boost/boost/geometry/strategies/spherical/envelope_segment.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2017-2018 Oracle and/or its affiliates.
   4 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
   5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   6
   7 // Use, modification and distribution is subject to the Boost Software License,
   8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   9 // http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_SEGMENT_HPP
  12 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_SEGMENT_HPP
  13
  14
  15 #include <cstddef>
  16 #include <utility>
  17
  18 #include <boost/core/ignore_unused.hpp>
  19 #include <boost/mpl/if.hpp>
  20 #include <boost/numeric/conversion/cast.hpp>
  21 #include <boost/type_traits/is_same.hpp>
  22
  23 #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp>
  24
  25 #include <boost/geometry/core/assert.hpp>
  26 #include <boost/geometry/core/coordinate_system.hpp>
  27 #include <boost/geometry/core/coordinate_type.hpp>
  28 #include <boost/geometry/core/cs.hpp>
  29 #include <boost/geometry/core/point_type.hpp>
  30 #include <boost/geometry/core/radian_access.hpp>
  31 #include <boost/geometry/core/tags.hpp>
  32
  33 #include <boost/geometry/formulas/meridian_segment.hpp>
  34 #include <boost/geometry/formulas/vertex_latitude.hpp>
  35
  36 #include <boost/geometry/geometries/helper_geometry.hpp>
  37
  38 #include <boost/geometry/strategies/cartesian/envelope_segment.hpp>
  39 #include <boost/geometry/strategies/envelope.hpp>
  40 #include <boost/geometry/strategies/normalize.hpp>
  41 #include <boost/geometry/strategies/spherical/azimuth.hpp>
  42 #include <boost/geometry/strategies/spherical/expand_box.hpp>
  43
  44 #include <boost/geometry/util/math.hpp>
  45
  46 namespace boost { namespace geometry { namespace strategy { namespace envelope
  47 {
  48
  49 #ifndef DOXYGEN_NO_DETAIL
  50 namespace detail
  51 {
  52
  53 template <typename CalculationType, typename CS_Tag>
  54 struct envelope_segment_call_vertex_latitude
  55 {
  56     template <typename T1, typename T2, typename Strategy>
  57     static inline CalculationType apply(T1 const& lat1,
  58                                         T2 const& alp1,
  59                                         Strategy const& )
  60     {
  61         return geometry::formula::vertex_latitude<CalculationType, CS_Tag>
  62             ::apply(lat1, alp1);
  63     }
  64 };
  65
  66 template <typename CalculationType>
  67 struct envelope_segment_call_vertex_latitude<CalculationType, geographic_tag>
  68 {
  69     template <typename T1, typename T2, typename Strategy>
  70     static inline CalculationType apply(T1 const& lat1,
  71                                         T2 const& alp1,
  72                                         Strategy const& strategy)
  73     {
  74         return geometry::formula::vertex_latitude<CalculationType, geographic_tag>
  75             ::apply(lat1, alp1, strategy.model());
  76     }
  77 };
  78
  79 template <typename Units, typename CS_Tag>
  80 struct envelope_segment_convert_polar
  81 {
  82     template <typename T>
  83     static inline void pre(T & , T & ) {}
  84
  85     template <typename T>
  86     static inline void post(T & , T & ) {}
  87 };
  88
  89 template <typename Units>
  90 struct envelope_segment_convert_polar<Units, spherical_polar_tag>
  91 {
  92     template <typename T>
  93     static inline void pre(T & lat1, T & lat2)
  94     {
  95         lat1 = math::latitude_convert_ep<Units>(lat1);
  96         lat2 = math::latitude_convert_ep<Units>(lat2);
  97     }
  98
  99     template <typename T>
 100     static inline void post(T & lat1, T & lat2)
 101     {
 102         lat1 = math::latitude_convert_ep<Units>(lat1);
 103         lat2 = math::latitude_convert_ep<Units>(lat2);
 104         std::swap(lat1, lat2);
 105     }
 106 };
 107
 108 template <typename CS_Tag>
 109 class envelope_segment_impl
 110 {
 111 private:
 112
 113     // degrees or radians
 114     template <typename CalculationType>
 115     static inline void swap(CalculationType& lon1,
 116                             CalculationType& lat1,
 117                             CalculationType& lon2,
 118                             CalculationType& lat2)
 119     {
 120         std::swap(lon1, lon2);
 121         std::swap(lat1, lat2);
 122     }
 123
 124     // radians
 125     template <typename CalculationType>
 126     static inline bool contains_pi_half(CalculationType const& a1,
 127                                         CalculationType const& a2)
 128     {
 129         // azimuths a1 and a2 are assumed to be in radians
 130         BOOST_GEOMETRY_ASSERT(! math::equals(a1, a2));
 131
 132         static CalculationType const pi_half = math::half_pi<CalculationType>();
 133
 134         return (a1 < a2)
 135                 ? (a1 < pi_half && pi_half < a2)
 136                 : (a1 > pi_half && pi_half > a2);
 137     }
 138
 139     // radians or degrees
 140     template <typename Units, typename CoordinateType>
 141     static inline bool crosses_antimeridian(CoordinateType const& lon1,
 142                                             CoordinateType const& lon2)
 143     {
 144         typedef math::detail::constants_on_spheroid
 145             <
 146                 CoordinateType, Units
 147             > constants;
 148
 149         return math::abs(lon1 - lon2) > constants::half_period(); // > pi
 150     }
 151
 152     // degrees or radians
 153     template <typename Units, typename CalculationType, typename Strategy>
 154     static inline void compute_box_corners(CalculationType& lon1,
 155                                            CalculationType& lat1,
 156                                            CalculationType& lon2,
 157                                            CalculationType& lat2,
 158                                            CalculationType a1,
 159                                            CalculationType a2,
 160                                            Strategy const& strategy)
 161     {
 162         // coordinates are assumed to be in radians
 163         BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
 164         boost::ignore_unused(lon1, lon2);
 165
 166         CalculationType lat1_rad = math::as_radian<Units>(lat1);
 167         CalculationType lat2_rad = math::as_radian<Units>(lat2);
 168
 169         if (math::equals(a1, a2))
 170         {
 171             // the segment must lie on the equator or is very short or is meridian
 172             return;
 173         }
 174
 175         if (lat1 > lat2)
 176         {
 177             std::swap(lat1, lat2);
 178             std::swap(lat1_rad, lat2_rad);
 179             std::swap(a1, a2);
 180         }
 181
 182         if (contains_pi_half(a1, a2))
 183         {
 184             CalculationType p_max = envelope_segment_call_vertex_latitude
 185                 <CalculationType, CS_Tag>::apply(lat1_rad, a1, strategy);
 186
 187             CalculationType const mid_lat = lat1 + lat2;
 188             if (mid_lat < 0)
 189             {
 190                 // update using min latitude
 191                 CalculationType const lat_min_rad = -p_max;
 192                 CalculationType const lat_min
 193                     = math::from_radian<Units>(lat_min_rad);
 194
 195                 if (lat1 > lat_min)
 196                 {
 197                     lat1 = lat_min;
 198                 }
 199             }
 200             else
 201             {
 202                 // update using max latitude
 203                 CalculationType const lat_max_rad = p_max;
 204                 CalculationType const lat_max
 205                     = math::from_radian<Units>(lat_max_rad);
 206
 207                 if (lat2 < lat_max)
 208                 {
 209                     lat2 = lat_max;
 210                 }
 211             }
 212         }
 213     }
 214
 215     template <typename Units, typename CalculationType>
 216     static inline void special_cases(CalculationType& lon1,
 217                                      CalculationType& lat1,
 218                                      CalculationType& lon2,
 219                                      CalculationType& lat2)
 220     {
 221         typedef math::detail::constants_on_spheroid
 222             <
 223                 CalculationType, Units
 224             > constants;
 225
 226         bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
 227         bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());
 228
 229         if (is_pole1 && is_pole2)
 230         {
 231             // both points are poles; nothing more to do:
 232             // longitudes are already normalized to 0
 233             // but just in case
 234             lon1 = 0;
 235             lon2 = 0;
 236         }
 237         else if (is_pole1 && !is_pole2)
 238         {
 239             // first point is a pole, second point is not:
 240             // make the longitude of the first point the same as that
 241             // of the second point
 242             lon1 = lon2;
 243         }
 244         else if (!is_pole1 && is_pole2)
 245         {
 246             // second point is a pole, first point is not:
 247             // make the longitude of the second point the same as that
 248             // of the first point
 249             lon2 = lon1;
 250         }
 251
 252         if (lon1 == lon2)
 253         {
 254             // segment lies on a meridian
 255             if (lat1 > lat2)
 256             {
 257                 std::swap(lat1, lat2);
 258             }
 259             return;
 260         }
 261
 262         BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);
 263
 264         if (lon1 > lon2)
 265         {
 266             swap(lon1, lat1, lon2, lat2);
 267         }
 268
 269         if (crosses_antimeridian<Units>(lon1, lon2))
 270         {
 271             lon1 += constants::period();
 272             swap(lon1, lat1, lon2, lat2);
 273         }
 274     }
 275
 276     template
 277     <
 278         typename Units,
 279         typename CalculationType,
 280         typename Box
 281     >
 282     static inline void create_box(CalculationType lon1,
 283                                   CalculationType lat1,
 284                                   CalculationType lon2,
 285                                   CalculationType lat2,
 286                                   Box& mbr)
 287     {
 288         typedef typename coordinate_type<Box>::type box_coordinate_type;
 289
 290         typedef typename helper_geometry
 291             <
 292                 Box, box_coordinate_type, Units
 293             >::type helper_box_type;
 294
 295         helper_box_type helper_mbr;
 296
 297         geometry::set
 298             <
 299                 min_corner, 0
 300             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lon1));
 301
 302         geometry::set
 303             <
 304                 min_corner, 1
 305             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lat1));
 306
 307         geometry::set
 308             <
 309                 max_corner, 0
 310             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lon2));
 311
 312         geometry::set
 313             <
 314                 max_corner, 1
 315             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lat2));
 316
 317         geometry::detail::envelope::transform_units(helper_mbr, mbr);
 318     }
 319
 320
 321     template <typename Units, typename CalculationType, typename Strategy>
 322     static inline void apply(CalculationType& lon1,
 323                              CalculationType& lat1,
 324                              CalculationType& lon2,
 325                              CalculationType& lat2,
 326                              Strategy const& strategy)
 327     {
 328         special_cases<Units>(lon1, lat1, lon2, lat2);
 329
 330         CalculationType lon1_rad = math::as_radian<Units>(lon1);
 331         CalculationType lat1_rad = math::as_radian<Units>(lat1);
 332         CalculationType lon2_rad = math::as_radian<Units>(lon2);
 333         CalculationType lat2_rad = math::as_radian<Units>(lat2);
 334         CalculationType alp1, alp2;
 335         strategy.apply(lon1_rad, lat1_rad, lon2_rad, lat2_rad, alp1, alp2);
 336
 337         compute_box_corners<Units>(lon1, lat1, lon2, lat2, alp1, alp2, strategy);
 338     }
 339
 340 public:
 341     template
 342     <
 343         typename Units,
 344         typename CalculationType,
 345         typename Box,
 346         typename Strategy
 347     >
 348     static inline void apply(CalculationType lon1,
 349                              CalculationType lat1,
 350                              CalculationType lon2,
 351                              CalculationType lat2,
 352                              Box& mbr,
 353                              Strategy const& strategy)
 354     {
 355         typedef envelope_segment_convert_polar<Units, typename cs_tag<Box>::type> convert_polar;
 356
 357         convert_polar::pre(lat1, lat2);
 358
 359         apply<Units>(lon1, lat1, lon2, lat2, strategy);
 360
 361         convert_polar::post(lat1, lat2);
 362
 363         create_box<Units>(lon1, lat1, lon2, lat2, mbr);
 364     }
 365
 366 };
 367
 368 } // namespace detail
 369 #endif // DOXYGEN_NO_DETAIL
 370
 371
 372 template
 373 <
 374     typename CalculationType = void
 375 >
 376 class spherical_segment
 377 {
 378 public:
 379     typedef strategy::expand::spherical_box box_expand_strategy_type;
 380     static inline box_expand_strategy_type get_box_expand_strategy()
 381     {
 382         return box_expand_strategy_type();
 383     }
 384
 385     template <typename Point, typename Box>
 386     static inline void apply(Point const& point1, Point const& point2,
 387                              Box& box)
 388     {
 389         Point p1_normalized, p2_normalized;
 390         strategy::normalize::spherical_point::apply(point1, p1_normalized);
 391         strategy::normalize::spherical_point::apply(point2, p2_normalized);
 392
 393         geometry::strategy::azimuth::spherical<CalculationType> azimuth_spherical;
 394
 395         typedef typename geometry::detail::cs_angular_units<Point>::type units_type;
 396
 397         // first compute the envelope range for the first two coordinates
 398         strategy::envelope::detail::envelope_segment_impl
 399             <
 400                 spherical_equatorial_tag
 401             >::template apply<units_type>(geometry::get<0>(p1_normalized),
 402                                           geometry::get<1>(p1_normalized),
 403                                           geometry::get<0>(p2_normalized),
 404                                           geometry::get<1>(p2_normalized),
 405                                           box,
 406                                           azimuth_spherical);
 407
 408         // now compute the envelope range for coordinates of
 409         // dimension 2 and higher
 410         strategy::envelope::detail::envelope_one_segment
 411             <
 412                 2, dimension<Point>::value
 413             >::apply(point1, point2, box);
 414   }
 415 };
 416
 417 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 418
 419 namespace services
 420 {
 421
 422 template <typename CalculationType>
 423 struct default_strategy<segment_tag, spherical_equatorial_tag, CalculationType>
 424 {
 425     typedef strategy::envelope::spherical_segment<CalculationType> type;
 426 };
 427
 428
 429 template <typename CalculationType>
 430 struct default_strategy<segment_tag, spherical_polar_tag, CalculationType>
 431 {
 432     typedef strategy::envelope::spherical_segment<CalculationType> type;
 433 };
 434
 435 }
 436
 437 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 438
 439
 440 }} // namespace strategy::envelope
 441
 442 }} //namepsace boost::geometry
 443
 444 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_SEGMENT_HPP
 445