algorithms/detail/extreme_points.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
0005 // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
0006 // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
0007
0008 // This file was modified by Oracle on 2017-2021.
0009 // Modifications copyright (c) 2017-2021 Oracle and/or its affiliates.
0010
0011 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0012
0013 // Use, modification and distribution is subject to the Boost Software License,
0014 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0015 // http://www.boost.org/LICENSE_1_0.txt)
0016
0017 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
0018 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
0019
0020
0021 #include <cstddef>
0022
0023 #include <boost/range/begin.hpp>
0024 #include <boost/range/end.hpp>
0025 #include <boost/range/size.hpp>
0026 #include <boost/range/value_type.hpp>
0027
0028 #include <boost/geometry/algorithms/detail/assign_box_corners.hpp>
0029
0030 #include <boost/geometry/arithmetic/arithmetic.hpp>
0031
0032 #include <boost/geometry/core/cs.hpp>
0033 #include <boost/geometry/core/point_order.hpp>
0034 #include <boost/geometry/core/point_type.hpp>
0035 #include <boost/geometry/core/ring_type.hpp>
0036 #include <boost/geometry/core/tags.hpp>
0037
0038 #include <boost/geometry/geometries/concepts/check.hpp>
0039 #include <boost/geometry/iterators/ever_circling_iterator.hpp>
0040
0041 #include <boost/geometry/strategies/side.hpp>
0042
0043 #include <boost/geometry/util/math.hpp>
0044
0045
0046 namespace boost { namespace geometry
0047 {
0048
0049
0050 #ifndef DOXYGEN_NO_DETAIL
0051 namespace detail { namespace extreme_points
0052 {
0053
0054 template <std::size_t Dimension>
0055 struct compare
0056 {
0057     template <typename Point>
0058     inline bool operator()(Point const& lhs, Point const& rhs)
0059     {
0060         return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs);
0061     }
0062 };
0063
0064
0065 template <std::size_t Dimension, typename PointType, typename CoordinateType>
0066 inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value)
0067 {
0068     // Moves a point along the vector (point, extreme) in the direction of the extreme  point
0069     // This adapts the possibly uneven legs of the triangle (or trapezium-like shape)
0070     //      _____extreme            _____
0071     //     /     \                 /     \                                    .
0072     //    /base   \         =>    /       \ point                             .
0073     //             \ point                                                    .
0074     //
0075     // For so-called intruders, it can be used to adapt both legs to the level of "base"
0076     // For the base, it can be used to adapt both legs to the level of the max-value of the intruders
0077     // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector
0078
0079     CoordinateType const value = geometry::get<Dimension>(point);
0080     //if (geometry::math::equals(value, base_value))
0081     if (value >= base_value)
0082     {
0083         return;
0084     }
0085
0086     PointType vector = point;
0087     subtract_point(vector, extreme);
0088
0089     CoordinateType const diff = geometry::get<Dimension>(vector);
0090
0091     // diff should never be zero
0092     // because of the way our triangle/trapezium is build.
0093     // We just return if it would be the case.
0094     if (geometry::math::equals(diff, 0))
0095     {
0096         return;
0097     }
0098
0099     CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme);
0100
0101     multiply_value(vector, base_diff);
0102     divide_value(vector, diff);
0103
0104     // The real move:
0105     point = extreme;
0106     add_point(point, vector);
0107 }
0108
0109
0110 template <std::size_t Dimension, typename Range, typename CoordinateType>
0111 inline void move_along_vector(Range& range, CoordinateType const& base_value)
0112 {
0113     if (range.size() >= 3)
0114     {
0115         move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value);
0116         move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value);
0117     }
0118 }
0119
0120
0121 template<typename Ring, std::size_t Dimension>
0122 struct extreme_points_on_ring
0123 {
0124
0125     typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
0126     typedef typename geometry::point_type<Ring>::type point_type;
0127
0128     template <typename CirclingIterator, typename Points>
0129     static inline bool extend(CirclingIterator& it,
0130             std::size_t n,
0131             coordinate_type max_coordinate_value,
0132             Points& points, int direction)
0133     {
0134         std::size_t safe_index = 0;
0135         do
0136         {
0137             it += direction;
0138
0139             points.push_back(*it);
0140
0141             if (safe_index++ >= n)
0142             {
0143                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
0144                 return false;
0145             }
0146         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
0147
0148         return true;
0149     }
0150
0151     // Overload without adding to poinst
0152     template <typename CirclingIterator>
0153     static inline bool extend(CirclingIterator& it,
0154             std::size_t n,
0155             coordinate_type max_coordinate_value,
0156             int direction)
0157     {
0158         std::size_t safe_index = 0;
0159         do
0160         {
0161             it += direction;
0162
0163             if (safe_index++ >= n)
0164             {
0165                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
0166                 return false;
0167             }
0168         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
0169
0170         return true;
0171     }
0172
0173     template <typename CirclingIterator>
0174     static inline bool extent_both_sides(Ring const& ring,
0175             point_type extreme,
0176             CirclingIterator& left,
0177             CirclingIterator& right)
0178     {
0179         std::size_t const n = boost::size(ring);
0180         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
0181
0182         if (! extend(left, n, max_coordinate_value, -1))
0183         {
0184             return false;
0185         }
0186         if (! extend(right, n, max_coordinate_value, +1))
0187         {
0188             return false;
0189         }
0190
0191         return true;
0192     }
0193
0194     template <typename Collection, typename CirclingIterator>
0195     static inline bool collect(Ring const& ring,
0196             point_type extreme,
0197             Collection& points,
0198             CirclingIterator& left,
0199             CirclingIterator& right)
0200     {
0201         std::size_t const n = boost::size(ring);
0202         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
0203
0204         // Collects first left, which is reversed (if more than one point) then adds the top itself, then right
0205         if (! extend(left, n, max_coordinate_value, points, -1))
0206         {
0207             return false;
0208         }
0209         std::reverse(points.begin(), points.end());
0210         points.push_back(extreme);
0211         if (! extend(right, n, max_coordinate_value, points, +1))
0212         {
0213             return false;
0214         }
0215
0216         return true;
0217     }
0218
0219     template <typename Extremes, typename Intruders, typename CirclingIterator, typename SideStrategy>
0220     static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right,
0221             Extremes const& extremes,
0222             Intruders& intruders,
0223             SideStrategy const& strategy)
0224     {
0225         if (boost::size(extremes) < 3)
0226         {
0227             return;
0228         }
0229         coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>()));
0230
0231         // Also select left/right (if Dimension=1)
0232         coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
0233         coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
0234
0235         std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too
0236         std::size_t const n = boost::size(ring);
0237         while (left != right && defensive_check_index < n)
0238         {
0239             coordinate_type const coordinate = geometry::get<Dimension>(*right);
0240             coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right);
0241             if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max)
0242             {
0243                 int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
0244                 int const first_side = strategy.apply(*right, extremes.front(), *(extremes.begin() + 1)) * factor;
0245                 int const last_side = strategy.apply(*right, *(extremes.rbegin() + 1), extremes.back()) * factor;
0246
0247                 // If not lying left from any of the extemes side
0248                 if (first_side != 1 && last_side != 1)
0249                 {
0250                     //std::cout << "first " << first_side << " last " << last_side << std::endl;
0251
0252                     // we start at this intrusion until it is handled, and don't affect our initial left iterator
0253                     CirclingIterator left_intrusion_it = right;
0254                     typename boost::range_value<Intruders>::type intruder;
0255                     collect(ring, *right, intruder, left_intrusion_it, right);
0256
0257                     // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies
0258                     // (we might postpone this action, it is often not necessary. However it is not time-consuming)
0259                     move_along_vector<Dimension>(intruder, min_value);
0260                     intruders.push_back(intruder);
0261                     --right;
0262                 }
0263             }
0264             ++right;
0265             defensive_check_index++;
0266         }
0267     }
0268
0269     template <typename Extremes, typename Intruders, typename SideStrategy>
0270     static inline void get_intruders(Ring const& ring,
0271             Extremes const& extremes,
0272             Intruders& intruders,
0273             SideStrategy const& strategy)
0274     {
0275         std::size_t const n = boost::size(ring);
0276         if (n >= 3)
0277         {
0278             geometry::ever_circling_range_iterator<Ring const> left(ring);
0279             geometry::ever_circling_range_iterator<Ring const> right(ring);
0280             ++right;
0281
0282             get_intruders(ring, left, right, extremes, intruders, strategy);
0283         }
0284     }
0285
0286     template <typename Iterator, typename SideStrategy>
0287     static inline bool right_turn(Ring const& ring, Iterator it, SideStrategy const& strategy)
0288     {
0289         auto const index = std::distance(boost::begin(ring), it);
0290         geometry::ever_circling_range_iterator<Ring const> left(ring);
0291         geometry::ever_circling_range_iterator<Ring const> right(ring);
0292         left += index;
0293         right += index;
0294
0295         if (! extent_both_sides(ring, *it, left, right))
0296         {
0297             return false;
0298         }
0299
0300         int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
0301         int const first_side = strategy.apply(*(right - 1), *right, *left) * factor;
0302         int const last_side = strategy.apply(*left, *(left + 1), *right) * factor;
0303
0304 //std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl;
0305
0306         // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases)
0307         return first_side != 1 && last_side != 1;
0308     }
0309
0310
0311     // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring
0312     template <typename Extremes, typename Intruders, typename SideStrategy>
0313     static inline bool apply(Ring const& ring,
0314                              Extremes& extremes,
0315                              Intruders& intruders,
0316                              SideStrategy const& strategy)
0317     {
0318         std::size_t const n = boost::size(ring);
0319         if (n < 3)
0320         {
0321             return false;
0322         }
0323
0324         // Get all maxima, usually one. In case of self-tangencies, or self-crossings,
0325         // the max might be is not valid. A valid max should make a right turn
0326         auto max_it = boost::begin(ring);
0327         compare<Dimension> smaller;
0328         for (auto it = max_it + 1; it != boost::end(ring); ++it)
0329         {
0330             if (smaller(*max_it, *it) && right_turn(ring, it, strategy))
0331             {
0332                 max_it = it;
0333             }
0334         }
0335
0336         if (max_it == boost::end(ring))
0337         {
0338             return false;
0339         }
0340
0341         auto const index = std::distance(boost::begin(ring), max_it);
0342 //std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl;
0343
0344         geometry::ever_circling_range_iterator<Ring const> left(ring);
0345         geometry::ever_circling_range_iterator<Ring const> right(ring);
0346         left += index;
0347         right += index;
0348
0349         // Collect all points (often 3) in a temporary vector
0350         std::vector<point_type> points;
0351         points.reserve(3);
0352         if (! collect(ring, *max_it, points, left, right))
0353         {
0354             return false;
0355         }
0356
0357 //std::cout << "Built vector of " << points.size() << std::endl;
0358
0359         coordinate_type const front_value = geometry::get<Dimension>(points.front());
0360         coordinate_type const back_value = geometry::get<Dimension>(points.back());
0361         coordinate_type const base_value = (std::max)(front_value, back_value);
0362         if (front_value < back_value)
0363         {
0364             move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value);
0365         }
0366         else
0367         {
0368             move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value);
0369         }
0370
0371         std::copy(points.begin(), points.end(), std::back_inserter(extremes));
0372
0373         get_intruders(ring, left, right, extremes, intruders, strategy);
0374
0375         return true;
0376     }
0377 };
0378
0379
0380
0381
0382
0383 }} // namespace detail::extreme_points
0384 #endif // DOXYGEN_NO_DETAIL
0385
0386
0387 #ifndef DOXYGEN_NO_DISPATCH
0388 namespace dispatch
0389 {
0390
0391
0392 template
0393 <
0394     typename Geometry,
0395     std::size_t Dimension,
0396     typename GeometryTag = typename tag<Geometry>::type
0397 >
0398 struct extreme_points
0399 {};
0400
0401
0402 template<typename Ring, std::size_t Dimension>
0403 struct extreme_points<Ring, Dimension, ring_tag>
0404     : detail::extreme_points::extreme_points_on_ring<Ring, Dimension>
0405 {};
0406
0407
0408 template<typename Polygon, std::size_t Dimension>
0409 struct extreme_points<Polygon, Dimension, polygon_tag>
0410 {
0411     template <typename Extremes, typename Intruders, typename SideStrategy>
0412     static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders,
0413                              SideStrategy const& strategy)
0414     {
0415         typedef typename geometry::ring_type<Polygon>::type ring_type;
0416         typedef detail::extreme_points::extreme_points_on_ring
0417             <
0418                 ring_type, Dimension
0419             > ring_implementation;
0420
0421         if (! ring_implementation::apply(geometry::exterior_ring(polygon),
0422                                          extremes, intruders, strategy))
0423         {
0424             return false;
0425         }
0426
0427         // For a polygon, its interior rings can contain intruders
0428         auto const& rings = interior_rings(polygon);
0429         for (auto it = boost::begin(rings); it != boost::end(rings); ++it)
0430         {
0431             ring_implementation::get_intruders(*it, extremes,  intruders, strategy);
0432         }
0433
0434         return true;
0435     }
0436 };
0437
0438 template<typename Box>
0439 struct extreme_points<Box, 1, box_tag>
0440 {
0441     template <typename Extremes, typename Intruders, typename SideStrategy>
0442     static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
0443                              SideStrategy const& )
0444     {
0445         extremes.resize(4);
0446         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
0447         // ll,ul,ur,lr, contains too exactly the right info
0448         return true;
0449     }
0450 };
0451
0452 template<typename Box>
0453 struct extreme_points<Box, 0, box_tag>
0454 {
0455     template <typename Extremes, typename Intruders, typename SideStrategy>
0456     static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
0457                              SideStrategy const& )
0458     {
0459         extremes.resize(4);
0460         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
0461         // ll,ul,ur,lr, rotate one to start with UL and end with LL
0462         std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end());
0463         return true;
0464     }
0465 };
0466
0467 template<typename MultiPolygon, std::size_t Dimension>
0468 struct extreme_points<MultiPolygon, Dimension, multi_polygon_tag>
0469 {
0470     template <typename Extremes, typename Intruders, typename SideStrategy>
0471     static inline bool apply(MultiPolygon const& multi, Extremes& extremes,
0472                              Intruders& intruders, SideStrategy const& strategy)
0473     {
0474         // Get one for the very first polygon, that is (for the moment) enough.
0475         // It is not guaranteed the "extreme" then, but for the current purpose
0476         // (point_on_surface) it can just be this point.
0477         if (boost::size(multi) >= 1)
0478         {
0479             return extreme_points
0480                 <
0481                     typename boost::range_value<MultiPolygon const>::type,
0482                     Dimension,
0483                     polygon_tag
0484                 >::apply(*boost::begin(multi), extremes, intruders, strategy);
0485         }
0486
0487         return false;
0488     }
0489 };
0490
0491 } // namespace dispatch
0492 #endif // DOXYGEN_NO_DISPATCH
0493
0494
0495 /*!
0496 \brief Returns extreme points (for Edge=1 in dimension 1, so the top,
0497        for Edge=0 in dimension 0, the right side)
0498 \note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that.
0499  */
0500 template
0501 <
0502     std::size_t Edge,
0503     typename Geometry,
0504     typename Extremes,
0505     typename Intruders,
0506     typename SideStrategy
0507 >
0508 inline bool extreme_points(Geometry const& geometry,
0509                            Extremes& extremes,
0510                            Intruders& intruders,
0511                            SideStrategy const& strategy)
0512 {
0513     concepts::check<Geometry const>();
0514
0515     // Extremes is not required to follow a geometry concept (but it should support an output iterator),
0516     // but its elements should fulfil the point-concept
0517     concepts::check<typename boost::range_value<Extremes>::type>();
0518
0519     // Intruders should contain collections which value type is point-concept
0520     // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept
0521     concepts::check
0522         <
0523             typename boost::range_value
0524                 <
0525                     typename boost::range_value<Intruders>::type
0526                 >::type
0527             const
0528         >();
0529
0530     return dispatch::extreme_points
0531             <
0532                 Geometry,
0533                 Edge
0534             >::apply(geometry, extremes, intruders, strategy);
0535 }
0536
0537
0538 template
0539 <
0540     std::size_t Edge,
0541     typename Geometry,
0542     typename Extremes,
0543     typename Intruders
0544 >
0545 inline bool extreme_points(Geometry const& geometry,
0546                            Extremes& extremes,
0547                            Intruders& intruders)
0548 {
0549     typedef typename strategy::side::services::default_strategy
0550             <
0551                 typename cs_tag<Geometry>::type
0552             >::type strategy_type;
0553
0554     return geometry::extreme_points<Edge>(geometry,extremes, intruders, strategy_type());
0555 }
0556
0557 }} // namespace boost::geometry
0558
0559
0560 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP