extreme_points.hpp 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
  3. // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
  4. // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
  5. // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
  6. // This file was modified by Oracle on 2017, 2018.
  7. // Modifications copyright (c) 2017-2018 Oracle and/or its affiliates.
  8. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  9. // Use, modification and distribution is subject to the Boost Software License,
  10. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  11. // http://www.boost.org/LICENSE_1_0.txt)
  12. #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
  13. #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
  14. #include <cstddef>
  15. #include <boost/range.hpp>
  16. #include <boost/geometry/algorithms/detail/interior_iterator.hpp>
  17. #include <boost/geometry/core/cs.hpp>
  18. #include <boost/geometry/core/point_order.hpp>
  19. #include <boost/geometry/core/point_type.hpp>
  20. #include <boost/geometry/core/ring_type.hpp>
  21. #include <boost/geometry/core/tags.hpp>
  22. #include <boost/geometry/geometries/concepts/check.hpp>
  23. #include <boost/geometry/iterators/ever_circling_iterator.hpp>
  24. #include <boost/geometry/algorithms/detail/assign_box_corners.hpp>
  25. #include <boost/geometry/strategies/side.hpp>
  26. #include <boost/geometry/util/math.hpp>
  27. namespace boost { namespace geometry
  28. {
  29. #ifndef DOXYGEN_NO_DETAIL
  30. namespace detail { namespace extreme_points
  31. {
  32. template <std::size_t Dimension>
  33. struct compare
  34. {
  35. template <typename Point>
  36. inline bool operator()(Point const& lhs, Point const& rhs)
  37. {
  38. return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs);
  39. }
  40. };
  41. template <std::size_t Dimension, typename PointType, typename CoordinateType>
  42. inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value)
  43. {
  44. // Moves a point along the vector (point, extreme) in the direction of the extreme point
  45. // This adapts the possibly uneven legs of the triangle (or trapezium-like shape)
  46. // _____extreme _____
  47. // / \ / \ .
  48. // /base \ => / \ point .
  49. // \ point .
  50. //
  51. // For so-called intruders, it can be used to adapt both legs to the level of "base"
  52. // For the base, it can be used to adapt both legs to the level of the max-value of the intruders
  53. // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector
  54. CoordinateType const value = geometry::get<Dimension>(point);
  55. //if (geometry::math::equals(value, base_value))
  56. if (value >= base_value)
  57. {
  58. return;
  59. }
  60. PointType vector = point;
  61. subtract_point(vector, extreme);
  62. CoordinateType const diff = geometry::get<Dimension>(vector);
  63. // diff should never be zero
  64. // because of the way our triangle/trapezium is build.
  65. // We just return if it would be the case.
  66. if (geometry::math::equals(diff, 0))
  67. {
  68. return;
  69. }
  70. CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme);
  71. multiply_value(vector, base_diff);
  72. divide_value(vector, diff);
  73. // The real move:
  74. point = extreme;
  75. add_point(point, vector);
  76. }
  77. template <std::size_t Dimension, typename Range, typename CoordinateType>
  78. inline void move_along_vector(Range& range, CoordinateType const& base_value)
  79. {
  80. if (range.size() >= 3)
  81. {
  82. move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value);
  83. move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value);
  84. }
  85. }
  86. template<typename Ring, std::size_t Dimension>
  87. struct extreme_points_on_ring
  88. {
  89. typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
  90. typedef typename boost::range_iterator<Ring const>::type range_iterator;
  91. typedef typename geometry::point_type<Ring>::type point_type;
  92. template <typename CirclingIterator, typename Points>
  93. static inline bool extend(CirclingIterator& it,
  94. std::size_t n,
  95. coordinate_type max_coordinate_value,
  96. Points& points, int direction)
  97. {
  98. std::size_t safe_index = 0;
  99. do
  100. {
  101. it += direction;
  102. points.push_back(*it);
  103. if (safe_index++ >= n)
  104. {
  105. // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
  106. return false;
  107. }
  108. } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
  109. return true;
  110. }
  111. // Overload without adding to poinst
  112. template <typename CirclingIterator>
  113. static inline bool extend(CirclingIterator& it,
  114. std::size_t n,
  115. coordinate_type max_coordinate_value,
  116. int direction)
  117. {
  118. std::size_t safe_index = 0;
  119. do
  120. {
  121. it += direction;
  122. if (safe_index++ >= n)
  123. {
  124. // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
  125. return false;
  126. }
  127. } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
  128. return true;
  129. }
  130. template <typename CirclingIterator>
  131. static inline bool extent_both_sides(Ring const& ring,
  132. point_type extreme,
  133. CirclingIterator& left,
  134. CirclingIterator& right)
  135. {
  136. std::size_t const n = boost::size(ring);
  137. coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
  138. if (! extend(left, n, max_coordinate_value, -1))
  139. {
  140. return false;
  141. }
  142. if (! extend(right, n, max_coordinate_value, +1))
  143. {
  144. return false;
  145. }
  146. return true;
  147. }
  148. template <typename Collection, typename CirclingIterator>
  149. static inline bool collect(Ring const& ring,
  150. point_type extreme,
  151. Collection& points,
  152. CirclingIterator& left,
  153. CirclingIterator& right)
  154. {
  155. std::size_t const n = boost::size(ring);
  156. coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
  157. // Collects first left, which is reversed (if more than one point) then adds the top itself, then right
  158. if (! extend(left, n, max_coordinate_value, points, -1))
  159. {
  160. return false;
  161. }
  162. std::reverse(points.begin(), points.end());
  163. points.push_back(extreme);
  164. if (! extend(right, n, max_coordinate_value, points, +1))
  165. {
  166. return false;
  167. }
  168. return true;
  169. }
  170. template <typename Extremes, typename Intruders, typename CirclingIterator, typename SideStrategy>
  171. static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right,
  172. Extremes const& extremes,
  173. Intruders& intruders,
  174. SideStrategy const& strategy)
  175. {
  176. if (boost::size(extremes) < 3)
  177. {
  178. return;
  179. }
  180. coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>()));
  181. // Also select left/right (if Dimension=1)
  182. coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
  183. coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
  184. std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too
  185. std::size_t const n = boost::size(ring);
  186. while (left != right && defensive_check_index < n)
  187. {
  188. coordinate_type const coordinate = geometry::get<Dimension>(*right);
  189. coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right);
  190. if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max)
  191. {
  192. int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
  193. int const first_side = strategy.apply(*right, extremes.front(), *(extremes.begin() + 1)) * factor;
  194. int const last_side = strategy.apply(*right, *(extremes.rbegin() + 1), extremes.back()) * factor;
  195. // If not lying left from any of the extemes side
  196. if (first_side != 1 && last_side != 1)
  197. {
  198. //std::cout << "first " << first_side << " last " << last_side << std::endl;
  199. // we start at this intrusion until it is handled, and don't affect our initial left iterator
  200. CirclingIterator left_intrusion_it = right;
  201. typename boost::range_value<Intruders>::type intruder;
  202. collect(ring, *right, intruder, left_intrusion_it, right);
  203. // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies
  204. // (we might postpone this action, it is often not necessary. However it is not time-consuming)
  205. move_along_vector<Dimension>(intruder, min_value);
  206. intruders.push_back(intruder);
  207. --right;
  208. }
  209. }
  210. ++right;
  211. defensive_check_index++;
  212. }
  213. }
  214. template <typename Extremes, typename Intruders, typename SideStrategy>
  215. static inline void get_intruders(Ring const& ring,
  216. Extremes const& extremes,
  217. Intruders& intruders,
  218. SideStrategy const& strategy)
  219. {
  220. std::size_t const n = boost::size(ring);
  221. if (n >= 3)
  222. {
  223. geometry::ever_circling_range_iterator<Ring const> left(ring);
  224. geometry::ever_circling_range_iterator<Ring const> right(ring);
  225. ++right;
  226. get_intruders(ring, left, right, extremes, intruders, strategy);
  227. }
  228. }
  229. template <typename Iterator, typename SideStrategy>
  230. static inline bool right_turn(Ring const& ring, Iterator it, SideStrategy const& strategy)
  231. {
  232. typename std::iterator_traits<Iterator>::difference_type const index
  233. = std::distance(boost::begin(ring), it);
  234. geometry::ever_circling_range_iterator<Ring const> left(ring);
  235. geometry::ever_circling_range_iterator<Ring const> right(ring);
  236. left += index;
  237. right += index;
  238. if (! extent_both_sides(ring, *it, left, right))
  239. {
  240. return false;
  241. }
  242. int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
  243. int const first_side = strategy.apply(*(right - 1), *right, *left) * factor;
  244. int const last_side = strategy.apply(*left, *(left + 1), *right) * factor;
  245. //std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl;
  246. // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases)
  247. return first_side != 1 && last_side != 1;
  248. }
  249. // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring
  250. template <typename Extremes, typename Intruders, typename SideStrategy>
  251. static inline bool apply(Ring const& ring,
  252. Extremes& extremes,
  253. Intruders& intruders,
  254. SideStrategy const& strategy)
  255. {
  256. std::size_t const n = boost::size(ring);
  257. if (n < 3)
  258. {
  259. return false;
  260. }
  261. // Get all maxima, usually one. In case of self-tangencies, or self-crossings,
  262. // the max might be is not valid. A valid max should make a right turn
  263. range_iterator max_it = boost::begin(ring);
  264. compare<Dimension> smaller;
  265. for (range_iterator it = max_it + 1; it != boost::end(ring); ++it)
  266. {
  267. if (smaller(*max_it, *it) && right_turn(ring, it, strategy))
  268. {
  269. max_it = it;
  270. }
  271. }
  272. if (max_it == boost::end(ring))
  273. {
  274. return false;
  275. }
  276. typename std::iterator_traits<range_iterator>::difference_type const
  277. index = std::distance(boost::begin(ring), max_it);
  278. //std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl;
  279. geometry::ever_circling_range_iterator<Ring const> left(ring);
  280. geometry::ever_circling_range_iterator<Ring const> right(ring);
  281. left += index;
  282. right += index;
  283. // Collect all points (often 3) in a temporary vector
  284. std::vector<point_type> points;
  285. points.reserve(3);
  286. if (! collect(ring, *max_it, points, left, right))
  287. {
  288. return false;
  289. }
  290. //std::cout << "Built vector of " << points.size() << std::endl;
  291. coordinate_type const front_value = geometry::get<Dimension>(points.front());
  292. coordinate_type const back_value = geometry::get<Dimension>(points.back());
  293. coordinate_type const base_value = (std::max)(front_value, back_value);
  294. if (front_value < back_value)
  295. {
  296. move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value);
  297. }
  298. else
  299. {
  300. move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value);
  301. }
  302. std::copy(points.begin(), points.end(), std::back_inserter(extremes));
  303. get_intruders(ring, left, right, extremes, intruders, strategy);
  304. return true;
  305. }
  306. };
  307. }} // namespace detail::extreme_points
  308. #endif // DOXYGEN_NO_DETAIL
  309. #ifndef DOXYGEN_NO_DISPATCH
  310. namespace dispatch
  311. {
  312. template
  313. <
  314. typename Geometry,
  315. std::size_t Dimension,
  316. typename GeometryTag = typename tag<Geometry>::type
  317. >
  318. struct extreme_points
  319. {};
  320. template<typename Ring, std::size_t Dimension>
  321. struct extreme_points<Ring, Dimension, ring_tag>
  322. : detail::extreme_points::extreme_points_on_ring<Ring, Dimension>
  323. {};
  324. template<typename Polygon, std::size_t Dimension>
  325. struct extreme_points<Polygon, Dimension, polygon_tag>
  326. {
  327. template <typename Extremes, typename Intruders, typename SideStrategy>
  328. static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders,
  329. SideStrategy const& strategy)
  330. {
  331. typedef typename geometry::ring_type<Polygon>::type ring_type;
  332. typedef detail::extreme_points::extreme_points_on_ring
  333. <
  334. ring_type, Dimension
  335. > ring_implementation;
  336. if (! ring_implementation::apply(geometry::exterior_ring(polygon),
  337. extremes, intruders, strategy))
  338. {
  339. return false;
  340. }
  341. // For a polygon, its interior rings can contain intruders
  342. typename interior_return_type<Polygon const>::type
  343. rings = interior_rings(polygon);
  344. for (typename detail::interior_iterator<Polygon const>::type
  345. it = boost::begin(rings); it != boost::end(rings); ++it)
  346. {
  347. ring_implementation::get_intruders(*it, extremes, intruders, strategy);
  348. }
  349. return true;
  350. }
  351. };
  352. template<typename Box>
  353. struct extreme_points<Box, 1, box_tag>
  354. {
  355. template <typename Extremes, typename Intruders, typename SideStrategy>
  356. static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
  357. SideStrategy const& )
  358. {
  359. extremes.resize(4);
  360. geometry::detail::assign_box_corners_oriented<false>(box, extremes);
  361. // ll,ul,ur,lr, contains too exactly the right info
  362. return true;
  363. }
  364. };
  365. template<typename Box>
  366. struct extreme_points<Box, 0, box_tag>
  367. {
  368. template <typename Extremes, typename Intruders, typename SideStrategy>
  369. static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
  370. SideStrategy const& )
  371. {
  372. extremes.resize(4);
  373. geometry::detail::assign_box_corners_oriented<false>(box, extremes);
  374. // ll,ul,ur,lr, rotate one to start with UL and end with LL
  375. std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end());
  376. return true;
  377. }
  378. };
  379. template<typename MultiPolygon, std::size_t Dimension>
  380. struct extreme_points<MultiPolygon, Dimension, multi_polygon_tag>
  381. {
  382. template <typename Extremes, typename Intruders, typename SideStrategy>
  383. static inline bool apply(MultiPolygon const& multi, Extremes& extremes,
  384. Intruders& intruders, SideStrategy const& strategy)
  385. {
  386. // Get one for the very first polygon, that is (for the moment) enough.
  387. // It is not guaranteed the "extreme" then, but for the current purpose
  388. // (point_on_surface) it can just be this point.
  389. if (boost::size(multi) >= 1)
  390. {
  391. return extreme_points
  392. <
  393. typename boost::range_value<MultiPolygon const>::type,
  394. Dimension,
  395. polygon_tag
  396. >::apply(*boost::begin(multi), extremes, intruders, strategy);
  397. }
  398. return false;
  399. }
  400. };
  401. } // namespace dispatch
  402. #endif // DOXYGEN_NO_DISPATCH
  403. /*!
  404. \brief Returns extreme points (for Edge=1 in dimension 1, so the top,
  405. for Edge=0 in dimension 0, the right side)
  406. \note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that.
  407. */
  408. template
  409. <
  410. std::size_t Edge,
  411. typename Geometry,
  412. typename Extremes,
  413. typename Intruders,
  414. typename SideStrategy
  415. >
  416. inline bool extreme_points(Geometry const& geometry,
  417. Extremes& extremes,
  418. Intruders& intruders,
  419. SideStrategy const& strategy)
  420. {
  421. concepts::check<Geometry const>();
  422. // Extremes is not required to follow a geometry concept (but it should support an output iterator),
  423. // but its elements should fulfil the point-concept
  424. concepts::check<typename boost::range_value<Extremes>::type>();
  425. // Intruders should contain collections which value type is point-concept
  426. // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept
  427. concepts::check
  428. <
  429. typename boost::range_value
  430. <
  431. typename boost::range_value<Intruders>::type
  432. >::type
  433. const
  434. >();
  435. return dispatch::extreme_points
  436. <
  437. Geometry,
  438. Edge
  439. >::apply(geometry, extremes, intruders, strategy);
  440. }
  441. template
  442. <
  443. std::size_t Edge,
  444. typename Geometry,
  445. typename Extremes,
  446. typename Intruders
  447. >
  448. inline bool extreme_points(Geometry const& geometry,
  449. Extremes& extremes,
  450. Intruders& intruders)
  451. {
  452. typedef typename strategy::side::services::default_strategy
  453. <
  454. typename cs_tag<Geometry>::type
  455. >::type strategy_type;
  456. return geometry::extreme_points<Edge>(geometry,extremes, intruders, strategy_type());
  457. }
  458. }} // namespace boost::geometry
  459. #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP