ceph/src/boost/boost/geometry/algorithms/point_on_surface.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 2014-2020.
   9 // Modifications copyright (c) 2014-2020 Oracle and/or its affiliates.
  10 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  11
  12 // Use, modification and distribution is subject to the Boost Software License,
  13 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  14 // http://www.boost.org/LICENSE_1_0.txt)
  15
  16 #ifndef BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
  17 #define BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
  18
  19
  20 #include <cstddef>
  21 #include <numeric>
  22
  23 #include <boost/concept_check.hpp>
  24 #include <boost/range/begin.hpp>
  25 #include <boost/range/end.hpp>
  26
  27 #include <boost/geometry/core/point_type.hpp>
  28 #include <boost/geometry/core/ring_type.hpp>
  29
  30 #include <boost/geometry/geometries/concepts/check.hpp>
  31
  32 #include <boost/geometry/algorithms/detail/extreme_points.hpp>
  33 #include <boost/geometry/algorithms/detail/signed_size_type.hpp>
  34
  35 #include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp>
  36 #include <boost/geometry/strategies/side.hpp>
  37
  38
  39 namespace boost { namespace geometry
  40 {
  41
  42
  43 #ifndef DOXYGEN_NO_DETAIL
  44 namespace detail { namespace point_on_surface
  45 {
  46
  47 template <typename CoordinateType, int Dimension>
  48 struct specific_coordinate_first
  49 {
  50     CoordinateType const m_value_to_be_first;
  51
  52     inline specific_coordinate_first(CoordinateType value_to_be_skipped)
  53         : m_value_to_be_first(value_to_be_skipped)
  54     {}
  55
  56     template <typename Point>
  57     inline bool operator()(Point const& lhs, Point const& rhs)
  58     {
  59         CoordinateType const lh = geometry::get<Dimension>(lhs);
  60         CoordinateType const rh = geometry::get<Dimension>(rhs);
  61
  62         // If both lhs and rhs equal m_value_to_be_first,
  63         // we should handle conform if lh < rh = FALSE
  64         // The first condition meets that, keep it first
  65         if (geometry::math::equals(rh, m_value_to_be_first))
  66         {
  67             // Handle conform lh < -INF, which is always false
  68             return false;
  69         }
  70
  71         if (geometry::math::equals(lh, m_value_to_be_first))
  72         {
  73             // Handle conform -INF < rh, which is always true
  74             return true;
  75         }
  76
  77         return lh < rh;
  78     }
  79 };
  80
  81 template <int Dimension, typename Collection, typename Value, typename Predicate>
  82 inline bool max_value(Collection const& collection, Value& the_max, Predicate const& predicate)
  83 {
  84     bool first = true;
  85     for (typename Collection::const_iterator it = collection.begin(); it != collection.end(); ++it)
  86     {
  87         if (! it->empty())
  88         {
  89             Value the_value = geometry::get<Dimension>(*std::max_element(it->begin(), it->end(), predicate));
  90             if (first || the_value > the_max)
  91             {
  92                 the_max = the_value;
  93                 first = false;
  94             }
  95         }
  96     }
  97     return ! first;
  98 }
  99
 100
 101 template <int Dimension, typename Value>
 102 struct select_below
 103 {
 104     Value m_value;
 105     inline select_below(Value const& v)
 106         : m_value(v)
 107     {}
 108
 109     template <typename Intruder>
 110     inline bool operator()(Intruder const& intruder) const
 111     {
 112         if (intruder.empty())
 113         {
 114             return true;
 115         }
 116         Value max = geometry::get<Dimension>(*std::max_element(intruder.begin(), intruder.end(), detail::extreme_points::compare<Dimension>()));
 117         return geometry::math::equals(max, m_value) || max < m_value;
 118     }
 119 };
 120
 121 template <int Dimension, typename Value>
 122 struct adapt_base
 123 {
 124     Value m_value;
 125     inline adapt_base(Value const& v)
 126         : m_value(v)
 127     {}
 128
 129     template <typename Intruder>
 130     inline void operator()(Intruder& intruder) const
 131     {
 132         if (intruder.size() >= 3)
 133         {
 134             detail::extreme_points::move_along_vector<Dimension>(intruder, m_value);
 135         }
 136     }
 137 };
 138
 139 template <int Dimension, typename Value>
 140 struct min_of_intruder
 141 {
 142     template <typename Intruder>
 143     inline bool operator()(Intruder const& lhs, Intruder const& rhs) const
 144     {
 145         Value lhs_min = geometry::get<Dimension>(*std::min_element(lhs.begin(), lhs.end(), detail::extreme_points::compare<Dimension>()));
 146         Value rhs_min = geometry::get<Dimension>(*std::min_element(rhs.begin(), rhs.end(), detail::extreme_points::compare<Dimension>()));
 147         return lhs_min < rhs_min;
 148     }
 149 };
 150
 151
 152 template <typename Point, typename P>
 153 inline void calculate_average(Point& point, std::vector<P> const& points)
 154 {
 155     typedef typename geometry::coordinate_type<Point>::type coordinate_type;
 156     typedef typename std::vector<P>::const_iterator iterator_type;
 157
 158     coordinate_type x = 0;
 159     coordinate_type y = 0;
 160
 161     iterator_type end = points.end();
 162     for ( iterator_type it = points.begin() ; it != end ; ++it)
 163     {
 164         x += geometry::get<0>(*it);
 165         y += geometry::get<1>(*it);
 166     }
 167
 168     signed_size_type const count = points.size();
 169     geometry::set<0>(point, x / count);
 170     geometry::set<1>(point, y / count);
 171 }
 172
 173
 174 template <int Dimension, typename Extremes, typename Intruders, typename CoordinateType>
 175 inline void replace_extremes_for_self_tangencies(Extremes& extremes, Intruders& intruders, CoordinateType const& max_intruder)
 176 {
 177     // This function handles self-tangencies.
 178     // Self-tangencies use, as usual, the major part of code...
 179
 180     //        ___ e
 181     //       /|\ \                                                            .
 182     //      / | \ \                                                           .
 183     //     /  |  \ \                                                          .
 184     //    /   |   \ \                                                         .
 185     //   / /\ |    \ \                                                        .
 186     //     i2    i1
 187
 188     // The picture above shows the extreme (outside, "e") and two intruders ("i1","i2")
 189     // Assume that "i1" is self-tangent with the extreme, in one point at the top
 190     // Now the "penultimate" value is searched, this is is the top of i2
 191     // Then everything including and below (this is "i2" here) is removed
 192     // Then the base of "i1" and of "e" is adapted to this penultimate value
 193     // It then looks like:
 194
 195     //      b ___ e
 196     //       /|\ \                                                            .
 197     //      / | \ \                                                           .
 198     //     /  |  \ \                                                          .
 199     //    /   |   \ \                                                         .
 200     //   a    c i1
 201
 202     // Then intruders (here "i1" but there may be more) are sorted from left to right
 203     // Finally points "a","b" and "c" (in this order) are selected as a new triangle.
 204     // This triangle will have a centroid which is inside (assumed that intruders left segment
 205     // is not equal to extremes left segment, but that polygon would be invalid)
 206
 207     // Find highest non-self tangent intrusion, if any
 208     CoordinateType penultimate_value;
 209     specific_coordinate_first<CoordinateType, Dimension> pu_compare(max_intruder);
 210     if (max_value<Dimension>(intruders, penultimate_value, pu_compare))
 211     {
 212         // Throw away all intrusions <= this value, and of the kept one set this as base.
 213         select_below<Dimension, CoordinateType> predicate(penultimate_value);
 214         intruders.erase
 215             (
 216                 std::remove_if(boost::begin(intruders), boost::end(intruders), predicate),
 217                 boost::end(intruders)
 218             );
 219         adapt_base<Dimension, CoordinateType> fe_predicate(penultimate_value);
 220         // Sort from left to right (or bottom to top if Dimension=0)
 221         std::for_each(boost::begin(intruders), boost::end(intruders), fe_predicate);
 222
 223         // Also adapt base of extremes
 224         detail::extreme_points::move_along_vector<Dimension>(extremes, penultimate_value);
 225     }
 226     // Then sort in 1-Dim. Take first to calc centroid.
 227     std::sort(boost::begin(intruders), boost::end(intruders), min_of_intruder<1 - Dimension, CoordinateType>());
 228
 229     Extremes triangle;
 230     triangle.reserve(3);
 231
 232     // Make a triangle of first two points of extremes (the ramp, from left to right), and last point of first intruder (which goes from right to left)
 233     std::copy(extremes.begin(), extremes.begin() + 2, std::back_inserter(triangle));
 234     triangle.push_back(intruders.front().back());
 235
 236     // (alternatively we could use the last two points of extremes, and first point of last intruder...):
 237     //// ALTERNATIVE: std::copy(extremes.rbegin(), extremes.rbegin() + 2, std::back_inserter(triangle));
 238     //// ALTERNATIVE: triangle.push_back(intruders.back().front());
 239
 240     // Now replace extremes with this smaller subset, a triangle, such that centroid calculation will result in a point inside
 241     extremes = triangle;
 242 }
 243
 244 template <int Dimension, typename Geometry, typename Point, typename SideStrategy>
 245 inline bool calculate_point_on_surface(Geometry const& geometry, Point& point,
 246                                        SideStrategy const& strategy)
 247 {
 248     typedef typename geometry::point_type<Geometry>::type point_type;
 249     typedef typename geometry::coordinate_type<Geometry>::type coordinate_type;
 250     std::vector<point_type> extremes;
 251
 252     typedef std::vector<std::vector<point_type> > intruders_type;
 253     intruders_type intruders;
 254     geometry::extreme_points<Dimension>(geometry, extremes, intruders, strategy);
 255
 256     if (extremes.size() < 3)
 257     {
 258         return false;
 259     }
 260
 261     // If there are intruders, find the max.
 262     if (! intruders.empty())
 263     {
 264         coordinate_type max_intruder;
 265         detail::extreme_points::compare<Dimension> compare;
 266         if (max_value<Dimension>(intruders, max_intruder, compare))
 267         {
 268             coordinate_type max_extreme = geometry::get<Dimension>(*std::max_element(extremes.begin(), extremes.end(), detail::extreme_points::compare<Dimension>()));
 269             if (max_extreme > max_intruder)
 270             {
 271                 detail::extreme_points::move_along_vector<Dimension>(extremes, max_intruder);
 272             }
 273             else
 274             {
 275                 replace_extremes_for_self_tangencies<Dimension>(extremes, intruders, max_intruder);
 276             }
 277         }
 278     }
 279
 280     // Now calculate the average/centroid of the (possibly adapted) extremes
 281     calculate_average(point, extremes);
 282
 283     return true;
 284 }
 285
 286 }} // namespace detail::point_on_surface
 287 #endif // DOXYGEN_NO_DETAIL
 288
 289
 290 /*!
 291 \brief Assigns a Point guaranteed to lie on the surface of the Geometry
 292 \tparam Geometry geometry type. This also defines the type of the output point
 293 \param geometry Geometry to take point from
 294 \param point Point to assign
 295 \param strategy side strategy
 296  */
 297 template <typename Geometry, typename Point, typename SideStrategy>
 298 inline void point_on_surface(Geometry const& geometry, Point & point,
 299                              SideStrategy const& strategy)
 300 {
 301     concepts::check<Point>();
 302     concepts::check<Geometry const>();
 303
 304     // First try in Y-direction (which should always succeed for valid polygons)
 305     if (! detail::point_on_surface::calculate_point_on_surface<1>(geometry, point, strategy))
 306     {
 307         // For invalid polygons, we might try X-direction
 308         detail::point_on_surface::calculate_point_on_surface<0>(geometry, point, strategy);
 309     }
 310 }
 311
 312 /*!
 313 \brief Assigns a Point guaranteed to lie on the surface of the Geometry
 314 \tparam Geometry geometry type. This also defines the type of the output point
 315 \param geometry Geometry to take point from
 316 \param point Point to assign
 317  */
 318 template <typename Geometry, typename Point>
 319 inline void point_on_surface(Geometry const& geometry, Point & point)
 320 {
 321     typedef typename strategy::side::services::default_strategy
 322         <
 323             typename cs_tag<Geometry>::type
 324         >::type strategy_type;
 325
 326     point_on_surface(geometry, point, strategy_type());
 327 }
 328
 329
 330 /*!
 331 \brief Returns point guaranteed to lie on the surface of the Geometry
 332 \tparam Geometry geometry type. This also defines the type of the output point
 333 \param geometry Geometry to take point from
 334 \param strategy side strategy
 335 \return The Point guaranteed to lie on the surface of the Geometry
 336  */
 337 template<typename Geometry, typename SideStrategy>
 338 inline typename geometry::point_type<Geometry>::type
 339 return_point_on_surface(Geometry const& geometry, SideStrategy const& strategy)
 340 {
 341     typename geometry::point_type<Geometry>::type result;
 342     geometry::point_on_surface(geometry, result, strategy);
 343     return result;
 344 }
 345
 346 /*!
 347 \brief Returns point guaranteed to lie on the surface of the Geometry
 348 \tparam Geometry geometry type. This also defines the type of the output point
 349 \param geometry Geometry to take point from
 350 \return The Point guaranteed to lie on the surface of the Geometry
 351  */
 352 template<typename Geometry>
 353 inline typename geometry::point_type<Geometry>::type
 354 return_point_on_surface(Geometry const& geometry)
 355 {
 356     typename geometry::point_type<Geometry>::type result;
 357     geometry::point_on_surface(geometry, result);
 358     return result;
 359 }
 360
 361 }} // namespace boost::geometry
 362
 363
 364 #endif // BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP