// Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP #define BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP #include #include #include #include #include #include #include #include #include #include #include namespace boost { namespace geometry { namespace strategy { namespace convex_hull { #ifndef DOXYGEN_NO_DETAIL namespace detail { template < typename InputRange, typename RangeIterator, typename StrategyLess, typename StrategyGreater > struct get_extremes { typedef typename point_type::type point_type; point_type left, right; bool first; StrategyLess less; StrategyGreater greater; inline get_extremes() : first(true) {} inline void apply(InputRange const& range) { if (boost::size(range) == 0) { return; } // First iterate through this range // (this two-stage approach avoids many point copies, // because iterators are kept in memory. Because iterators are // not persistent (in MSVC) this approach is not applicable // for more ranges together) RangeIterator left_it = boost::begin(range); RangeIterator right_it = boost::begin(range); for (RangeIterator it = boost::begin(range) + 1; it != boost::end(range); ++it) { if (less(*it, *left_it)) { left_it = it; } if (greater(*it, *right_it)) { right_it = it; } } // Then compare with earlier if (first) { // First time, assign left/right left = *left_it; right = *right_it; first = false; } else { // Next time, check if this range was left/right from // the extremes already collected if (less(*left_it, left)) { left = *left_it; } if (greater(*right_it, right)) { right = *right_it; } } } }; template < typename InputRange, typename RangeIterator, typename Container, typename SideStrategy > struct assign_range { Container lower_points, upper_points; typedef typename point_type::type point_type; point_type const& most_left; point_type const& most_right; inline assign_range(point_type const& left, point_type const& right) : most_left(left) , most_right(right) {} inline void apply(InputRange const& range) { typedef SideStrategy side; // Put points in one of the two output sequences for (RangeIterator it = boost::begin(range); it != boost::end(range); ++it) { // check if it is lying most_left or most_right from the line int dir = side::apply(most_left, most_right, *it); switch(dir) { case 1 : // left side upper_points.push_back(*it); break; case -1 : // right side lower_points.push_back(*it); break; // 0: on line most_left-most_right, // or most_left, or most_right, // -> all never part of hull } } } }; template static inline void sort(Range& range) { typedef typename boost::range_value::type point_type; typedef geometry::less comparator; std::sort(boost::begin(range), boost::end(range), comparator()); } } // namespace detail #endif // DOXYGEN_NO_DETAIL /*! \brief Graham scan strategy to calculate convex hull \ingroup strategies \note Completely reworked version inspired on the sources listed below \see http://www.ddj.com/architect/201806315 \see http://marknelson.us/2007/08/22/convex */ template class graham_andrew { public : typedef OutputPoint point_type; typedef InputGeometry geometry_type; private: typedef typename cs_tag::type cs_tag; typedef typename std::vector container_type; typedef typename std::vector::const_iterator iterator; typedef typename std::vector::const_reverse_iterator rev_iterator; class partitions { friend class graham_andrew; container_type m_lower_hull; container_type m_upper_hull; container_type m_copied_input; }; public: typedef partitions state_type; inline void apply(InputGeometry const& geometry, partitions& state) const { // First pass. // Get min/max (in most cases left / right) points // This makes use of the geometry::less/greater predicates // For the left boundary it is important that multiple points // are sorted from bottom to top. Therefore the less predicate // does not take the x-only template parameter (this fixes ticket #6019. // For the right boundary it is not necessary (though also not harmful), // because points are sorted from bottom to top in a later stage. // For symmetry and to get often more balanced lower/upper halves // we keep it. typedef typename geometry::detail::range_type::type range_type; typedef typename boost::range_iterator < range_type const >::type range_iterator; detail::get_extremes < range_type, range_iterator, geometry::less, geometry::greater > extremes; geometry::detail::for_each_range(geometry, extremes); // Bounding left/right points // Second pass, now that extremes are found, assign all points // in either lower, either upper detail::assign_range < range_type, range_iterator, container_type, typename strategy::side::services::default_strategy::type > assigner(extremes.left, extremes.right); geometry::detail::for_each_range(geometry, assigner); // Sort both collections, first on x(, then on y) detail::sort(assigner.lower_points); detail::sort(assigner.upper_points); //std::cout << boost::size(assigner.lower_points) << std::endl; //std::cout << boost::size(assigner.upper_points) << std::endl; // And decide which point should be in the final hull build_half_hull<-1>(assigner.lower_points, state.m_lower_hull, extremes.left, extremes.right); build_half_hull<1>(assigner.upper_points, state.m_upper_hull, extremes.left, extremes.right); } template inline void result(partitions const& state, OutputIterator out, bool clockwise) const { if (clockwise) { output_range(state.m_upper_hull, out, false); output_range(state.m_lower_hull, out, true); } else { output_range(state.m_lower_hull, out, false); output_range(state.m_upper_hull, out, true); } } private: template static inline void build_half_hull(container_type const& input, container_type& output, point_type const& left, point_type const& right) { output.push_back(left); for(iterator it = input.begin(); it != input.end(); ++it) { add_to_hull(*it, output); } add_to_hull(right, output); } template static inline void add_to_hull(point_type const& p, container_type& output) { typedef typename strategy::side::services::default_strategy::type side; output.push_back(p); register std::size_t output_size = output.size(); while (output_size >= 3) { rev_iterator rit = output.rbegin(); point_type const& last = *rit++; point_type const& last2 = *rit++; if (Factor * side::apply(*rit, last, last2) <= 0) { // Remove last two points from stack, and add last again // This is much faster then erasing the one but last. output.pop_back(); output.pop_back(); output.push_back(last); output_size--; } else { return; } } } template static inline void output_range(container_type const& range, OutputIterator out, bool skip_first) { typedef typename reversible_view::type view_type; view_type view(range); bool first = true; for (typename boost::range_iterator::type it = boost::begin(view); it != boost::end(view); ++it) { if (first && skip_first) { first = false; } else { *out = *it; ++out; } } } }; }} // namespace strategy::convex_hull #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS template struct strategy_convex_hull { typedef strategy::convex_hull::graham_andrew type; }; #endif }} // namespace boost::geometry #endif // BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_CONVEX_GRAHAM_ANDREW_HPP