hull_graham_andrew.hpp 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
  3. // This file was modified by Oracle on 2014.
  4. // Modifications copyright (c) 2014 Oracle and/or its affiliates.
  5. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  6. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  7. // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  8. // Use, modification and distribution is subject to the Boost Software License,
  9. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  10. // http://www.boost.org/LICENSE_1_0.txt)
  11. #ifndef BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP
  12. #define BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP
  13. #include <cstddef>
  14. #include <algorithm>
  15. #include <vector>
  16. #include <boost/range.hpp>
  17. #include <boost/geometry/core/assert.hpp>
  18. #include <boost/geometry/core/cs.hpp>
  19. #include <boost/geometry/core/point_type.hpp>
  20. #include <boost/geometry/strategies/convex_hull.hpp>
  21. #include <boost/geometry/views/detail/range_type.hpp>
  22. #include <boost/geometry/policies/compare.hpp>
  23. #include <boost/geometry/algorithms/detail/for_each_range.hpp>
  24. #include <boost/geometry/views/reversible_view.hpp>
  25. namespace boost { namespace geometry
  26. {
  27. namespace strategy { namespace convex_hull
  28. {
  29. #ifndef DOXYGEN_NO_DETAIL
  30. namespace detail
  31. {
  32. template
  33. <
  34. typename InputRange,
  35. typename RangeIterator,
  36. typename StrategyLess,
  37. typename StrategyGreater
  38. >
  39. struct get_extremes
  40. {
  41. typedef typename point_type<InputRange>::type point_type;
  42. point_type left, right;
  43. bool first;
  44. StrategyLess less;
  45. StrategyGreater greater;
  46. inline get_extremes()
  47. : first(true)
  48. {}
  49. inline void apply(InputRange const& range)
  50. {
  51. if (boost::size(range) == 0)
  52. {
  53. return;
  54. }
  55. // First iterate through this range
  56. // (this two-stage approach avoids many point copies,
  57. // because iterators are kept in memory. Because iterators are
  58. // not persistent (in MSVC) this approach is not applicable
  59. // for more ranges together)
  60. RangeIterator left_it = boost::begin(range);
  61. RangeIterator right_it = boost::begin(range);
  62. for (RangeIterator it = boost::begin(range) + 1;
  63. it != boost::end(range);
  64. ++it)
  65. {
  66. if (less(*it, *left_it))
  67. {
  68. left_it = it;
  69. }
  70. if (greater(*it, *right_it))
  71. {
  72. right_it = it;
  73. }
  74. }
  75. // Then compare with earlier
  76. if (first)
  77. {
  78. // First time, assign left/right
  79. left = *left_it;
  80. right = *right_it;
  81. first = false;
  82. }
  83. else
  84. {
  85. // Next time, check if this range was left/right from
  86. // the extremes already collected
  87. if (less(*left_it, left))
  88. {
  89. left = *left_it;
  90. }
  91. if (greater(*right_it, right))
  92. {
  93. right = *right_it;
  94. }
  95. }
  96. }
  97. };
  98. template
  99. <
  100. typename InputRange,
  101. typename RangeIterator,
  102. typename Container,
  103. typename SideStrategy
  104. >
  105. struct assign_range
  106. {
  107. Container lower_points, upper_points;
  108. typedef typename point_type<InputRange>::type point_type;
  109. point_type const& most_left;
  110. point_type const& most_right;
  111. inline assign_range(point_type const& left, point_type const& right)
  112. : most_left(left)
  113. , most_right(right)
  114. {}
  115. inline void apply(InputRange const& range)
  116. {
  117. typedef SideStrategy side;
  118. // Put points in one of the two output sequences
  119. for (RangeIterator it = boost::begin(range);
  120. it != boost::end(range);
  121. ++it)
  122. {
  123. // check if it is lying most_left or most_right from the line
  124. int dir = side::apply(most_left, most_right, *it);
  125. switch(dir)
  126. {
  127. case 1 : // left side
  128. upper_points.push_back(*it);
  129. break;
  130. case -1 : // right side
  131. lower_points.push_back(*it);
  132. break;
  133. // 0: on line most_left-most_right,
  134. // or most_left, or most_right,
  135. // -> all never part of hull
  136. }
  137. }
  138. }
  139. };
  140. template <typename Range>
  141. static inline void sort(Range& range)
  142. {
  143. typedef typename boost::range_value<Range>::type point_type;
  144. typedef geometry::less<point_type> comparator;
  145. std::sort(boost::begin(range), boost::end(range), comparator());
  146. }
  147. } // namespace detail
  148. #endif // DOXYGEN_NO_DETAIL
  149. /*!
  150. \brief Graham scan strategy to calculate convex hull
  151. \ingroup strategies
  152. \note Completely reworked version inspired on the sources listed below
  153. \see http://www.ddj.com/architect/201806315
  154. \see http://marknelson.us/2007/08/22/convex
  155. */
  156. template <typename InputGeometry, typename OutputPoint>
  157. class graham_andrew
  158. {
  159. public :
  160. typedef OutputPoint point_type;
  161. typedef InputGeometry geometry_type;
  162. private:
  163. typedef typename cs_tag<point_type>::type cs_tag;
  164. typedef typename std::vector<point_type> container_type;
  165. typedef typename std::vector<point_type>::const_iterator iterator;
  166. typedef typename std::vector<point_type>::const_reverse_iterator rev_iterator;
  167. class partitions
  168. {
  169. friend class graham_andrew;
  170. container_type m_lower_hull;
  171. container_type m_upper_hull;
  172. container_type m_copied_input;
  173. };
  174. public:
  175. typedef partitions state_type;
  176. inline void apply(InputGeometry const& geometry, partitions& state) const
  177. {
  178. // First pass.
  179. // Get min/max (in most cases left / right) points
  180. // This makes use of the geometry::less/greater predicates
  181. // For the left boundary it is important that multiple points
  182. // are sorted from bottom to top. Therefore the less predicate
  183. // does not take the x-only template parameter (this fixes ticket #6019.
  184. // For the right boundary it is not necessary (though also not harmful),
  185. // because points are sorted from bottom to top in a later stage.
  186. // For symmetry and to get often more balanced lower/upper halves
  187. // we keep it.
  188. typedef typename geometry::detail::range_type<InputGeometry>::type range_type;
  189. typedef typename boost::range_iterator
  190. <
  191. range_type const
  192. >::type range_iterator;
  193. detail::get_extremes
  194. <
  195. range_type,
  196. range_iterator,
  197. geometry::less<point_type>,
  198. geometry::greater<point_type>
  199. > extremes;
  200. geometry::detail::for_each_range(geometry, extremes);
  201. // Bounding left/right points
  202. // Second pass, now that extremes are found, assign all points
  203. // in either lower, either upper
  204. detail::assign_range
  205. <
  206. range_type,
  207. range_iterator,
  208. container_type,
  209. typename strategy::side::services::default_strategy<cs_tag>::type
  210. > assigner(extremes.left, extremes.right);
  211. geometry::detail::for_each_range(geometry, assigner);
  212. // Sort both collections, first on x(, then on y)
  213. detail::sort(assigner.lower_points);
  214. detail::sort(assigner.upper_points);
  215. //std::cout << boost::size(assigner.lower_points) << std::endl;
  216. //std::cout << boost::size(assigner.upper_points) << std::endl;
  217. // And decide which point should be in the final hull
  218. build_half_hull<-1>(assigner.lower_points, state.m_lower_hull,
  219. extremes.left, extremes.right);
  220. build_half_hull<1>(assigner.upper_points, state.m_upper_hull,
  221. extremes.left, extremes.right);
  222. }
  223. template <typename OutputIterator>
  224. inline void result(partitions const& state,
  225. OutputIterator out,
  226. bool clockwise,
  227. bool closed) const
  228. {
  229. if (clockwise)
  230. {
  231. output_ranges(state.m_upper_hull, state.m_lower_hull, out, closed);
  232. }
  233. else
  234. {
  235. output_ranges(state.m_lower_hull, state.m_upper_hull, out, closed);
  236. }
  237. }
  238. private:
  239. template <int Factor>
  240. static inline void build_half_hull(container_type const& input,
  241. container_type& output,
  242. point_type const& left, point_type const& right)
  243. {
  244. output.push_back(left);
  245. for(iterator it = input.begin(); it != input.end(); ++it)
  246. {
  247. add_to_hull<Factor>(*it, output);
  248. }
  249. add_to_hull<Factor>(right, output);
  250. }
  251. template <int Factor>
  252. static inline void add_to_hull(point_type const& p, container_type& output)
  253. {
  254. typedef typename strategy::side::services::default_strategy<cs_tag>::type side;
  255. output.push_back(p);
  256. std::size_t output_size = output.size();
  257. while (output_size >= 3)
  258. {
  259. rev_iterator rit = output.rbegin();
  260. point_type const last = *rit++;
  261. point_type const& last2 = *rit++;
  262. if (Factor * side::apply(*rit, last, last2) <= 0)
  263. {
  264. // Remove last two points from stack, and add last again
  265. // This is much faster then erasing the one but last.
  266. output.pop_back();
  267. output.pop_back();
  268. output.push_back(last);
  269. output_size--;
  270. }
  271. else
  272. {
  273. return;
  274. }
  275. }
  276. }
  277. template <typename OutputIterator>
  278. static inline void output_ranges(container_type const& first, container_type const& second,
  279. OutputIterator out, bool closed)
  280. {
  281. std::copy(boost::begin(first), boost::end(first), out);
  282. BOOST_GEOMETRY_ASSERT(closed ? !boost::empty(second) : boost::size(second) > 1);
  283. std::copy(++boost::rbegin(second), // skip the first Point
  284. closed ? boost::rend(second) : --boost::rend(second), // skip the last Point if open
  285. out);
  286. typedef typename boost::range_size<container_type>::type size_type;
  287. size_type const count = boost::size(first) + boost::size(second) - 1;
  288. // count describes a closed case but comparison with min size of closed
  289. // gives the result compatible also with open
  290. // here core_detail::closure::minimum_ring_size<closed> could be used
  291. if (count < 4)
  292. {
  293. // there should be only one missing
  294. *out++ = *boost::begin(first);
  295. }
  296. }
  297. };
  298. }} // namespace strategy::convex_hull
  299. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  300. template <typename InputGeometry, typename OutputPoint>
  301. struct strategy_convex_hull<InputGeometry, OutputPoint, cartesian_tag>
  302. {
  303. typedef strategy::convex_hull::graham_andrew<InputGeometry, OutputPoint> type;
  304. };
  305. #endif
  306. }} // namespace boost::geometry
  307. #endif // BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP