strategy/spherical/envelope_range.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2021-2022, Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0006
0007 // Licensed under the Boost Software License version 1.0.
0008 // http://www.boost.org/users/license.html
0009
0010 #ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_RANGE_HPP
0011 #define BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_RANGE_HPP
0012
0013 #include <boost/range/size.hpp>
0014
0015 #include <boost/geometry/algorithms/assign.hpp>
0016 #include <boost/geometry/algorithms/detail/envelope/initialize.hpp>
0017 #include <boost/geometry/geometries/segment.hpp>
0018 #include <boost/geometry/strategy/spherical/envelope_point.hpp>
0019 #include <boost/geometry/strategy/spherical/envelope_segment.hpp>
0020 #include <boost/geometry/strategy/spherical/expand_segment.hpp>
0021 #include <boost/geometry/views/closeable_view.hpp>
0022
0023 // Get rid of this dependency?
0024 #include <boost/geometry/strategies/spherical/point_in_poly_winding.hpp>
0025
0026
0027 namespace boost { namespace geometry
0028 {
0029
0030 namespace strategy { namespace envelope
0031 {
0032
0033 #ifndef DOXYGEN_NO_DETAIL
0034 namespace detail
0035 {
0036
0037 template <typename Range, typename Box, typename EnvelopeStrategy, typename ExpandStrategy>
0038 inline void spheroidal_linestring(Range const& range, Box& mbr,
0039                                   EnvelopeStrategy const& envelope_strategy,
0040                                   ExpandStrategy const& expand_strategy)
0041 {
0042     auto it = boost::begin(range);
0043     auto const end = boost::end(range);
0044     if (it == end)
0045     {
0046         // initialize box (assign inverse)
0047         geometry::detail::envelope::initialize<Box>::apply(mbr);
0048         return;
0049     }
0050
0051     auto prev = it;
0052     ++it;
0053     if (it == end)
0054     {
0055         // initialize box with the first point
0056         envelope::spherical_point::apply(*prev, mbr);
0057         return;
0058     }
0059
0060     // initialize box with the first segment
0061     envelope_strategy.apply(*prev, *it, mbr);
0062
0063     // consider now the remaining segments in the range (if any)
0064     prev = it;
0065     ++it;
0066     while (it != end)
0067     {
0068         using point_t = typename boost::range_value<Range>::type;
0069         geometry::model::referring_segment<point_t const> const seg(*prev, *it);
0070         expand_strategy.apply(mbr, seg);
0071         prev = it;
0072         ++it;
0073     }
0074 }
0075
0076
0077 // This strategy is intended to be used together with winding strategy to check
0078 // if ring/polygon has a pole in its interior or exterior. It is not intended
0079 // for checking if the pole is on the boundary.
0080 template <typename CalculationType = void>
0081 struct side_of_pole
0082 {
0083     typedef spherical_tag cs_tag;
0084
0085     template <typename P>
0086     static inline int apply(P const& p1, P const& p2, P const& pole)
0087     {
0088         using calc_t = typename promote_floating_point
0089             <
0090                 typename select_calculation_type_alt
0091                     <
0092                         CalculationType, P
0093                     >::type
0094             >::type;
0095
0096         using units_t = typename geometry::detail::cs_angular_units<P>::type;
0097         using constants = math::detail::constants_on_spheroid<calc_t, units_t>;
0098
0099         calc_t const c0 = 0;
0100         calc_t const pi = constants::half_period();
0101
0102         calc_t const lon1 = get<0>(p1);
0103         calc_t const lat1 = get<1>(p1);
0104         calc_t const lon2 = get<0>(p2);
0105         calc_t const lat2 = get<1>(p2);
0106         calc_t const lat_pole = get<1>(pole);
0107
0108         calc_t const s_lon_diff = math::longitude_distance_signed<units_t>(lon1, lon2);
0109         bool const s_vertical = math::equals(s_lon_diff, c0)
0110                              || math::equals(s_lon_diff, pi);
0111         // Side of vertical segment is 0 for both poles.
0112         if (s_vertical)
0113         {
0114             return 0;
0115         }
0116
0117         // This strategy shouldn't be called in this case but just in case
0118         // check if segment starts at a pole
0119         if (math::equals(lat_pole, lat1) || math::equals(lat_pole, lat2))
0120         {
0121             return 0;
0122         }
0123
0124         // -1 is rhs
0125         //  1 is lhs
0126         if (lat_pole >= c0) // north pole
0127         {
0128             return s_lon_diff < c0 ? -1 : 1;
0129         }
0130         else // south pole
0131         {
0132             return s_lon_diff > c0 ? -1 : 1;
0133         }
0134     }
0135 };
0136
0137
0138 template <typename Point, typename Range, typename Strategy>
0139 inline int point_in_range(Point const& point, Range const& range, Strategy const& strategy)
0140 {
0141     typename Strategy::state_type state;
0142
0143     auto it = boost::begin(range);
0144     auto const end = boost::end(range);
0145     for (auto previous = it++ ; it != end ; ++previous, ++it )
0146     {
0147         if (! strategy.apply(point, *previous, *it, state))
0148         {
0149             break;
0150         }
0151     }
0152
0153     return strategy.result(state);
0154 }
0155
0156
0157 template <typename T, typename Ring, typename PoleWithinStrategy>
0158 inline bool pole_within(T const& lat_pole, Ring const& ring,
0159                         PoleWithinStrategy const& pole_within_strategy)
0160 {
0161     if (boost::size(ring) < core_detail::closure::minimum_ring_size
0162                                 <
0163                                     geometry::closure<Ring>::value
0164                                 >::value)
0165     {
0166         return false;
0167     }
0168
0169     using point_t = typename geometry::point_type<Ring>::type;
0170     point_t point;
0171     geometry::assign_zero(point);
0172     geometry::set<1>(point, lat_pole);
0173     geometry::detail::closed_clockwise_view<Ring const> view(ring);
0174     return point_in_range(point, view, pole_within_strategy) > 0;
0175 }
0176
0177 template
0178 <
0179     typename Range,
0180     typename Box,
0181     typename EnvelopeStrategy,
0182     typename ExpandStrategy,
0183     typename PoleWithinStrategy
0184 >
0185 inline void spheroidal_ring(Range const& range, Box& mbr,
0186                             EnvelopeStrategy const& envelope_strategy,
0187                             ExpandStrategy const& expand_strategy,
0188                             PoleWithinStrategy const& pole_within_strategy)
0189 {
0190     geometry::detail::closed_view<Range const> closed_range(range);
0191
0192     spheroidal_linestring(closed_range, mbr, envelope_strategy, expand_strategy);
0193
0194     using coord_t = typename geometry::coordinate_type<Box>::type;
0195     using point_t = typename geometry::point_type<Box>::type;
0196     using units_t = typename geometry::detail::cs_angular_units<point_t>::type;
0197     using constants_t = math::detail::constants_on_spheroid<coord_t, units_t>;
0198     coord_t const two_pi = constants_t::period();
0199     coord_t const lon_min = geometry::get<0, 0>(mbr);
0200     coord_t const lon_max = geometry::get<1, 0>(mbr);
0201     // If box covers the whole longitude range it is possible that the ring contains
0202     // one of the poles.
0203     // Technically it is possible that a reversed ring may cover more than
0204     // half of the globe and mbr of it's linear ring may be small and not cover the
0205     // longitude range. We currently don't support such rings.
0206     if (lon_max - lon_min >= two_pi)
0207     {
0208         coord_t const lat_n_pole = constants_t::max_latitude();
0209         coord_t const lat_s_pole = constants_t::min_latitude();
0210         coord_t lat_min = geometry::get<0, 1>(mbr);
0211         coord_t lat_max = geometry::get<1, 1>(mbr);
0212         // Normalize box latitudes, just in case
0213         if (math::equals(lat_min, lat_s_pole))
0214         {
0215             lat_min = lat_s_pole;
0216         }
0217         if (math::equals(lat_max, lat_n_pole))
0218         {
0219             lat_max = lat_n_pole;
0220         }
0221
0222         if (lat_max < lat_n_pole)
0223         {
0224             if (pole_within(lat_n_pole, range, pole_within_strategy))
0225             {
0226                 lat_max = lat_n_pole;
0227             }
0228         }
0229         if (lat_min > lat_s_pole)
0230         {
0231             if (pole_within(lat_s_pole, range, pole_within_strategy))
0232             {
0233                 lat_min = lat_s_pole;
0234             }
0235         }
0236
0237         geometry::set<0, 1>(mbr, lat_min);
0238         geometry::set<1, 1>(mbr, lat_max);
0239     }
0240 }
0241
0242 } // namespace detail
0243 #endif // DOXYGEN_NO_DETAIL
0244
0245 template <typename CalculationType = void>
0246 class spherical_linestring
0247 {
0248 public:
0249     template <typename Range, typename Box>
0250     static inline void apply(Range const& range, Box& mbr)
0251     {
0252         detail::spheroidal_linestring(range, mbr,
0253                                       envelope::spherical_segment<CalculationType>(),
0254                                       expand::spherical_segment<CalculationType>());
0255     }
0256 };
0257
0258 template <typename CalculationType = void>
0259 class spherical_ring
0260 {
0261 public:
0262     template <typename Range, typename Box>
0263     static inline void apply(Range const& range, Box& mbr)
0264     {
0265         detail::spheroidal_ring(range, mbr,
0266                                 envelope::spherical_segment<CalculationType>(),
0267                                 expand::spherical_segment<CalculationType>(),
0268                                 within::detail::spherical_winding_base
0269                                     <
0270                                         envelope::detail::side_of_pole<CalculationType>,
0271                                         CalculationType
0272                                     >());
0273     }
0274 };
0275
0276 }} // namespace strategy::envelope
0277
0278 }} //namepsace boost::geometry
0279
0280 #endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_RANGE_HPP