]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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, 2016. | |
8 | // Modifications copyright (c) 2015-2016, Oracle and/or its affiliates. | |
9 | ||
10 | // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle | |
11 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
12 | ||
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) | |
16 | ||
17 | #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP | |
18 | #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP | |
19 | ||
20 | #include <cstddef> | |
21 | #include <utility> | |
22 | ||
23 | #include <boost/numeric/conversion/cast.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/util/math.hpp> | |
34 | ||
35 | #include <boost/geometry/geometries/helper_geometry.hpp> | |
36 | ||
37 | #include <boost/geometry/strategies/compare.hpp> | |
38 | ||
39 | #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> | |
40 | #include <boost/geometry/algorithms/detail/normalize.hpp> | |
41 | ||
42 | #include <boost/geometry/algorithms/detail/envelope/point.hpp> | |
43 | #include <boost/geometry/algorithms/detail/envelope/transform_units.hpp> | |
44 | ||
45 | #include <boost/geometry/algorithms/detail/expand/point.hpp> | |
46 | ||
47 | #include <boost/geometry/algorithms/dispatch/envelope.hpp> | |
48 | ||
49 | ||
50 | namespace boost { namespace geometry | |
51 | { | |
52 | ||
53 | #ifndef DOXYGEN_NO_DETAIL | |
54 | namespace detail { namespace envelope | |
55 | { | |
56 | ||
57 | ||
58 | template <std::size_t Dimension, std::size_t DimensionCount> | |
59 | struct envelope_one_segment | |
60 | { | |
61 | template<typename Point, typename Box> | |
62 | static inline void apply(Point const& p1, Point const& p2, Box& mbr) | |
63 | { | |
64 | envelope_one_point<Dimension, DimensionCount>::apply(p1, mbr); | |
65 | detail::expand::point_loop | |
66 | < | |
67 | strategy::compare::default_strategy, | |
68 | strategy::compare::default_strategy, | |
69 | Dimension, | |
70 | DimensionCount | |
71 | >::apply(mbr, p2); | |
72 | } | |
73 | }; | |
74 | ||
75 | ||
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 | |
78 | // segment. | |
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 | |
82 | // those of the box. | |
83 | class compute_mbr_of_segment | |
84 | { | |
85 | private: | |
86 | // computes the azimuths of the segment with endpoints (lon1, lat1) | |
87 | // and (lon2, lat2) | |
88 | // radians | |
89 | template <typename CalculationType> | |
90 | static inline void azimuths(CalculationType const& lon1, | |
91 | CalculationType const& lat1, | |
92 | CalculationType const& lon2, | |
93 | CalculationType const& lat2, | |
94 | CalculationType& a1, | |
95 | CalculationType& a2) | |
96 | { | |
97 | BOOST_GEOMETRY_ASSERT(lon1 <= lon2); | |
98 | ||
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); | |
106 | ||
107 | a1 = atan2(sin_dlon * cos_lat2, | |
108 | cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon); | |
109 | ||
110 | a2 = atan2(-sin_dlon * cos_lat1, | |
111 | cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon); | |
112 | a2 += math::pi<CalculationType>(); | |
113 | } | |
114 | ||
115 | // degrees or radians | |
116 | template <typename CalculationType> | |
117 | static inline void swap(CalculationType& lon1, | |
118 | CalculationType& lat1, | |
119 | CalculationType& lon2, | |
120 | CalculationType& lat2) | |
121 | { | |
122 | std::swap(lon1, lon2); | |
123 | std::swap(lat1, lat2); | |
124 | } | |
125 | ||
126 | // radians | |
127 | template <typename CalculationType> | |
128 | static inline bool contains_pi_half(CalculationType const& a1, | |
129 | CalculationType const& a2) | |
130 | { | |
131 | // azimuths a1 and a2 are assumed to be in radians | |
132 | BOOST_GEOMETRY_ASSERT(! math::equals(a1, a2)); | |
133 | ||
134 | static CalculationType const pi_half = math::half_pi<CalculationType>(); | |
135 | ||
136 | return (a1 < a2) | |
137 | ? (a1 < pi_half && pi_half < a2) | |
138 | : (a1 > pi_half && pi_half > a2); | |
139 | } | |
140 | ||
141 | // radians or degrees | |
142 | template <typename Units, typename CoordinateType> | |
143 | static inline bool crosses_antimeridian(CoordinateType const& lon1, | |
144 | CoordinateType const& lon2) | |
145 | { | |
146 | typedef math::detail::constants_on_spheroid | |
147 | < | |
148 | CoordinateType, Units | |
149 | > constants; | |
150 | ||
151 | return math::abs(lon1 - lon2) > constants::half_period(); // > pi | |
152 | } | |
153 | ||
154 | // radians | |
155 | template <typename CalculationType> | |
156 | static inline CalculationType max_latitude(CalculationType const& azimuth, | |
157 | CalculationType const& latitude) | |
158 | { | |
159 | // azimuth and latitude are assumed to be in radians | |
160 | return acos( math::abs(cos(latitude) * sin(azimuth)) ); | |
161 | } | |
162 | ||
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) | |
169 | { | |
170 | // coordinates are assumed to be in radians | |
171 | BOOST_GEOMETRY_ASSERT(lon1 <= lon2); | |
172 | ||
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); | |
177 | ||
178 | CalculationType a1 = 0, a2 = 0; | |
179 | azimuths(lon1_rad, lat1_rad, lon2_rad, lat2_rad, a1, a2); | |
180 | ||
181 | if (lat1 > lat2) | |
182 | { | |
183 | std::swap(lat1, lat2); | |
184 | std::swap(lat1_rad, lat2_rad); | |
185 | } | |
186 | ||
187 | if (math::equals(a1, a2)) | |
188 | { | |
189 | // the segment must lie on the equator or is very short | |
190 | return; | |
191 | } | |
192 | ||
193 | if (contains_pi_half(a1, a2)) | |
194 | { | |
195 | CalculationType const mid_lat = lat1 + lat2; | |
196 | if (mid_lat < 0) | |
197 | { | |
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); | |
201 | ||
202 | if (lat1 > lat_min) | |
203 | { | |
204 | lat1 = lat_min; | |
205 | } | |
206 | } | |
207 | else if (mid_lat > 0) | |
208 | { | |
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); | |
212 | ||
213 | if (lat2 < lat_max) | |
214 | { | |
215 | lat2 = lat_max; | |
216 | } | |
217 | } | |
218 | } | |
219 | } | |
220 | ||
221 | template <typename Units, typename CalculationType> | |
222 | static inline void apply(CalculationType& lon1, | |
223 | CalculationType& lat1, | |
224 | CalculationType& lon2, | |
225 | CalculationType& lat2) | |
226 | { | |
227 | typedef math::detail::constants_on_spheroid | |
228 | < | |
229 | CalculationType, Units | |
230 | > constants; | |
231 | ||
232 | bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude()); | |
233 | bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude()); | |
234 | ||
235 | if (is_pole1 && is_pole2) | |
236 | { | |
237 | // both points are poles; nothing more to do: | |
238 | // longitudes are already normalized to 0 | |
239 | // but just in case | |
240 | lon1 = 0; | |
241 | lon2 = 0; | |
242 | } | |
243 | else if (is_pole1 && !is_pole2) | |
244 | { | |
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 | |
248 | lon1 = lon2; | |
249 | } | |
250 | else if (!is_pole1 && is_pole2) | |
251 | { | |
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 | |
255 | lon2 = lon1; | |
256 | } | |
257 | ||
258 | if (lon1 == lon2) | |
259 | { | |
260 | // segment lies on a meridian | |
261 | if (lat1 > lat2) | |
262 | { | |
263 | std::swap(lat1, lat2); | |
264 | } | |
265 | return; | |
266 | } | |
267 | ||
268 | BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2); | |
269 | ||
270 | if (lon1 > lon2) | |
271 | { | |
272 | swap(lon1, lat1, lon2, lat2); | |
273 | } | |
274 | ||
275 | if (crosses_antimeridian<Units>(lon1, lon2)) | |
276 | { | |
277 | lon1 += constants::period(); | |
278 | swap(lon1, lat1, lon2, lat2); | |
279 | } | |
280 | ||
281 | compute_box_corners<Units>(lon1, lat1, lon2, lat2); | |
282 | } | |
283 | ||
284 | public: | |
285 | template <typename Units, typename CalculationType, typename Box> | |
286 | static inline void apply(CalculationType lon1, | |
287 | CalculationType lat1, | |
288 | CalculationType lon2, | |
289 | CalculationType lat2, | |
290 | Box& mbr) | |
291 | { | |
292 | typedef typename coordinate_type<Box>::type box_coordinate_type; | |
293 | ||
294 | typedef typename helper_geometry | |
295 | < | |
296 | Box, box_coordinate_type, Units | |
297 | >::type helper_box_type; | |
298 | ||
299 | helper_box_type radian_mbr; | |
300 | ||
301 | apply<Units>(lon1, lat1, lon2, lat2); | |
302 | ||
303 | geometry::set | |
304 | < | |
305 | min_corner, 0 | |
306 | >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon1)); | |
307 | ||
308 | geometry::set | |
309 | < | |
310 | min_corner, 1 | |
311 | >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat1)); | |
312 | ||
313 | geometry::set | |
314 | < | |
315 | max_corner, 0 | |
316 | >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lon2)); | |
317 | ||
318 | geometry::set | |
319 | < | |
320 | max_corner, 1 | |
321 | >(radian_mbr, boost::numeric_cast<box_coordinate_type>(lat2)); | |
322 | ||
323 | transform_units(radian_mbr, mbr); | |
324 | } | |
325 | }; | |
326 | ||
327 | ||
328 | template <std::size_t DimensionCount> | |
329 | struct envelope_segment_on_sphere | |
330 | { | |
331 | template <typename Point, typename Box> | |
332 | static inline void apply(Point const& p1, Point const& p2, Box& mbr) | |
333 | { | |
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); | |
337 | ||
338 | typedef typename coordinate_system<Point>::type::units units_type; | |
339 | ||
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), | |
345 | mbr); | |
346 | ||
347 | // now compute the envelope range for coordinates of | |
348 | // dimension 2 and higher | |
349 | envelope_one_segment<2, DimensionCount>::apply(p1, p2, mbr); | |
350 | } | |
351 | ||
352 | template <typename Segment, typename Box> | |
353 | static inline void apply(Segment const& segment, Box& mbr) | |
354 | { | |
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); | |
359 | } | |
360 | }; | |
361 | ||
362 | ||
363 | ||
364 | template <std::size_t DimensionCount, typename CS_Tag> | |
365 | struct envelope_segment | |
366 | : envelope_one_segment<0, DimensionCount> | |
367 | {}; | |
368 | ||
369 | ||
370 | template <std::size_t DimensionCount> | |
371 | struct envelope_segment<DimensionCount, spherical_equatorial_tag> | |
372 | : envelope_segment_on_sphere<DimensionCount> | |
373 | {}; | |
374 | ||
375 | ||
376 | ||
377 | }} // namespace detail::envelope | |
378 | #endif // DOXYGEN_NO_DETAIL | |
379 | ||
380 | ||
381 | #ifndef DOXYGEN_NO_DISPATCH | |
382 | namespace dispatch | |
383 | { | |
384 | ||
385 | ||
386 | template <typename Segment, typename CS_Tag> | |
387 | struct envelope<Segment, segment_tag, CS_Tag> | |
388 | { | |
389 | template <typename Box> | |
390 | static inline void apply(Segment const& segment, Box& mbr) | |
391 | { | |
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 | |
396 | < | |
397 | dimension<Segment>::value, CS_Tag | |
398 | >::apply(p[0], p[1], mbr); | |
399 | } | |
400 | }; | |
401 | ||
402 | ||
403 | } // namespace dispatch | |
404 | #endif // DOXYGEN_NO_DISPATCH | |
405 | ||
406 | }} // namespace boost::geometry | |
407 | ||
408 | #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ENVELOPE_SEGMENT_HPP |