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