1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 
3 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
4 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
5 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
6 
7 // This file was modified by Oracle on 2015-2017.
8 // Modifications copyright (c) 2015-2017, Oracle and/or its affiliates.
9 
10 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
11 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
12 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
13 
14 // Distributed under the Boost Software License, Version 1.0.
15 // (See accompanying file LICENSE_1_0.txt or copy at
16 // http://www.boost.org/LICENSE_1_0.txt)
17 
18 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
19 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
20 
21 #include <cstddef>
22 #include <utility>
23 
24 #include <boost/numeric/conversion/cast.hpp>
25 
26 #include <boost/geometry/core/assert.hpp>
27 #include <boost/geometry/core/coordinate_system.hpp>
28 #include <boost/geometry/core/coordinate_type.hpp>
29 #include <boost/geometry/core/cs.hpp>
30 #include <boost/geometry/core/srs.hpp>
31 #include <boost/geometry/core/point_type.hpp>
32 #include <boost/geometry/core/radian_access.hpp>
33 #include <boost/geometry/core/tags.hpp>
34 
35 #include <boost/geometry/util/math.hpp>
36 
37 #include <boost/geometry/geometries/helper_geometry.hpp>
38 
39 #include <boost/geometry/formulas/vertex_latitude.hpp>
40 
41 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
42 
43 #include <boost/geometry/algorithms/detail/envelope/point.hpp>
44 #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp>
45 
46 #include <boost/geometry/algorithms/detail/expand/point.hpp>
47 
48 #include <boost/geometry/algorithms/dispatch/envelope.hpp>
49 
50 namespace boost { namespace geometry
51 {
52 
53 #ifndef DOXYGEN_NO_DETAIL
54 namespace detail { namespace envelope
55 {
56 
57 template <typename CalculationType, typename CS_Tag>
58 struct envelope_segment_call_vertex_latitude
59 {
60     template <typename T1, typename T2, typename Strategy>
applyboost::geometry::detail::envelope::envelope_segment_call_vertex_latitude61     static inline CalculationType apply(T1 const& lat1,
62                                         T2 const& alp1,
63                                         Strategy const& )
64     {
65         return geometry::formula::vertex_latitude<CalculationType, CS_Tag>
66             ::apply(lat1, alp1);
67     }
68 };
69 
70 template <typename CalculationType>
71 struct envelope_segment_call_vertex_latitude<CalculationType, geographic_tag>
72 {
73     template <typename T1, typename T2, typename Strategy>
applyboost::geometry::detail::envelope::envelope_segment_call_vertex_latitude74     static inline CalculationType apply(T1 const& lat1,
75                                         T2 const& alp1,
76                                         Strategy const& strategy)
77     {
78         return geometry::formula::vertex_latitude<CalculationType, geographic_tag>
79             ::apply(lat1, alp1, strategy.model());
80     }
81 };
82 
83 template <typename CS_Tag>
84 class envelope_segment_impl
85 {
86 private:
87 
88     // degrees or radians
89     template <typename CalculationType>
swap(CalculationType & lon1,CalculationType & lat1,CalculationType & lon2,CalculationType & lat2)90     static inline void swap(CalculationType& lon1,
91                             CalculationType& lat1,
92                             CalculationType& lon2,
93                             CalculationType& lat2)
94     {
95         std::swap(lon1, lon2);
96         std::swap(lat1, lat2);
97     }
98 
99     // radians
100     template <typename CalculationType>
contains_pi_half(CalculationType const & a1,CalculationType const & a2)101     static inline bool contains_pi_half(CalculationType const& a1,
102                                         CalculationType const& a2)
103     {
104         // azimuths a1 and a2 are assumed to be in radians
105         BOOST_GEOMETRY_ASSERT(! math::equals(a1, a2));
106 
107         static CalculationType const pi_half = math::half_pi<CalculationType>();
108 
109         return (a1 < a2)
110                 ? (a1 < pi_half && pi_half < a2)
111                 : (a1 > pi_half && pi_half > a2);
112     }
113 
114     // radians or degrees
115     template <typename Units, typename CoordinateType>
crosses_antimeridian(CoordinateType const & lon1,CoordinateType const & lon2)116     static inline bool crosses_antimeridian(CoordinateType const& lon1,
117                                             CoordinateType const& lon2)
118     {
119         typedef math::detail::constants_on_spheroid
120             <
121                 CoordinateType, Units
122             > constants;
123 
124         return math::abs(lon1 - lon2) > constants::half_period(); // > pi
125     }
126 
127     // degrees or radians
128     template <typename Units, typename CalculationType, typename Strategy>
compute_box_corners(CalculationType & lon1,CalculationType & lat1,CalculationType & lon2,CalculationType & lat2,CalculationType a1,CalculationType a2,Strategy const & strategy)129     static inline void compute_box_corners(CalculationType& lon1,
130                                            CalculationType& lat1,
131                                            CalculationType& lon2,
132                                            CalculationType& lat2,
133                                            CalculationType a1,
134                                            CalculationType a2,
135                                            Strategy const& strategy)
136     {
137         // coordinates are assumed to be in radians
138         BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
139 
140         CalculationType lat1_rad = math::as_radian<Units>(lat1);
141         CalculationType lat2_rad = math::as_radian<Units>(lat2);
142 
143         if (lat1 > lat2)
144         {
145             std::swap(lat1, lat2);
146             std::swap(lat1_rad, lat2_rad);
147             std::swap(a1, a2);
148         }
149 
150         if (math::equals(a1, a2))
151         {
152             // the segment must lie on the equator or is very short
153             return;
154         }
155 
156         if (contains_pi_half(a1, a2))
157         {
158             CalculationType p_max = envelope_segment_call_vertex_latitude
159                 <CalculationType, CS_Tag>::apply(lat1_rad, a1, strategy);
160 
161             CalculationType const mid_lat = lat1 + lat2;
162             if (mid_lat < 0)
163             {
164                 // update using min latitude
165                 CalculationType const lat_min_rad = -p_max;
166                 CalculationType const lat_min
167                     = math::from_radian<Units>(lat_min_rad);
168 
169                 if (lat1 > lat_min)
170                 {
171                     lat1 = lat_min;
172                 }
173             }
174             else if (mid_lat > 0)
175             {
176                 // update using max latitude
177                 CalculationType const lat_max_rad = p_max;
178                 CalculationType const lat_max
179                     = math::from_radian<Units>(lat_max_rad);
180 
181                 if (lat2 < lat_max)
182                 {
183                     lat2 = lat_max;
184                 }
185             }
186         }
187     }
188 
189     template <typename Units, typename CalculationType>
special_cases(CalculationType & lon1,CalculationType & lat1,CalculationType & lon2,CalculationType & lat2)190     static inline void special_cases(CalculationType& lon1,
191                                      CalculationType& lat1,
192                                      CalculationType& lon2,
193                                      CalculationType& lat2)
194     {
195         typedef math::detail::constants_on_spheroid
196             <
197                 CalculationType, Units
198             > constants;
199 
200         bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
201         bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());
202 
203         if (is_pole1 && is_pole2)
204         {
205             // both points are poles; nothing more to do:
206             // longitudes are already normalized to 0
207             // but just in case
208             lon1 = 0;
209             lon2 = 0;
210         }
211         else if (is_pole1 && !is_pole2)
212         {
213             // first point is a pole, second point is not:
214             // make the longitude of the first point the same as that
215             // of the second point
216             lon1 = lon2;
217         }
218         else if (!is_pole1 && is_pole2)
219         {
220             // second point is a pole, first point is not:
221             // make the longitude of the second point the same as that
222             // of the first point
223             lon2 = lon1;
224         }
225 
226         if (lon1 == lon2)
227         {
228             // segment lies on a meridian
229             if (lat1 > lat2)
230             {
231                 std::swap(lat1, lat2);
232             }
233             return;
234         }
235 
236         BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);
237 
238         if (lon1 > lon2)
239         {
240             swap(lon1, lat1, lon2, lat2);
241         }
242 
243         if (crosses_antimeridian<Units>(lon1, lon2))
244         {
245             lon1 += constants::period();
246             swap(lon1, lat1, lon2, lat2);
247         }
248     }
249 
250     template
251     <
252         typename Units,
253         typename CalculationType,
254         typename Box
255     >
create_box(CalculationType lon1,CalculationType lat1,CalculationType lon2,CalculationType lat2,Box & mbr)256     static inline void create_box(CalculationType lon1,
257                                   CalculationType lat1,
258                                   CalculationType lon2,
259                                   CalculationType lat2,
260                                   Box& mbr)
261     {
262         typedef typename coordinate_type<Box>::type box_coordinate_type;
263 
264         typedef typename helper_geometry
265             <
266                 Box, box_coordinate_type, Units
267             >::type helper_box_type;
268 
269         helper_box_type radian_mbr;
270 
271         geometry::set
272             <
273                 min_corner, 0
274             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon1));
275 
276         geometry::set
277             <
278                 min_corner, 1
279             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat1));
280 
281         geometry::set
282             <
283                 max_corner, 0
284             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon2));
285 
286         geometry::set
287             <
288                 max_corner, 1
289             >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat2));
290 
291         transform_units(radian_mbr, mbr);
292     }
293 
294 
295     template <typename Units, typename CalculationType, typename Strategy>
apply(CalculationType & lon1,CalculationType & lat1,CalculationType & lon2,CalculationType & lat2,Strategy const & strategy)296     static inline void apply(CalculationType& lon1,
297                              CalculationType& lat1,
298                              CalculationType& lon2,
299                              CalculationType& lat2,
300                              Strategy const& strategy)
301     {
302         special_cases<Units>(lon1, lat1, lon2, lat2);
303 
304         CalculationType lon1_rad = math::as_radian<Units>(lon1);
305         CalculationType lat1_rad = math::as_radian<Units>(lat1);
306         CalculationType lon2_rad = math::as_radian<Units>(lon2);
307         CalculationType lat2_rad = math::as_radian<Units>(lat2);
308         CalculationType alp1, alp2;
309         strategy.apply(lon1_rad, lat1_rad, lon2_rad, lat2_rad, alp1, alp2);
310 
311         compute_box_corners<Units>(lon1, lat1, lon2, lat2, alp1, alp2, strategy);
312     }
313 
314     template <typename Units, typename CalculationType, typename Strategy>
apply(CalculationType & lon1,CalculationType & lat1,CalculationType & lon2,CalculationType & lat2,Strategy const & strategy,CalculationType alp1)315     static inline void apply(CalculationType& lon1,
316                              CalculationType& lat1,
317                              CalculationType& lon2,
318                              CalculationType& lat2,
319                              Strategy const& strategy,
320                              CalculationType alp1)
321     {
322         special_cases<Units>(lon1, lat1, lon2, lat2);
323 
324         CalculationType lon1_rad = math::as_radian<Units>(lon1);
325         CalculationType lat1_rad = math::as_radian<Units>(lat1);
326         CalculationType lon2_rad = math::as_radian<Units>(lon2);
327         CalculationType lat2_rad = math::as_radian<Units>(lat2);
328         CalculationType alp2;
329         strategy.apply(lon2_rad, lat2_rad, lon1_rad, lat1_rad, alp2);
330         alp2 += math::pi<CalculationType>();
331 
332         compute_box_corners<Units>(lon1, lat1, lon2, lat2, alp1, alp2, strategy);
333     }
334 
335 public:
336     template
337     <
338         typename Units,
339         typename CalculationType,
340         typename Box,
341         typename Strategy
342     >
apply(CalculationType lon1,CalculationType lat1,CalculationType lon2,CalculationType lat2,Box & mbr,Strategy const & strategy)343     static inline void apply(CalculationType lon1,
344                              CalculationType lat1,
345                              CalculationType lon2,
346                              CalculationType lat2,
347                              Box& mbr,
348                              Strategy const& strategy)
349     {
350         apply<Units>(lon1, lat1, lon2, lat2, strategy);
351         create_box<Units>(lon1, lat1, lon2, lat2, mbr);
352     }
353 
354     template
355     <
356         typename Units,
357         typename CalculationType,
358         typename Box,
359         typename Strategy
360     >
apply(CalculationType lon1,CalculationType lat1,CalculationType lon2,CalculationType lat2,Box & mbr,Strategy const & strategy,CalculationType alp1)361     static inline void apply(CalculationType lon1,
362                              CalculationType lat1,
363                              CalculationType lon2,
364                              CalculationType lat2,
365                              Box& mbr,
366                              Strategy const& strategy,
367                              CalculationType alp1)
368     {
369         apply<Units>(lon1, lat1, lon2, lat2, strategy, alp1);
370         create_box<Units>(lon1, lat1, lon2, lat2, mbr);
371     }
372 };
373 
374 template <std::size_t Dimension, std::size_t DimensionCount>
375 struct envelope_one_segment
376 {
377     template<typename Point, typename Box, typename Strategy>
applyboost::geometry::detail::envelope::envelope_one_segment378     static inline void apply(Point const& p1,
379                              Point const& p2,
380                              Box& mbr,
381                              Strategy const& strategy)
382     {
383         envelope_one_point<Dimension, DimensionCount>::apply(p1, mbr, strategy);
384         detail::expand::point_loop
385             <
386                 strategy::compare::default_strategy,
387                 strategy::compare::default_strategy,
388                 Dimension,
389                 DimensionCount
390             >::apply(mbr, p2, strategy);
391     }
392 };
393 
394 
395 template <std::size_t DimensionCount>
396 struct envelope_segment
397 {
398     template <typename Point, typename Box, typename Strategy>
applyboost::geometry::detail::envelope::envelope_segment399     static inline void apply(Point const& p1,
400                              Point const& p2,
401                              Box& mbr,
402                              Strategy const& strategy)
403     {
404         // first compute the envelope range for the first two coordinates
405         strategy.apply(p1, p2, mbr);
406 
407         // now compute the envelope range for coordinates of
408         // dimension 2 and higher
409         envelope_one_segment<2, DimensionCount>::apply(p1, p2, mbr, strategy);
410     }
411 
412     template <typename Segment, typename Box>
applyboost::geometry::detail::envelope::envelope_segment413     static inline void apply(Segment const& segment, Box& mbr)
414     {
415         typename point_type<Segment>::type p[2];
416         detail::assign_point_from_index<0>(segment, p[0]);
417         detail::assign_point_from_index<1>(segment, p[1]);
418         apply(p[0], p[1], mbr);
419     }
420 };
421 
422 }} // namespace detail::envelope
423 #endif // DOXYGEN_NO_DETAIL
424 
425 
426 #ifndef DOXYGEN_NO_DISPATCH
427 namespace dispatch
428 {
429 
430 
431 template <typename Segment>
432 struct envelope<Segment, segment_tag>
433 {
434     template <typename Box, typename Strategy>
applyboost::geometry::dispatch::envelope435     static inline void apply(Segment const& segment,
436                              Box& mbr,
437                              Strategy const& strategy)
438     {
439         typename point_type<Segment>::type p[2];
440         detail::assign_point_from_index<0>(segment, p[0]);
441         detail::assign_point_from_index<1>(segment, p[1]);
442         detail::envelope::envelope_segment
443             <
444                dimension<Segment>::value
445             >::apply(p[0], p[1], mbr, strategy);
446     }
447 };
448 
449 } // namespace dispatch
450 #endif // DOXYGEN_NO_DISPATCH
451 
452 }} // namespace boost::geometry
453 
454 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
455