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