1 // Boost.Geometry (aka GGL, Generic Geometry Library) 2 3 // Copyright (c) 2015-2018, Oracle and/or its affiliates. 4 5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle 6 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle 7 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle 8 9 // Distributed under the Boost Software License, Version 1.0. 10 // (See accompanying file LICENSE_1_0.txt or copy at 11 // http://www.boost.org/LICENSE_1_0.txt) 12 13 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_MULTIPOINT_HPP 14 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_MULTIPOINT_HPP 15 16 #include <cstddef> 17 #include <algorithm> 18 #include <utility> 19 #include <vector> 20 21 #include <boost/algorithm/minmax_element.hpp> 22 #include <boost/range.hpp> 23 24 #include <boost/geometry/core/access.hpp> 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/tags.hpp> 29 30 #include <boost/geometry/util/math.hpp> 31 #include <boost/geometry/util/range.hpp> 32 33 #include <boost/geometry/geometries/helper_geometry.hpp> 34 35 #include <boost/geometry/algorithms/detail/envelope/box.hpp> 36 #include <boost/geometry/algorithms/detail/envelope/initialize.hpp> 37 #include <boost/geometry/algorithms/detail/envelope/range.hpp> 38 #include <boost/geometry/algorithms/detail/expand/point.hpp> 39 40 #include <boost/geometry/strategies/cartesian/envelope_point.hpp> 41 #include <boost/geometry/strategies/normalize.hpp> 42 #include <boost/geometry/strategies/spherical/envelope_box.hpp> 43 #include <boost/geometry/strategies/spherical/envelope_point.hpp> 44 45 46 namespace boost { namespace geometry 47 { 48 49 namespace strategy { namespace envelope 50 { 51 52 class spherical_multipoint 53 { 54 private: 55 template <std::size_t Dim> 56 struct coordinate_less 57 { 58 template <typename Point> operator ()boost::geometry::strategy::envelope::spherical_multipoint::coordinate_less59 inline bool operator()(Point const& point1, Point const& point2) const 60 { 61 return math::smaller(geometry::get<Dim>(point1), 62 geometry::get<Dim>(point2)); 63 } 64 }; 65 66 template <typename Constants, typename MultiPoint, typename OutputIterator> analyze_point_coordinates(MultiPoint const & multipoint,bool & has_south_pole,bool & has_north_pole,OutputIterator oit)67 static inline void analyze_point_coordinates(MultiPoint const& multipoint, 68 bool& has_south_pole, 69 bool& has_north_pole, 70 OutputIterator oit) 71 { 72 typedef typename boost::range_value<MultiPoint>::type point_type; 73 typedef typename boost::range_iterator 74 < 75 MultiPoint const 76 >::type iterator_type; 77 78 // analyze point coordinates: 79 // (1) normalize point coordinates 80 // (2) check if any point is the north or the south pole 81 // (3) put all non-pole points in a container 82 // 83 // notice that at this point in the algorithm, we have at 84 // least two points on the spheroid 85 has_south_pole = false; 86 has_north_pole = false; 87 88 for (iterator_type it = boost::begin(multipoint); 89 it != boost::end(multipoint); 90 ++it) 91 { 92 point_type point; 93 normalize::spherical_point::apply(*it, point); 94 95 if (math::equals(geometry::get<1>(point), 96 Constants::min_latitude())) 97 { 98 has_south_pole = true; 99 } 100 else if (math::equals(geometry::get<1>(point), 101 Constants::max_latitude())) 102 { 103 has_north_pole = true; 104 } 105 else 106 { 107 *oit++ = point; 108 } 109 } 110 } 111 112 template <typename SortedRange, typename Value> maximum_gap(SortedRange const & sorted_range,Value & max_gap_left,Value & max_gap_right)113 static inline Value maximum_gap(SortedRange const& sorted_range, 114 Value& max_gap_left, 115 Value& max_gap_right) 116 { 117 typedef typename boost::range_iterator 118 < 119 SortedRange const 120 >::type iterator_type; 121 122 iterator_type it1 = boost::begin(sorted_range), it2 = it1; 123 ++it2; 124 max_gap_left = geometry::get<0>(*it1); 125 max_gap_right = geometry::get<0>(*it2); 126 127 Value max_gap = max_gap_right - max_gap_left; 128 for (++it1, ++it2; it2 != boost::end(sorted_range); ++it1, ++it2) 129 { 130 Value gap = geometry::get<0>(*it2) - geometry::get<0>(*it1); 131 if (math::larger(gap, max_gap)) 132 { 133 max_gap_left = geometry::get<0>(*it1); 134 max_gap_right = geometry::get<0>(*it2); 135 max_gap = gap; 136 } 137 } 138 139 return max_gap; 140 } 141 142 template 143 < 144 typename Constants, 145 typename PointRange, 146 typename LongitudeLess, 147 typename CoordinateType 148 > get_min_max_longitudes(PointRange & range,LongitudeLess const & lon_less,CoordinateType & lon_min,CoordinateType & lon_max)149 static inline void get_min_max_longitudes(PointRange& range, 150 LongitudeLess const& lon_less, 151 CoordinateType& lon_min, 152 CoordinateType& lon_max) 153 { 154 typedef typename boost::range_iterator 155 < 156 PointRange const 157 >::type iterator_type; 158 159 // compute min and max longitude values 160 std::pair<iterator_type, iterator_type> min_max_longitudes 161 = boost::minmax_element(boost::begin(range), 162 boost::end(range), 163 lon_less); 164 165 lon_min = geometry::get<0>(*min_max_longitudes.first); 166 lon_max = geometry::get<0>(*min_max_longitudes.second); 167 168 // if the longitude span is "large" compute the true maximum gap 169 if (math::larger(lon_max - lon_min, Constants::half_period())) 170 { 171 std::sort(boost::begin(range), boost::end(range), lon_less); 172 173 CoordinateType max_gap_left = 0, max_gap_right = 0; 174 CoordinateType max_gap 175 = maximum_gap(range, max_gap_left, max_gap_right); 176 177 CoordinateType complement_gap 178 = Constants::period() + lon_min - lon_max; 179 180 if (math::larger(max_gap, complement_gap)) 181 { 182 lon_min = max_gap_right; 183 lon_max = max_gap_left + Constants::period(); 184 } 185 } 186 } 187 188 template 189 < 190 typename Constants, 191 typename Iterator, 192 typename LatitudeLess, 193 typename CoordinateType 194 > get_min_max_latitudes(Iterator const first,Iterator const last,LatitudeLess const & lat_less,bool has_south_pole,bool has_north_pole,CoordinateType & lat_min,CoordinateType & lat_max)195 static inline void get_min_max_latitudes(Iterator const first, 196 Iterator const last, 197 LatitudeLess const& lat_less, 198 bool has_south_pole, 199 bool has_north_pole, 200 CoordinateType& lat_min, 201 CoordinateType& lat_max) 202 { 203 if (has_south_pole && has_north_pole) 204 { 205 lat_min = Constants::min_latitude(); 206 lat_max = Constants::max_latitude(); 207 } 208 else if (has_south_pole) 209 { 210 lat_min = Constants::min_latitude(); 211 lat_max 212 = geometry::get<1>(*std::max_element(first, last, lat_less)); 213 } 214 else if (has_north_pole) 215 { 216 lat_min 217 = geometry::get<1>(*std::min_element(first, last, lat_less)); 218 lat_max = Constants::max_latitude(); 219 } 220 else 221 { 222 std::pair<Iterator, Iterator> min_max_latitudes 223 = boost::minmax_element(first, last, lat_less); 224 225 lat_min = geometry::get<1>(*min_max_latitudes.first); 226 lat_max = geometry::get<1>(*min_max_latitudes.second); 227 } 228 } 229 230 public: 231 template <typename MultiPoint, typename Box> apply(MultiPoint const & multipoint,Box & mbr)232 static inline void apply(MultiPoint const& multipoint, Box& mbr) 233 { 234 typedef typename point_type<MultiPoint>::type point_type; 235 typedef typename coordinate_type<MultiPoint>::type coordinate_type; 236 typedef typename boost::range_iterator 237 < 238 MultiPoint const 239 >::type iterator_type; 240 241 typedef math::detail::constants_on_spheroid 242 < 243 coordinate_type, 244 typename geometry::detail::cs_angular_units<MultiPoint>::type 245 > constants; 246 247 if (boost::empty(multipoint)) 248 { 249 geometry::detail::envelope::initialize<Box, 0, dimension<Box>::value>::apply(mbr); 250 return; 251 } 252 253 geometry::detail::envelope::initialize<Box, 0, 2>::apply(mbr); 254 255 if (boost::size(multipoint) == 1) 256 { 257 return dispatch::envelope 258 < 259 typename boost::range_value<MultiPoint>::type 260 >::apply(range::front(multipoint), mbr, strategy::envelope::spherical_point()); 261 } 262 263 // analyze the points and put the non-pole ones in the 264 // points vector 265 std::vector<point_type> points; 266 bool has_north_pole = false, has_south_pole = false; 267 268 analyze_point_coordinates<constants>(multipoint, 269 has_south_pole, has_north_pole, 270 std::back_inserter(points)); 271 272 coordinate_type lon_min, lat_min, lon_max, lat_max; 273 if (points.size() == 1) 274 { 275 // we have one non-pole point and at least one pole point 276 lon_min = geometry::get<0>(range::front(points)); 277 lon_max = geometry::get<0>(range::front(points)); 278 lat_min = has_south_pole 279 ? constants::min_latitude() 280 : constants::max_latitude(); 281 lat_max = has_north_pole 282 ? constants::max_latitude() 283 : constants::min_latitude(); 284 } 285 else if (points.empty()) 286 { 287 // all points are pole points 288 BOOST_GEOMETRY_ASSERT(has_south_pole || has_north_pole); 289 lon_min = coordinate_type(0); 290 lon_max = coordinate_type(0); 291 lat_min = has_south_pole 292 ? constants::min_latitude() 293 : constants::max_latitude(); 294 lat_max = (has_north_pole) 295 ? constants::max_latitude() 296 : constants::min_latitude(); 297 } 298 else 299 { 300 get_min_max_longitudes<constants>(points, 301 coordinate_less<0>(), 302 lon_min, 303 lon_max); 304 305 get_min_max_latitudes<constants>(points.begin(), 306 points.end(), 307 coordinate_less<1>(), 308 has_south_pole, 309 has_north_pole, 310 lat_min, 311 lat_max); 312 } 313 314 typedef typename helper_geometry 315 < 316 Box, 317 coordinate_type, 318 typename geometry::detail::cs_angular_units<MultiPoint>::type 319 >::type helper_box_type; 320 321 helper_box_type helper_mbr; 322 323 geometry::set<min_corner, 0>(helper_mbr, lon_min); 324 geometry::set<min_corner, 1>(helper_mbr, lat_min); 325 geometry::set<max_corner, 0>(helper_mbr, lon_max); 326 geometry::set<max_corner, 1>(helper_mbr, lat_max); 327 328 // now transform to output MBR (per index) 329 geometry::detail::envelope::envelope_indexed_box_on_spheroid<min_corner, 2>::apply(helper_mbr, mbr); 330 geometry::detail::envelope::envelope_indexed_box_on_spheroid<max_corner, 2>::apply(helper_mbr, mbr); 331 332 // compute envelope for higher coordinates 333 iterator_type it = boost::begin(multipoint); 334 geometry::detail::envelope::envelope_one_point<2, dimension<Box>::value>::apply(*it, mbr); 335 336 for (++it; it != boost::end(multipoint); ++it) 337 { 338 strategy::expand::detail::point_loop 339 < 340 2, dimension<Box>::value 341 >::apply(mbr, *it); 342 } 343 } 344 }; 345 346 347 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 348 349 namespace services 350 { 351 352 template <typename CalculationType> 353 struct default_strategy<multi_point_tag, spherical_equatorial_tag, CalculationType> 354 { 355 typedef strategy::envelope::spherical_multipoint type; 356 }; 357 358 template <typename CalculationType> 359 struct default_strategy<multi_point_tag, spherical_polar_tag, CalculationType> 360 { 361 typedef strategy::envelope::spherical_multipoint type; 362 }; 363 364 template <typename CalculationType> 365 struct default_strategy<multi_point_tag, geographic_tag, CalculationType> 366 { 367 typedef strategy::envelope::spherical_multipoint type; 368 }; 369 370 } // namespace services 371 372 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS 373 374 375 }} // namespace strategy::envelope 376 377 }} // namespace boost::geometry 378 379 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_MULTIPOINT_HPP 380