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