1 // Boost.Geometry (aka GGL, Generic Geometry Library)
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.
7 // This file was modified by Oracle on 2015, 2016.
8 // Modifications copyright (c) 2015-2016, Oracle and/or its affiliates.
10 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
13 // Distributed under the Boost Software License, Version 1.0.
14 // (See accompanying file LICENSE_1_0.txt or copy at
15 // http://www.boost.org/LICENSE_1_0.txt)
17 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
18 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP
23 #include <boost/numeric/conversion/cast.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/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>
33 #include <boost/geometry/util/math.hpp>
35 #include <boost/geometry/geometries/helper_geometry.hpp>
37 #include <boost/geometry/strategies/compare.hpp>
39 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
40 #include <boost/geometry/algorithms/detail/normalize.hpp>
42 #include <boost/geometry/algorithms/detail/envelope/point.hpp>
43 #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp>
45 #include <boost/geometry/algorithms/detail/expand/point.hpp>
47 #include <boost/geometry/algorithms/dispatch/envelope.hpp>
50 namespace boost { namespace geometry
53 #ifndef DOXYGEN_NO_DETAIL
54 namespace detail { namespace envelope
58 template <std::size_t Dimension, std::size_t DimensionCount>
59 struct envelope_one_segment
61 template<typename Point, typename Box>
62 static inline void apply(Point const& p1, Point const& p2, Box& mbr)
64 envelope_one_point<Dimension, DimensionCount>::apply(p1, mbr);
65 detail::expand::point_loop
67 strategy::compare::default_strategy,
68 strategy::compare::default_strategy,
76 // Computes the MBR of a segment given by (lon1, lat1) and (lon2,
77 // lat2), and with azimuths a1 and a2 at the two endpoints of the
79 // It is assumed that the spherical coordinates of the segment are
80 // normalized and in radians.
81 // The longitudes and latitudes of the endpoints are overridden by
83 class compute_mbr_of_segment
86 // computes the azimuths of the segment with endpoints (lon1, lat1)
89 template <typename CalculationType>
90 static inline void azimuths(CalculationType const& lon1,
91 CalculationType const& lat1,
92 CalculationType const& lon2,
93 CalculationType const& lat2,
97 BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
99 CalculationType dlon = lon2 - lon1;
100 CalculationType sin_dlon = sin(dlon);
101 CalculationType cos_dlon = cos(dlon);
102 CalculationType cos_lat1 = cos(lat1);
103 CalculationType cos_lat2 = cos(lat2);
104 CalculationType sin_lat1 = sin(lat1);
105 CalculationType sin_lat2 = sin(lat2);
107 a1 = atan2(sin_dlon * cos_lat2,
108 cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon);
110 a2 = atan2(-sin_dlon * cos_lat1,
111 cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon);
112 a2 += math::pi<CalculationType>();
115 // degrees or radians
116 template <typename CalculationType>
117 static inline void swap(CalculationType& lon1,
118 CalculationType& lat1,
119 CalculationType& lon2,
120 CalculationType& lat2)
122 std::swap(lon1, lon2);
123 std::swap(lat1, lat2);
127 template <typename CalculationType>
128 static inline bool contains_pi_half(CalculationType const& a1,
129 CalculationType const& a2)
131 // azimuths a1 and a2 are assumed to be in radians
132 BOOST_GEOMETRY_ASSERT(! math::equals(a1, a2));
134 static CalculationType const pi_half = math::half_pi<CalculationType>();
137 ? (a1 < pi_half && pi_half < a2)
138 : (a1 > pi_half && pi_half > a2);
141 // radians or degrees
142 template <typename Units, typename CoordinateType>
143 static inline bool crosses_antimeridian(CoordinateType const& lon1,
144 CoordinateType const& lon2)
146 typedef math::detail::constants_on_spheroid
148 CoordinateType, Units
151 return math::abs(lon1 - lon2) > constants::half_period(); // > pi
155 template <typename CalculationType>
156 static inline CalculationType max_latitude(CalculationType const& azimuth,
157 CalculationType const& latitude)
159 // azimuth and latitude are assumed to be in radians
160 return acos( math::abs(cos(latitude) * sin(azimuth)) );
163 // degrees or radians
164 template <typename Units, typename CalculationType>
165 static inline void compute_box_corners(CalculationType& lon1,
166 CalculationType& lat1,
167 CalculationType& lon2,
168 CalculationType& lat2)
170 // coordinates are assumed to be in radians
171 BOOST_GEOMETRY_ASSERT(lon1 <= lon2);
173 CalculationType lon1_rad = math::as_radian<Units>(lon1);
174 CalculationType lat1_rad = math::as_radian<Units>(lat1);
175 CalculationType lon2_rad = math::as_radian<Units>(lon2);
176 CalculationType lat2_rad = math::as_radian<Units>(lat2);
178 CalculationType a1 = 0, a2 = 0;
179 azimuths(lon1_rad, lat1_rad, lon2_rad, lat2_rad, a1, a2);
183 std::swap(lat1, lat2);
184 std::swap(lat1_rad, lat2_rad);
187 if (math::equals(a1, a2))
189 // the segment must lie on the equator or is very short
193 if (contains_pi_half(a1, a2))
195 CalculationType const mid_lat = lat1 + lat2;
198 // update using min latitude
199 CalculationType const lat_min_rad = -max_latitude(a1, lat1_rad);
200 CalculationType const lat_min = math::from_radian<Units>(lat_min_rad);
207 else if (mid_lat > 0)
209 // update using max latitude
210 CalculationType const lat_max_rad = max_latitude(a1, lat1_rad);
211 CalculationType const lat_max = math::from_radian<Units>(lat_max_rad);
221 template <typename Units, typename CalculationType>
222 static inline void apply(CalculationType& lon1,
223 CalculationType& lat1,
224 CalculationType& lon2,
225 CalculationType& lat2)
227 typedef math::detail::constants_on_spheroid
229 CalculationType, Units
232 bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude());
233 bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude());
235 if (is_pole1 && is_pole2)
237 // both points are poles; nothing more to do:
238 // longitudes are already normalized to 0
243 else if (is_pole1 && !is_pole2)
245 // first point is a pole, second point is not:
246 // make the longitude of the first point the same as that
247 // of the second point
250 else if (!is_pole1 && is_pole2)
252 // second point is a pole, first point is not:
253 // make the longitude of the second point the same as that
254 // of the first point
260 // segment lies on a meridian
263 std::swap(lat1, lat2);
268 BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2);
272 swap(lon1, lat1, lon2, lat2);
275 if (crosses_antimeridian<Units>(lon1, lon2))
277 lon1 += constants::period();
278 swap(lon1, lat1, lon2, lat2);
281 compute_box_corners<Units>(lon1, lat1, lon2, lat2);
285 template <typename Units, typename CalculationType, typename Box>
286 static inline void apply(CalculationType lon1,
287 CalculationType lat1,
288 CalculationType lon2,
289 CalculationType lat2,
292 typedef typename coordinate_type<Box>::type box_coordinate_type;
294 typedef typename helper_geometry
296 Box, box_coordinate_type, Units
297 >::type helper_box_type;
299 helper_box_type radian_mbr;
301 apply<Units>(lon1, lat1, lon2, lat2);
306 >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon1));
311 >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat1));
316 >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon2));
321 >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat2));
323 transform_units(radian_mbr, mbr);
328 template <std::size_t DimensionCount>
329 struct envelope_segment_on_sphere
331 template <typename Point, typename Box>
332 static inline void apply(Point const& p1, Point const& p2, Box& mbr)
334 // first compute the envelope range for the first two coordinates
335 Point p1_normalized = detail::return_normalized<Point>(p1);
336 Point p2_normalized = detail::return_normalized<Point>(p2);
338 typedef typename coordinate_system<Point>::type::units units_type;
340 compute_mbr_of_segment::template apply<units_type>(
341 geometry::get<0>(p1_normalized),
342 geometry::get<1>(p1_normalized),
343 geometry::get<0>(p2_normalized),
344 geometry::get<1>(p2_normalized),
347 // now compute the envelope range for coordinates of
348 // dimension 2 and higher
349 envelope_one_segment<2, DimensionCount>::apply(p1, p2, mbr);
352 template <typename Segment, typename Box>
353 static inline void apply(Segment const& segment, Box& mbr)
355 typename point_type<Segment>::type p[2];
356 detail::assign_point_from_index<0>(segment, p[0]);
357 detail::assign_point_from_index<1>(segment, p[1]);
358 apply(p[0], p[1], mbr);
364 template <std::size_t DimensionCount, typename CS_Tag>
365 struct envelope_segment
366 : envelope_one_segment<0, DimensionCount>
370 template <std::size_t DimensionCount>
371 struct envelope_segment<DimensionCount, spherical_equatorial_tag>
372 : envelope_segment_on_sphere<DimensionCount>
377 }} // namespace detail::envelope
378 #endif // DOXYGEN_NO_DETAIL
381 #ifndef DOXYGEN_NO_DISPATCH
386 template <typename Segment, typename CS_Tag>
387 struct envelope<Segment, segment_tag, CS_Tag>
389 template <typename Box>
390 static inline void apply(Segment const& segment, Box& mbr)
392 typename point_type<Segment>::type p[2];
393 detail::assign_point_from_index<0>(segment, p[0]);
394 detail::assign_point_from_index<1>(segment, p[1]);
395 detail::envelope::envelope_segment
397 dimension<Segment>::value, CS_Tag
398 >::apply(p[0], p[1], mbr);
403 } // namespace dispatch
404 #endif // DOXYGEN_NO_DISPATCH
406 }} // namespace boost::geometry
408 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP