hull_graham_andrew.hpp 11 KB

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