ceph/src/boost/boost/geometry/algorithms/detail/extreme_points.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
   4 // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
   5 // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
   6 // Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
   7
   8 // This file was modified by Oracle on 2017-2020.
   9 // Modifications copyright (c) 2017-2020 Oracle and/or its affiliates.
  10
  11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  12
  13 // Use, modification and distribution is subject to the Boost Software License,
  14 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  15 // http://www.boost.org/LICENSE_1_0.txt)
  16
  17 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
  18 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
  19
  20
  21 #include <cstddef>
  22
  23 #include <boost/range/begin.hpp>
  24 #include <boost/range/end.hpp>
  25 #include <boost/range/size.hpp>
  26 #include <boost/range/value_type.hpp>
  27
  28 #include <boost/geometry/algorithms/detail/interior_iterator.hpp>
  29
  30 #include <boost/geometry/core/cs.hpp>
  31 #include <boost/geometry/core/point_order.hpp>
  32 #include <boost/geometry/core/point_type.hpp>
  33 #include <boost/geometry/core/ring_type.hpp>
  34 #include <boost/geometry/core/tags.hpp>
  35
  36 #include <boost/geometry/geometries/concepts/check.hpp>
  37 #include <boost/geometry/iterators/ever_circling_iterator.hpp>
  38
  39 #include <boost/geometry/algorithms/detail/assign_box_corners.hpp>
  40
  41 #include <boost/geometry/strategies/side.hpp>
  42
  43 #include <boost/geometry/util/math.hpp>
  44
  45
  46 namespace boost { namespace geometry
  47 {
  48
  49
  50 #ifndef DOXYGEN_NO_DETAIL
  51 namespace detail { namespace extreme_points
  52 {
  53
  54 template <std::size_t Dimension>
  55 struct compare
  56 {
  57     template <typename Point>
  58     inline bool operator()(Point const& lhs, Point const& rhs)
  59     {
  60         return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs);
  61     }
  62 };
  63
  64
  65 template <std::size_t Dimension, typename PointType, typename CoordinateType>
  66 inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value)
  67 {
  68     // Moves a point along the vector (point, extreme) in the direction of the extreme  point
  69     // This adapts the possibly uneven legs of the triangle (or trapezium-like shape)
  70     //      _____extreme            _____
  71     //     /     \                 /     \                                    .
  72     //    /base   \         =>    /       \ point                             .
  73     //             \ point                                                    .
  74     //
  75     // For so-called intruders, it can be used to adapt both legs to the level of "base"
  76     // For the base, it can be used to adapt both legs to the level of the max-value of the intruders
  77     // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector
  78
  79     CoordinateType const value = geometry::get<Dimension>(point);
  80     //if (geometry::math::equals(value, base_value))
  81     if (value >= base_value)
  82     {
  83         return;
  84     }
  85
  86     PointType vector = point;
  87     subtract_point(vector, extreme);
  88
  89     CoordinateType const diff = geometry::get<Dimension>(vector);
  90
  91     // diff should never be zero
  92     // because of the way our triangle/trapezium is build.
  93     // We just return if it would be the case.
  94     if (geometry::math::equals(diff, 0))
  95     {
  96         return;
  97     }
  98
  99     CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme);
 100
 101     multiply_value(vector, base_diff);
 102     divide_value(vector, diff);
 103
 104     // The real move:
 105     point = extreme;
 106     add_point(point, vector);
 107 }
 108
 109
 110 template <std::size_t Dimension, typename Range, typename CoordinateType>
 111 inline void move_along_vector(Range& range, CoordinateType const& base_value)
 112 {
 113     if (range.size() >= 3)
 114     {
 115         move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value);
 116         move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value);
 117     }
 118 }
 119
 120
 121 template<typename Ring, std::size_t Dimension>
 122 struct extreme_points_on_ring
 123 {
 124
 125     typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
 126     typedef typename boost::range_iterator<Ring const>::type range_iterator;
 127     typedef typename geometry::point_type<Ring>::type point_type;
 128
 129     template <typename CirclingIterator, typename Points>
 130     static inline bool extend(CirclingIterator& it,
 131             std::size_t n,
 132             coordinate_type max_coordinate_value,
 133             Points& points, int direction)
 134     {
 135         std::size_t safe_index = 0;
 136         do
 137         {
 138             it += direction;
 139
 140             points.push_back(*it);
 141
 142             if (safe_index++ >= n)
 143             {
 144                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
 145                 return false;
 146             }
 147         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
 148
 149         return true;
 150     }
 151
 152     // Overload without adding to poinst
 153     template <typename CirclingIterator>
 154     static inline bool extend(CirclingIterator& it,
 155             std::size_t n,
 156             coordinate_type max_coordinate_value,
 157             int direction)
 158     {
 159         std::size_t safe_index = 0;
 160         do
 161         {
 162             it += direction;
 163
 164             if (safe_index++ >= n)
 165             {
 166                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
 167                 return false;
 168             }
 169         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
 170
 171         return true;
 172     }
 173
 174     template <typename CirclingIterator>
 175     static inline bool extent_both_sides(Ring const& ring,
 176             point_type extreme,
 177             CirclingIterator& left,
 178             CirclingIterator& right)
 179     {
 180         std::size_t const n = boost::size(ring);
 181         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
 182
 183         if (! extend(left, n, max_coordinate_value, -1))
 184         {
 185             return false;
 186         }
 187         if (! extend(right, n, max_coordinate_value, +1))
 188         {
 189             return false;
 190         }
 191
 192         return true;
 193     }
 194
 195     template <typename Collection, typename CirclingIterator>
 196     static inline bool collect(Ring const& ring,
 197             point_type extreme,
 198             Collection& points,
 199             CirclingIterator& left,
 200             CirclingIterator& right)
 201     {
 202         std::size_t const n = boost::size(ring);
 203         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
 204
 205         // Collects first left, which is reversed (if more than one point) then adds the top itself, then right
 206         if (! extend(left, n, max_coordinate_value, points, -1))
 207         {
 208             return false;
 209         }
 210         std::reverse(points.begin(), points.end());
 211         points.push_back(extreme);
 212         if (! extend(right, n, max_coordinate_value, points, +1))
 213         {
 214             return false;
 215         }
 216
 217         return true;
 218     }
 219
 220     template <typename Extremes, typename Intruders, typename CirclingIterator, typename SideStrategy>
 221     static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right,
 222             Extremes const& extremes,
 223             Intruders& intruders,
 224             SideStrategy const& strategy)
 225     {
 226         if (boost::size(extremes) < 3)
 227         {
 228             return;
 229         }
 230         coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>()));
 231
 232         // Also select left/right (if Dimension=1)
 233         coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
 234         coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
 235
 236         std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too
 237         std::size_t const n = boost::size(ring);
 238         while (left != right && defensive_check_index < n)
 239         {
 240             coordinate_type const coordinate = geometry::get<Dimension>(*right);
 241             coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right);
 242             if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max)
 243             {
 244                 int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
 245                 int const first_side = strategy.apply(*right, extremes.front(), *(extremes.begin() + 1)) * factor;
 246                 int const last_side = strategy.apply(*right, *(extremes.rbegin() + 1), extremes.back()) * factor;
 247
 248                 // If not lying left from any of the extemes side
 249                 if (first_side != 1 && last_side != 1)
 250                 {
 251                     //std::cout << "first " << first_side << " last " << last_side << std::endl;
 252
 253                     // we start at this intrusion until it is handled, and don't affect our initial left iterator
 254                     CirclingIterator left_intrusion_it = right;
 255                     typename boost::range_value<Intruders>::type intruder;
 256                     collect(ring, *right, intruder, left_intrusion_it, right);
 257
 258                     // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies
 259                     // (we might postpone this action, it is often not necessary. However it is not time-consuming)
 260                     move_along_vector<Dimension>(intruder, min_value);
 261                     intruders.push_back(intruder);
 262                     --right;
 263                 }
 264             }
 265             ++right;
 266             defensive_check_index++;
 267         }
 268     }
 269
 270     template <typename Extremes, typename Intruders, typename SideStrategy>
 271     static inline void get_intruders(Ring const& ring,
 272             Extremes const& extremes,
 273             Intruders& intruders,
 274             SideStrategy const& strategy)
 275     {
 276         std::size_t const n = boost::size(ring);
 277         if (n >= 3)
 278         {
 279             geometry::ever_circling_range_iterator<Ring const> left(ring);
 280             geometry::ever_circling_range_iterator<Ring const> right(ring);
 281             ++right;
 282
 283             get_intruders(ring, left, right, extremes, intruders, strategy);
 284         }
 285     }
 286
 287     template <typename Iterator, typename SideStrategy>
 288     static inline bool right_turn(Ring const& ring, Iterator it, SideStrategy const& strategy)
 289     {
 290         typename std::iterator_traits<Iterator>::difference_type const index
 291             = std::distance(boost::begin(ring), it);
 292         geometry::ever_circling_range_iterator<Ring const> left(ring);
 293         geometry::ever_circling_range_iterator<Ring const> right(ring);
 294         left += index;
 295         right += index;
 296
 297         if (! extent_both_sides(ring, *it, left, right))
 298         {
 299             return false;
 300         }
 301
 302         int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
 303         int const first_side = strategy.apply(*(right - 1), *right, *left) * factor;
 304         int const last_side = strategy.apply(*left, *(left + 1), *right) * factor;
 305
 306 //std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl;
 307
 308         // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases)
 309         return first_side != 1 && last_side != 1;
 310     }
 311
 312
 313     // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring
 314     template <typename Extremes, typename Intruders, typename SideStrategy>
 315     static inline bool apply(Ring const& ring,
 316                              Extremes& extremes,
 317                              Intruders& intruders,
 318                              SideStrategy const& strategy)
 319     {
 320         std::size_t const n = boost::size(ring);
 321         if (n < 3)
 322         {
 323             return false;
 324         }
 325
 326         // Get all maxima, usually one. In case of self-tangencies, or self-crossings,
 327         // the max might be is not valid. A valid max should make a right turn
 328         range_iterator max_it = boost::begin(ring);
 329         compare<Dimension> smaller;
 330         for (range_iterator it = max_it + 1; it != boost::end(ring); ++it)
 331         {
 332             if (smaller(*max_it, *it) && right_turn(ring, it, strategy))
 333             {
 334                 max_it = it;
 335             }
 336         }
 337
 338         if (max_it == boost::end(ring))
 339         {
 340             return false;
 341         }
 342
 343         typename std::iterator_traits<range_iterator>::difference_type const
 344             index = std::distance(boost::begin(ring), max_it);
 345 //std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl;
 346
 347         geometry::ever_circling_range_iterator<Ring const> left(ring);
 348         geometry::ever_circling_range_iterator<Ring const> right(ring);
 349         left += index;
 350         right += index;
 351
 352         // Collect all points (often 3) in a temporary vector
 353         std::vector<point_type> points;
 354         points.reserve(3);
 355         if (! collect(ring, *max_it, points, left, right))
 356         {
 357             return false;
 358         }
 359
 360 //std::cout << "Built vector of " << points.size() << std::endl;
 361
 362         coordinate_type const front_value = geometry::get<Dimension>(points.front());
 363         coordinate_type const back_value = geometry::get<Dimension>(points.back());
 364         coordinate_type const base_value = (std::max)(front_value, back_value);
 365         if (front_value < back_value)
 366         {
 367             move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value);
 368         }
 369         else
 370         {
 371             move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value);
 372         }
 373
 374         std::copy(points.begin(), points.end(), std::back_inserter(extremes));
 375
 376         get_intruders(ring, left, right, extremes, intruders, strategy);
 377
 378         return true;
 379     }
 380 };
 381
 382
 383
 384
 385
 386 }} // namespace detail::extreme_points
 387 #endif // DOXYGEN_NO_DETAIL
 388
 389
 390 #ifndef DOXYGEN_NO_DISPATCH
 391 namespace dispatch
 392 {
 393
 394
 395 template
 396 <
 397     typename Geometry,
 398     std::size_t Dimension,
 399     typename GeometryTag = typename tag<Geometry>::type
 400 >
 401 struct extreme_points
 402 {};
 403
 404
 405 template<typename Ring, std::size_t Dimension>
 406 struct extreme_points<Ring, Dimension, ring_tag>
 407     : detail::extreme_points::extreme_points_on_ring<Ring, Dimension>
 408 {};
 409
 410
 411 template<typename Polygon, std::size_t Dimension>
 412 struct extreme_points<Polygon, Dimension, polygon_tag>
 413 {
 414     template <typename Extremes, typename Intruders, typename SideStrategy>
 415     static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders,
 416                              SideStrategy const& strategy)
 417     {
 418         typedef typename geometry::ring_type<Polygon>::type ring_type;
 419         typedef detail::extreme_points::extreme_points_on_ring
 420             <
 421                 ring_type, Dimension
 422             > ring_implementation;
 423
 424         if (! ring_implementation::apply(geometry::exterior_ring(polygon),
 425                                          extremes, intruders, strategy))
 426         {
 427             return false;
 428         }
 429
 430         // For a polygon, its interior rings can contain intruders
 431         typename interior_return_type<Polygon const>::type
 432             rings = interior_rings(polygon);
 433         for (typename detail::interior_iterator<Polygon const>::type
 434                 it = boost::begin(rings); it != boost::end(rings); ++it)
 435         {
 436             ring_implementation::get_intruders(*it, extremes,  intruders, strategy);
 437         }
 438
 439         return true;
 440     }
 441 };
 442
 443 template<typename Box>
 444 struct extreme_points<Box, 1, box_tag>
 445 {
 446     template <typename Extremes, typename Intruders, typename SideStrategy>
 447     static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
 448                              SideStrategy const& )
 449     {
 450         extremes.resize(4);
 451         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
 452         // ll,ul,ur,lr, contains too exactly the right info
 453         return true;
 454     }
 455 };
 456
 457 template<typename Box>
 458 struct extreme_points<Box, 0, box_tag>
 459 {
 460     template <typename Extremes, typename Intruders, typename SideStrategy>
 461     static inline bool apply(Box const& box, Extremes& extremes, Intruders& ,
 462                              SideStrategy const& )
 463     {
 464         extremes.resize(4);
 465         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
 466         // ll,ul,ur,lr, rotate one to start with UL and end with LL
 467         std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end());
 468         return true;
 469     }
 470 };
 471
 472 template<typename MultiPolygon, std::size_t Dimension>
 473 struct extreme_points<MultiPolygon, Dimension, multi_polygon_tag>
 474 {
 475     template <typename Extremes, typename Intruders, typename SideStrategy>
 476     static inline bool apply(MultiPolygon const& multi, Extremes& extremes,
 477                              Intruders& intruders, SideStrategy const& strategy)
 478     {
 479         // Get one for the very first polygon, that is (for the moment) enough.
 480         // It is not guaranteed the "extreme" then, but for the current purpose
 481         // (point_on_surface) it can just be this point.
 482         if (boost::size(multi) >= 1)
 483         {
 484             return extreme_points
 485                 <
 486                     typename boost::range_value<MultiPolygon const>::type,
 487                     Dimension,
 488                     polygon_tag
 489                 >::apply(*boost::begin(multi), extremes, intruders, strategy);
 490         }
 491
 492         return false;
 493     }
 494 };
 495
 496 } // namespace dispatch
 497 #endif // DOXYGEN_NO_DISPATCH
 498
 499
 500 /*!
 501 \brief Returns extreme points (for Edge=1 in dimension 1, so the top,
 502        for Edge=0 in dimension 0, the right side)
 503 \note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that.
 504  */
 505 template
 506 <
 507     std::size_t Edge,
 508     typename Geometry,
 509     typename Extremes,
 510     typename Intruders,
 511     typename SideStrategy
 512 >
 513 inline bool extreme_points(Geometry const& geometry,
 514                            Extremes& extremes,
 515                            Intruders& intruders,
 516                            SideStrategy const& strategy)
 517 {
 518     concepts::check<Geometry const>();
 519
 520     // Extremes is not required to follow a geometry concept (but it should support an output iterator),
 521     // but its elements should fulfil the point-concept
 522     concepts::check<typename boost::range_value<Extremes>::type>();
 523
 524     // Intruders should contain collections which value type is point-concept
 525     // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept
 526     concepts::check
 527         <
 528             typename boost::range_value
 529                 <
 530                     typename boost::range_value<Intruders>::type
 531                 >::type
 532             const
 533         >();
 534
 535     return dispatch::extreme_points
 536             <
 537                 Geometry,
 538                 Edge
 539             >::apply(geometry, extremes, intruders, strategy);
 540 }
 541
 542
 543 template
 544 <
 545     std::size_t Edge,
 546     typename Geometry,
 547     typename Extremes,
 548     typename Intruders
 549 >
 550 inline bool extreme_points(Geometry const& geometry,
 551                            Extremes& extremes,
 552                            Intruders& intruders)
 553 {
 554     typedef typename strategy::side::services::default_strategy
 555             <
 556                 typename cs_tag<Geometry>::type
 557             >::type strategy_type;
 558
 559     return geometry::extreme_points<Edge>(geometry,extremes, intruders, strategy_type());
 560 }
 561
 562 }} // namespace boost::geometry
 563
 564
 565 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP