strategy/spherical/envelope_segment.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2017-2020 Oracle and/or its affiliates.
0004 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
0005 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0006
0007 // Use, modification and distribution is subject to the Boost Software License,
0008 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0009 // http://www.boost.org/LICENSE_1_0.txt)
0010
0011 #ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_SEGMENT_HPP
0012 #define BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_SEGMENT_HPP
0013
0014
0015 #include <cstddef>
0016 #include <utility>
0017
0018 #include <boost/core/ignore_unused.hpp>
0019 #include <boost/numeric/conversion/cast.hpp>
0020
0021 #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp>
0022
0023 #include <boost/geometry/core/assert.hpp>
0024 #include <boost/geometry/core/coordinate_system.hpp>
0025 #include <boost/geometry/core/coordinate_type.hpp>
0026 #include <boost/geometry/core/cs.hpp>
0027 #include <boost/geometry/core/point_type.hpp>
0028 #include <boost/geometry/core/radian_access.hpp>
0029 #include <boost/geometry/core/tags.hpp>
0030
0031 #include <boost/geometry/formulas/meridian_segment.hpp>
0032 #include <boost/geometry/formulas/vertex_latitude.hpp>
0033
0034 #include <boost/geometry/geometries/helper_geometry.hpp>
0035
0036 #include <boost/geometry/strategy/cartesian/envelope_segment.hpp>
0037 #include <boost/geometry/strategy/envelope.hpp>
0038 #include <boost/geometry/strategies/normalize.hpp>
0039 #include <boost/geometry/strategies/spherical/azimuth.hpp>
0040 #include <boost/geometry/strategy/spherical/expand_box.hpp>
0041
0042 #include <boost/geometry/util/math.hpp>
0043
0044 namespace boost { namespace geometry { namespace strategy { namespace envelope
0045 {
0046
0047 #ifndef DOXYGEN_NO_DETAIL
0048 namespace detail
0049 {
0050
0051 template <typename CalculationType, typename CS_Tag>
0052 struct envelope_segment_call_vertex_latitude
0053 {
0054     template <typename T1, typename T2, typename Strategy>
0055     static inline CalculationType apply(T1 const& lat1,
0056                                         T2 const& alp1,
0057                                         Strategy const& )
0058     {
0059         return geometry::formula::vertex_latitude<CalculationType, CS_Tag>
0060             ::apply(lat1, alp1);
0061     }
0062 };
0063
0064 template <typename CalculationType>
0065 struct envelope_segment_call_vertex_latitude<CalculationType, geographic_tag>
0066 {
0067     template <typename T1, typename T2, typename Strategy>
0068     static inline CalculationType apply(T1 const& lat1,
0069                                         T2 const& alp1,
0070                                         Strategy const& strategy)
0071     {
0072         return geometry::formula::vertex_latitude<CalculationType, geographic_tag>
0073             ::apply(lat1, alp1, strategy.model());
0074     }
0075 };
0076
0077 template <typename Units, typename CS_Tag>
0078 struct envelope_segment_convert_polar
0079 {
0080     template <typename T>
0081     static inline void pre(T & , T & ) {}
0082
0083     template <typename T>
0084     static inline void post(T & , T & ) {}
0085 };
0086
0087 template <typename Units>
0088 struct envelope_segment_convert_polar<Units, spherical_polar_tag>
0089 {
0090     template <typename T>
0091     static inline void pre(T & lat1, T & lat2)
0092     {
0093         lat1 = math::latitude_convert_ep<Units>(lat1);
0094         lat2 = math::latitude_convert_ep<Units>(lat2);
0095     }
0096
0097     template <typename T>
0098     static inline void post(T & lat1, T & lat2)
0099     {
0100         lat1 = math::latitude_convert_ep<Units>(lat1);
0101         lat2 = math::latitude_convert_ep<Units>(lat2);
0102         std::swap(lat1, lat2);
0103     }
0104 };
0105
0106 template <typename CS_Tag>
0107 class envelope_segment_impl
0108 {
0109 private:
0110
0111     // degrees or radians
0112     template <typename CalculationType>
0113     static inline void swap(CalculationType& lon1,
0114                             CalculationType& lat1,
0115                             CalculationType& lon2,
0116                             CalculationType& lat2)
0117     {
0118         std::swap(lon1, lon2);
0119         std::swap(lat1, lat2);
0120     }
0121
0122     // radians
0123     template <typename CalculationType>
0124     static inline bool contains_pi_half(CalculationType const& a1,
0125                                         CalculationType const& a2)
0126     {
0127         // azimuths a1 and a2 are assumed to be in radians
0128
0129         static CalculationType const pi_half = math::half_pi<CalculationType>();
0130
0131         return (a1 < a2)
0132                 ? (a1 < pi_half && pi_half < a2)
0133                 : (a1 > pi_half && pi_half > a2);
0134     }
0135
0136     // radians or degrees
0137     template <typename Units, typename CoordinateType>
0138     static inline bool crosses_antimeridian(CoordinateType const& lon1,
0139                                             CoordinateType const& lon2)
0140     {
0141         typedef math::detail::constants_on_spheroid
0142             <
0143                 CoordinateType, Units
0144             > constants;
0145
0146         return math::abs(lon1 - lon2) > constants::half_period(); // > pi
0147     }
0148
0149     // degrees or radians
0150     template <typename Units, typename CalculationType, typename Strategy>
0151     static inline void compute_box_corners(CalculationType& lon1,
0152                                            CalculationType& lat1,
0153                                            CalculationType& lon2,
0154                                            CalculationType& lat2,
0155                                            CalculationType a1,
0156                                            CalculationType a2,
0157                                            Strategy const& strategy)
0158     {
0159         // coordinates are assumed to be in radians
0160         BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
0161         boost::ignore_unused(lon1, lon2);
0162
0163         CalculationType lat1_rad = math::as_radian<Units>(lat1);
0164         CalculationType lat2_rad = math::as_radian<Units>(lat2);
0165
0166         if (lat1 > lat2)
0167         {
0168             std::swap(lat1, lat2);
0169             std::swap(lat1_rad, lat2_rad);
0170             std::swap(a1, a2);
0171         }
0172
0173         if (contains_pi_half(a1, a2))
0174         {
0175             CalculationType p_max = envelope_segment_call_vertex_latitude
0176                 <CalculationType, CS_Tag>::apply(lat1_rad, a1, strategy);
0177
0178             CalculationType const mid_lat = lat1 + lat2;
0179             if (mid_lat < 0)
0180             {
0181                 // update using min latitude
0182                 CalculationType const lat_min_rad = -p_max;
0183                 CalculationType const lat_min
0184                     = math::from_radian<Units>(lat_min_rad);
0185
0186                 if (lat1 > lat_min)
0187                 {
0188                     lat1 = lat_min;
0189                 }
0190             }
0191             else
0192             {
0193                 // update using max latitude
0194                 CalculationType const lat_max_rad = p_max;
0195                 CalculationType const lat_max
0196                     = math::from_radian<Units>(lat_max_rad);
0197
0198                 if (lat2 < lat_max)
0199                 {
0200                     lat2 = lat_max;
0201                 }
0202             }
0203         }
0204     }
0205
0206     template <typename Units, typename CalculationType>
0207     static inline void special_cases(CalculationType& lon1,
0208                                      CalculationType& lat1,
0209                                      CalculationType& lon2,
0210                                      CalculationType& lat2)
0211     {
0212         typedef math::detail::constants_on_spheroid
0213             <
0214                 CalculationType, Units
0215             > constants;
0216
0217         bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
0218         bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());
0219
0220         if (is_pole1 && is_pole2)
0221         {
0222             // both points are poles; nothing more to do:
0223             // longitudes are already normalized to 0
0224             // but just in case
0225             lon1 = 0;
0226             lon2 = 0;
0227         }
0228         else if (is_pole1 && !is_pole2)
0229         {
0230             // first point is a pole, second point is not:
0231             // make the longitude of the first point the same as that
0232             // of the second point
0233             lon1 = lon2;
0234         }
0235         else if (!is_pole1 && is_pole2)
0236         {
0237             // second point is a pole, first point is not:
0238             // make the longitude of the second point the same as that
0239             // of the first point
0240             lon2 = lon1;
0241         }
0242
0243         if (lon1 == lon2)
0244         {
0245             // segment lies on a meridian
0246             if (lat1 > lat2)
0247             {
0248                 std::swap(lat1, lat2);
0249             }
0250             return;
0251         }
0252
0253         BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);
0254
0255         if (lon1 > lon2)
0256         {
0257             swap(lon1, lat1, lon2, lat2);
0258         }
0259
0260         if (crosses_antimeridian<Units>(lon1, lon2))
0261         {
0262             lon1 += constants::period();
0263             swap(lon1, lat1, lon2, lat2);
0264         }
0265     }
0266
0267     template
0268     <
0269         typename Units,
0270         typename CalculationType,
0271         typename Box
0272     >
0273     static inline void create_box(CalculationType lon1,
0274                                   CalculationType lat1,
0275                                   CalculationType lon2,
0276                                   CalculationType lat2,
0277                                   Box& mbr)
0278     {
0279         typedef typename coordinate_type<Box>::type box_coordinate_type;
0280
0281         typedef typename helper_geometry
0282             <
0283                 Box, box_coordinate_type, Units
0284             >::type helper_box_type;
0285
0286         helper_box_type helper_mbr;
0287
0288         geometry::set
0289             <
0290                 min_corner, 0
0291             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lon1));
0292
0293         geometry::set
0294             <
0295                 min_corner, 1
0296             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lat1));
0297
0298         geometry::set
0299             <
0300                 max_corner, 0
0301             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lon2));
0302
0303         geometry::set
0304             <
0305                 max_corner, 1
0306             >(helper_mbr, boost::numeric_cast<box_coordinate_type>(lat2));
0307
0308         geometry::detail::envelope::transform_units(helper_mbr, mbr);
0309     }
0310
0311
0312     template <typename Units, typename CalculationType, typename Strategy>
0313     static inline void apply(CalculationType& lon1,
0314                              CalculationType& lat1,
0315                              CalculationType& lon2,
0316                              CalculationType& lat2,
0317                              Strategy const& strategy)
0318     {
0319         special_cases<Units>(lon1, lat1, lon2, lat2);
0320
0321         CalculationType lon1_rad = math::as_radian<Units>(lon1);
0322         CalculationType lat1_rad = math::as_radian<Units>(lat1);
0323         CalculationType lon2_rad = math::as_radian<Units>(lon2);
0324         CalculationType lat2_rad = math::as_radian<Units>(lat2);
0325         CalculationType alp1, alp2;
0326         strategy.apply(lon1_rad, lat1_rad, lon2_rad, lat2_rad, alp1, alp2);
0327
0328         compute_box_corners<Units>(lon1, lat1, lon2, lat2, alp1, alp2, strategy);
0329     }
0330
0331 public:
0332     template
0333     <
0334         typename Units,
0335         typename CalculationType,
0336         typename Box,
0337         typename Strategy
0338     >
0339     static inline void apply(CalculationType lon1,
0340                              CalculationType lat1,
0341                              CalculationType lon2,
0342                              CalculationType lat2,
0343                              Box& mbr,
0344                              Strategy const& strategy)
0345     {
0346         typedef envelope_segment_convert_polar<Units, typename cs_tag<Box>::type> convert_polar;
0347
0348         convert_polar::pre(lat1, lat2);
0349
0350         apply<Units>(lon1, lat1, lon2, lat2, strategy);
0351
0352         convert_polar::post(lat1, lat2);
0353
0354         create_box<Units>(lon1, lat1, lon2, lat2, mbr);
0355     }
0356
0357 };
0358
0359 } // namespace detail
0360 #endif // DOXYGEN_NO_DETAIL
0361
0362
0363 template
0364 <
0365     typename CalculationType = void
0366 >
0367 class spherical_segment
0368 {
0369 public:
0370     template <typename Point, typename Box>
0371     static inline void apply(Point const& point1, Point const& point2,
0372                              Box& box)
0373     {
0374         Point p1_normalized, p2_normalized;
0375         strategy::normalize::spherical_point::apply(point1, p1_normalized);
0376         strategy::normalize::spherical_point::apply(point2, p2_normalized);
0377
0378         geometry::strategy::azimuth::spherical<CalculationType> azimuth_spherical;
0379
0380         typedef typename geometry::detail::cs_angular_units<Point>::type units_type;
0381
0382         // first compute the envelope range for the first two coordinates
0383         strategy::envelope::detail::envelope_segment_impl
0384             <
0385                 spherical_equatorial_tag
0386             >::template apply<units_type>(geometry::get<0>(p1_normalized),
0387                                           geometry::get<1>(p1_normalized),
0388                                           geometry::get<0>(p2_normalized),
0389                                           geometry::get<1>(p2_normalized),
0390                                           box,
0391                                           azimuth_spherical);
0392
0393         // now compute the envelope range for coordinates of
0394         // dimension 2 and higher
0395         strategy::envelope::detail::envelope_one_segment
0396             <
0397                 2, dimension<Point>::value
0398             >::apply(point1, point2, box);
0399   }
0400 };
0401
0402 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0403
0404 namespace services
0405 {
0406
0407 template <typename CalculationType>
0408 struct default_strategy<segment_tag, spherical_equatorial_tag, CalculationType>
0409 {
0410     typedef strategy::envelope::spherical_segment<CalculationType> type;
0411 };
0412
0413
0414 template <typename CalculationType>
0415 struct default_strategy<segment_tag, spherical_polar_tag, CalculationType>
0416 {
0417     typedef strategy::envelope::spherical_segment<CalculationType> type;
0418 };
0419
0420 }
0421
0422 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0423
0424
0425 }} // namespace strategy::envelope
0426
0427 }} //namepsace boost::geometry
0428
0429 #endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_SEGMENT_HPP
0430