segment.hpp 16 KB

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