detail/convex_hull/graham_andrew.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2014-2023.
0006 // Modifications copyright (c) 2014-2023 Oracle and/or its affiliates.
0007
0008 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0009 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0010
0011 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
0012 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 #ifndef BOOST_GEOMETRY_ALGORITHMS_CONVEX_HULL_GRAHAM_ANDREW_HPP
0019 #define BOOST_GEOMETRY_ALGORITHMS_CONVEX_HULL_GRAHAM_ANDREW_HPP
0020
0021
0022 #include <cstddef>
0023 #include <algorithm>
0024 #include <vector>
0025
0026 #include <boost/range/size.hpp>
0027
0028 #include <boost/geometry/algorithms/detail/for_each_range.hpp>
0029 #include <boost/geometry/core/assert.hpp>
0030 #include <boost/geometry/core/closure.hpp>
0031 #include <boost/geometry/core/cs.hpp>
0032 #include <boost/geometry/core/point_type.hpp>
0033 #include <boost/geometry/core/point_order.hpp>
0034 #include <boost/geometry/policies/compare.hpp>
0035 #include <boost/geometry/strategies/convex_hull/cartesian.hpp>
0036 #include <boost/geometry/strategies/convex_hull/geographic.hpp>
0037 #include <boost/geometry/strategies/convex_hull/spherical.hpp>
0038
0039 #include <boost/geometry/util/range.hpp>
0040
0041
0042 namespace boost { namespace geometry
0043 {
0044
0045 #ifndef DOXYGEN_NO_DETAIL
0046 namespace detail { namespace convex_hull
0047 {
0048
0049 // TODO: All of the copies could be avoided if this function stored pointers to points.
0050 //       But would it be possible considering that a range can return proxy reference?
0051 template <typename InputProxy, typename Point, typename Less>
0052 inline void get_extremes(InputProxy const& in_proxy,
0053                          Point& left, Point& right,
0054                          Less const& less)
0055 {
0056     bool first = true;
0057     in_proxy.for_each_range([&](auto const& range)
0058     {
0059         if (boost::empty(range))
0060         {
0061             return;
0062         }
0063
0064         // First iterate through this range
0065         // (this two-stage approach avoids many point copies,
0066         //  because iterators are kept in memory. Because iterators are
0067         //  not persistent (in MSVC) this approach is not applicable
0068         //  for more ranges together)
0069
0070         auto left_it = boost::begin(range);
0071         auto right_it = boost::begin(range);
0072
0073         auto it = boost::begin(range);
0074         for (++it; it != boost::end(range); ++it)
0075         {
0076             if (less(*it, *left_it))
0077             {
0078                 left_it = it;
0079             }
0080
0081             if (less(*right_it, *it))
0082             {
0083                 right_it = it;
0084             }
0085         }
0086
0087         // Then compare with earlier
0088         if (first)
0089         {
0090             // First time, assign left/right
0091             left = *left_it;
0092             right = *right_it;
0093             first = false;
0094         }
0095         else
0096         {
0097             // Next time, check if this range was left/right from
0098             // the extremes already collected
0099             if (less(*left_it, left))
0100             {
0101                 left = *left_it;
0102             }
0103
0104             if (less(right, *right_it))
0105             {
0106                 right = *right_it;
0107             }
0108         }
0109     });
0110 }
0111
0112
0113 template <typename InputProxy, typename Point, typename Container, typename SideStrategy>
0114 inline void assign_ranges(InputProxy const& in_proxy,
0115                           Point const& most_left, Point const& most_right,
0116                           Container& lower_points, Container& upper_points,
0117                           SideStrategy const& side)
0118 {
0119     in_proxy.for_each_range([&](auto const& range)
0120     {
0121         // Put points in one of the two output sequences
0122         for (auto it = boost::begin(range); it != boost::end(range); ++it)
0123         {
0124             // check if it is lying most_left or most_right from the line
0125
0126             int dir = side.apply(most_left, most_right, *it);
0127             switch(dir)
0128             {
0129                 case 1 : // left side
0130                     upper_points.push_back(*it);
0131                     break;
0132                 case -1 : // right side
0133                     lower_points.push_back(*it);
0134                     break;
0135
0136                 // 0: on line most_left-most_right,
0137                 //    or most_left, or most_right,
0138                 //    -> all never part of hull
0139             }
0140         }
0141     });
0142 }
0143
0144
0145 /*!
0146 \brief Graham scan algorithm to calculate convex hull
0147  */
0148 template <typename InputPoint>
0149 class graham_andrew
0150 {
0151     typedef InputPoint point_type;
0152     typedef typename std::vector<point_type> container_type;
0153
0154     class partitions
0155     {
0156         friend class graham_andrew;
0157
0158         container_type m_lower_hull;
0159         container_type m_upper_hull;
0160         container_type m_copied_input;
0161     };
0162
0163 public:
0164     template <typename InputProxy, typename OutputRing, typename Strategy>
0165     static void apply(InputProxy const& in_proxy, OutputRing & out_ring, Strategy& strategy)
0166     {
0167         partitions state;
0168
0169         apply(in_proxy, state, strategy);
0170
0171         result(state,
0172                range::back_inserter(out_ring),
0173                geometry::point_order<OutputRing>::value == clockwise,
0174                geometry::closure<OutputRing>::value != open);
0175     }
0176
0177 private:
0178     template <typename InputProxy, typename Strategy>
0179     static void apply(InputProxy const& in_proxy, partitions& state, Strategy& strategy)
0180     {
0181         // First pass.
0182         // Get min/max (in most cases left / right) points
0183         // This makes use of the geometry::less/greater predicates
0184
0185         // For the left boundary it is important that multiple points
0186         // are sorted from bottom to top. Therefore the less predicate
0187         // does not take the x-only template parameter (this fixes ticket #6019.
0188         // For the right boundary it is not necessary (though also not harmful),
0189         // because points are sorted from bottom to top in a later stage.
0190         // For symmetry and to get often more balanced lower/upper halves
0191         // we keep it.
0192
0193         point_type most_left, most_right;
0194
0195         geometry::less_exact<point_type, -1, Strategy> less;
0196
0197         detail::convex_hull::get_extremes(in_proxy, most_left, most_right, less);
0198
0199         container_type lower_points, upper_points;
0200
0201         auto const side_strategy = strategy.side();
0202
0203         // Bounding left/right points
0204         // Second pass, now that extremes are found, assign all points
0205         // in either lower, either upper
0206         detail::convex_hull::assign_ranges(in_proxy, most_left, most_right,
0207                               lower_points, upper_points,
0208                               side_strategy);
0209
0210         // Sort both collections, first on x(, then on y)
0211         std::sort(boost::begin(lower_points), boost::end(lower_points), less);
0212         std::sort(boost::begin(upper_points), boost::end(upper_points), less);
0213
0214         // And decide which point should be in the final hull
0215         build_half_hull<-1>(lower_points, state.m_lower_hull,
0216                             most_left, most_right,
0217                             side_strategy);
0218         build_half_hull<1>(upper_points, state.m_upper_hull,
0219                            most_left, most_right,
0220                            side_strategy);
0221     }
0222
0223     template <int Factor, typename SideStrategy>
0224     static inline void build_half_hull(container_type const& input,
0225             container_type& output,
0226             point_type const& left, point_type const& right,
0227             SideStrategy const& side)
0228     {
0229         output.push_back(left);
0230         for (auto const& i : input)
0231         {
0232             add_to_hull<Factor>(i, output, side);
0233         }
0234         add_to_hull<Factor>(right, output, side);
0235     }
0236
0237
0238     template <int Factor, typename SideStrategy>
0239     static inline void add_to_hull(point_type const& p, container_type& output,
0240                                    SideStrategy const& side)
0241     {
0242         output.push_back(p);
0243         std::size_t output_size = output.size();
0244         while (output_size >= 3)
0245         {
0246             auto rit = output.rbegin();
0247             point_type const last = *rit++;
0248             point_type const& last2 = *rit++;
0249
0250             if (Factor * side.apply(*rit, last, last2) <= 0)
0251             {
0252                 // Remove last two points from stack, and add last again
0253                 // This is much faster then erasing the one but last.
0254                 output.pop_back();
0255                 output.pop_back();
0256                 output.push_back(last);
0257                 output_size--;
0258             }
0259             else
0260             {
0261                 return;
0262             }
0263         }
0264     }
0265
0266
0267     template <typename OutputIterator>
0268     static void result(partitions const& state, OutputIterator out, bool clockwise, bool closed)
0269     {
0270         if (clockwise)
0271         {
0272             output_ranges(state.m_upper_hull, state.m_lower_hull, out, closed);
0273         }
0274         else
0275         {
0276             output_ranges(state.m_lower_hull, state.m_upper_hull, out, closed);
0277         }
0278     }
0279
0280     template <typename OutputIterator>
0281     static inline void output_ranges(container_type const& first,
0282                                      container_type const& second,
0283                                      OutputIterator out,
0284                                      bool closed)
0285     {
0286         std::copy(boost::begin(first), boost::end(first), out);
0287
0288         BOOST_GEOMETRY_ASSERT(closed ? !boost::empty(second) : boost::size(second) > 1);
0289         std::copy(++boost::rbegin(second), // skip the first Point
0290                   closed ? boost::rend(second) : --boost::rend(second), // skip the last Point if open
0291                   out);
0292
0293         typedef typename boost::range_size<container_type>::type size_type;
0294         size_type const count = boost::size(first) + boost::size(second) - 1;
0295         // count describes a closed case but comparison with min size of closed
0296         // gives the result compatible also with open
0297         // here core_detail::closure::minimum_ring_size<closed> could be used
0298         if (count < 4)
0299         {
0300             // there should be only one missing
0301             *out++ = *boost::begin(first);
0302         }
0303     }
0304 };
0305
0306
0307 }} // namespace detail::convex_hull
0308 #endif // DOXYGEN_NO_DETAIL
0309
0310 }} // namespace boost::geometry
0311
0312 #endif // BOOST_GEOMETRY_ALGORITHMS_CONVEX_HULL_GRAHAM_ANDREW_HPP