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