strategy/spherical/envelope_multipoint.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2015-2023, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0007 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0009
0010 // Distributed under the Boost Software License, Version 1.0.
0011 // (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 #ifndef BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_MULTIPOINT_HPP
0015 #define BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_MULTIPOINT_HPP
0016
0017 #include <algorithm>
0018 #include <cstddef>
0019 #include <utility>
0020 #include <vector>
0021
0022 #include <boost/range/begin.hpp>
0023 #include <boost/range/empty.hpp>
0024 #include <boost/range/end.hpp>
0025 #include <boost/range/size.hpp>
0026 #include <boost/range/value_type.hpp>
0027
0028 #include <boost/geometry/core/access.hpp>
0029 #include <boost/geometry/core/assert.hpp>
0030 #include <boost/geometry/core/coordinate_system.hpp>
0031 #include <boost/geometry/core/coordinate_type.hpp>
0032 #include <boost/geometry/core/tags.hpp>
0033
0034 #include <boost/geometry/util/math.hpp>
0035 #include <boost/geometry/util/range.hpp>
0036
0037 #include <boost/geometry/geometries/helper_geometry.hpp>
0038
0039 #include <boost/geometry/algorithms/detail/envelope/box.hpp>
0040 #include <boost/geometry/algorithms/detail/envelope/initialize.hpp>
0041 #include <boost/geometry/algorithms/detail/envelope/range.hpp>
0042 #include <boost/geometry/algorithms/detail/expand/point.hpp>
0043
0044 #include <boost/geometry/strategy/cartesian/envelope_point.hpp>
0045 #include <boost/geometry/strategy/cartesian/expand_point.hpp>
0046 #include <boost/geometry/strategies/normalize.hpp>
0047 #include <boost/geometry/strategy/spherical/envelope_box.hpp>
0048 #include <boost/geometry/strategy/spherical/envelope_point.hpp>
0049
0050
0051 namespace boost { namespace geometry
0052 {
0053
0054 namespace strategy { namespace envelope
0055 {
0056
0057 class spherical_multipoint
0058 {
0059 private:
0060     template <std::size_t Dim>
0061     struct coordinate_less
0062     {
0063         template <typename Point>
0064         inline bool operator()(Point const& point1, Point const& point2) const
0065         {
0066             return math::smaller(geometry::get<Dim>(point1), geometry::get<Dim>(point2));
0067         }
0068     };
0069
0070     template <typename Constants, typename MultiPoint, typename OutputIterator>
0071     static inline void analyze_point_coordinates(MultiPoint const& multipoint,
0072                                                  bool& has_south_pole,
0073                                                  bool& has_north_pole,
0074                                                  OutputIterator oit)
0075     {
0076         // analyze point coordinates:
0077         // (1) normalize point coordinates
0078         // (2) check if any point is the north or the south pole
0079         // (3) put all non-pole points in a container
0080         //
0081         // notice that at this point in the algorithm, we have at
0082         // least two points on the spheroid
0083         has_south_pole = false;
0084         has_north_pole = false;
0085
0086         for (auto it = boost::begin(multipoint); it != boost::end(multipoint); ++it)
0087         {
0088             typename boost::range_value<MultiPoint>::type point;
0089             normalize::spherical_point::apply(*it, point);
0090
0091             if (math::equals(geometry::get<1>(point), Constants::min_latitude()))
0092             {
0093                 has_south_pole = true;
0094             }
0095             else if (math::equals(geometry::get<1>(point), Constants::max_latitude()))
0096             {
0097                 has_north_pole = true;
0098             }
0099             else
0100             {
0101                 *oit++ = point;
0102             }
0103         }
0104     }
0105
0106     template <typename SortedRange, typename Value>
0107     static inline Value maximum_gap(SortedRange const& sorted_range,
0108                                     Value& max_gap_left,
0109                                     Value& max_gap_right)
0110     {
0111         auto it1 = boost::begin(sorted_range);
0112         auto it2 = it1;
0113         ++it2;
0114         max_gap_left = geometry::get<0>(*it1);
0115         max_gap_right = geometry::get<0>(*it2);
0116
0117         Value max_gap = max_gap_right - max_gap_left;
0118         for (++it1, ++it2; it2 != boost::end(sorted_range); ++it1, ++it2)
0119         {
0120             Value gap = geometry::get<0>(*it2) - geometry::get<0>(*it1);
0121             if (math::larger(gap, max_gap))
0122             {
0123                 max_gap_left = geometry::get<0>(*it1);
0124                 max_gap_right = geometry::get<0>(*it2);
0125                 max_gap = gap;
0126             }
0127         }
0128
0129         return max_gap;
0130     }
0131
0132     template
0133     <
0134         typename Constants,
0135         typename PointRange,
0136         typename LongitudeLess,
0137         typename CoordinateType
0138     >
0139     static inline void get_min_max_longitudes(PointRange& range,
0140                                               LongitudeLess const& lon_less,
0141                                               CoordinateType& lon_min,
0142                                               CoordinateType& lon_max)
0143     {
0144         // compute min and max longitude values
0145         auto const min_max_longitudes
0146             = std::minmax_element(boost::begin(range), boost::end(range), lon_less);
0147
0148         lon_min = geometry::get<0>(*min_max_longitudes.first);
0149         lon_max = geometry::get<0>(*min_max_longitudes.second);
0150
0151         // if the longitude span is "large" compute the true maximum gap
0152         if (math::larger(lon_max - lon_min, Constants::half_period()))
0153         {
0154             std::sort(boost::begin(range), boost::end(range), lon_less);
0155
0156             CoordinateType max_gap_left = 0, max_gap_right = 0;
0157             CoordinateType max_gap
0158                 = maximum_gap(range, max_gap_left, max_gap_right);
0159
0160             CoordinateType complement_gap
0161                 = Constants::period() + lon_min - lon_max;
0162
0163             if (math::larger(max_gap, complement_gap))
0164             {
0165                 lon_min = max_gap_right;
0166                 lon_max = max_gap_left + Constants::period();
0167             }
0168         }
0169     }
0170
0171     template
0172     <
0173         typename Constants,
0174         typename Iterator,
0175         typename LatitudeLess,
0176         typename CoordinateType
0177     >
0178     static inline void get_min_max_latitudes(Iterator const first,
0179                                              Iterator const last,
0180                                              LatitudeLess const& lat_less,
0181                                              bool has_south_pole,
0182                                              bool has_north_pole,
0183                                              CoordinateType& lat_min,
0184                                              CoordinateType& lat_max)
0185     {
0186         if (has_south_pole && has_north_pole)
0187         {
0188             lat_min = Constants::min_latitude();
0189             lat_max = Constants::max_latitude();
0190         }
0191         else if (has_south_pole)
0192         {
0193             lat_min = Constants::min_latitude();
0194             lat_max = geometry::get<1>(*std::max_element(first, last, lat_less));
0195         }
0196         else if (has_north_pole)
0197         {
0198             lat_min = geometry::get<1>(*std::min_element(first, last, lat_less));
0199             lat_max = Constants::max_latitude();
0200         }
0201         else
0202         {
0203             auto const min_max_latitudes = std::minmax_element(first, last, lat_less);
0204             lat_min = geometry::get<1>(*min_max_latitudes.first);
0205             lat_max = geometry::get<1>(*min_max_latitudes.second);
0206         }
0207     }
0208
0209 public:
0210     template <typename MultiPoint, typename Box>
0211     static inline void apply(MultiPoint const& multipoint, Box& mbr)
0212     {
0213         typedef typename point_type<MultiPoint>::type point_type;
0214         typedef typename coordinate_type<MultiPoint>::type coordinate_type;
0215         typedef math::detail::constants_on_spheroid
0216             <
0217                 coordinate_type,
0218                 typename geometry::detail::cs_angular_units<MultiPoint>::type
0219             > constants;
0220
0221         if (boost::empty(multipoint))
0222         {
0223             geometry::detail::envelope::initialize<Box, 0, dimension<Box>::value>::apply(mbr);
0224             return;
0225         }
0226
0227         geometry::detail::envelope::initialize<Box, 0, 2>::apply(mbr);
0228
0229         if (boost::size(multipoint) == 1)
0230         {
0231             spherical_point::apply(range::front(multipoint), mbr);
0232             return;
0233         }
0234
0235         // analyze the points and put the non-pole ones in the
0236         // points vector
0237         std::vector<point_type> points;
0238         bool has_north_pole = false, has_south_pole = false;
0239
0240         analyze_point_coordinates<constants>(multipoint,
0241                                              has_south_pole, has_north_pole,
0242                                              std::back_inserter(points));
0243
0244         coordinate_type lon_min, lat_min, lon_max, lat_max;
0245         if (points.size() == 1)
0246         {
0247             // we have one non-pole point and at least one pole point
0248             lon_min = geometry::get<0>(range::front(points));
0249             lon_max = geometry::get<0>(range::front(points));
0250             lat_min = has_south_pole
0251                 ? constants::min_latitude()
0252                 : constants::max_latitude();
0253             lat_max = has_north_pole
0254                 ? constants::max_latitude()
0255                 : constants::min_latitude();
0256         }
0257         else if (points.empty())
0258         {
0259             // all points are pole points
0260             BOOST_GEOMETRY_ASSERT(has_south_pole || has_north_pole);
0261             lon_min = coordinate_type(0);
0262             lon_max = coordinate_type(0);
0263             lat_min = has_south_pole
0264                 ? constants::min_latitude()
0265                 : constants::max_latitude();
0266             lat_max = (has_north_pole)
0267                 ? constants::max_latitude()
0268                 : constants::min_latitude();
0269         }
0270         else
0271         {
0272             get_min_max_longitudes<constants>(points,
0273                                               coordinate_less<0>(),
0274                                               lon_min,
0275                                               lon_max);
0276
0277             get_min_max_latitudes<constants>(points.begin(),
0278                                              points.end(),
0279                                              coordinate_less<1>(),
0280                                              has_south_pole,
0281                                              has_north_pole,
0282                                              lat_min,
0283                                              lat_max);
0284         }
0285
0286         typedef typename helper_geometry
0287             <
0288                 Box,
0289                 coordinate_type,
0290                 typename geometry::detail::cs_angular_units<MultiPoint>::type
0291             >::type helper_box_type;
0292
0293         helper_box_type helper_mbr;
0294
0295         geometry::set<min_corner, 0>(helper_mbr, lon_min);
0296         geometry::set<min_corner, 1>(helper_mbr, lat_min);
0297         geometry::set<max_corner, 0>(helper_mbr, lon_max);
0298         geometry::set<max_corner, 1>(helper_mbr, lat_max);
0299
0300         // now transform to output MBR (per index)
0301         geometry::detail::envelope::envelope_indexed_box_on_spheroid<min_corner, 2>::apply(helper_mbr, mbr);
0302         geometry::detail::envelope::envelope_indexed_box_on_spheroid<max_corner, 2>::apply(helper_mbr, mbr);
0303
0304         // compute envelope for higher coordinates
0305         auto it = boost::begin(multipoint);
0306         geometry::detail::envelope::envelope_one_point<2, dimension<Box>::value>::apply(*it, mbr);
0307
0308         for (++it; it != boost::end(multipoint); ++it)
0309         {
0310             strategy::expand::detail::point_loop
0311                 <
0312                     2, dimension<Box>::value
0313                 >::apply(mbr, *it);
0314         }
0315     }
0316 };
0317
0318
0319 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0320
0321 namespace services
0322 {
0323
0324 template <typename CalculationType>
0325 struct default_strategy<multi_point_tag, spherical_equatorial_tag, CalculationType>
0326 {
0327     typedef strategy::envelope::spherical_multipoint type;
0328 };
0329
0330 template <typename CalculationType>
0331 struct default_strategy<multi_point_tag, spherical_polar_tag, CalculationType>
0332 {
0333     typedef strategy::envelope::spherical_multipoint type;
0334 };
0335
0336 template <typename CalculationType>
0337 struct default_strategy<multi_point_tag, geographic_tag, CalculationType>
0338 {
0339     typedef strategy::envelope::spherical_multipoint type;
0340 };
0341
0342 } // namespace services
0343
0344 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0345
0346
0347 }} // namespace strategy::envelope
0348
0349 }} // namespace boost::geometry
0350
0351 #endif // BOOST_GEOMETRY_STRATEGY_SPHERICAL_ENVELOPE_MULTIPOINT_HPP