geometry/algorithms/point_on_surface.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 2014-2023.
0009 // Modifications copyright (c) 2014-2023 Oracle and/or its affiliates.
0010 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
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_POINT_ON_SURFACE_HPP
0018 #define BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
0019
0020 #include <boost/range/begin.hpp>
0021 #include <boost/range/end.hpp>
0022
0023 #include <boost/geometry/core/point_type.hpp>
0024
0025 #include <boost/geometry/geometries/concepts/check.hpp>
0026
0027 #include <boost/geometry/algorithms/detail/extreme_points.hpp>
0028 #include <boost/geometry/algorithms/detail/signed_size_type.hpp>
0029
0030 #include <boost/geometry/strategies/side.hpp>
0031
0032
0033 namespace boost { namespace geometry
0034 {
0035
0036
0037 #ifndef DOXYGEN_NO_DETAIL
0038 namespace detail { namespace point_on_surface
0039 {
0040
0041 template <typename CoordinateType, int Dimension>
0042 struct specific_coordinate_first
0043 {
0044     CoordinateType const m_value_to_be_first;
0045
0046     inline specific_coordinate_first(CoordinateType value_to_be_skipped)
0047         : m_value_to_be_first(value_to_be_skipped)
0048     {}
0049
0050     template <typename Point>
0051     inline bool operator()(Point const& lhs, Point const& rhs)
0052     {
0053         CoordinateType const lh = geometry::get<Dimension>(lhs);
0054         CoordinateType const rh = geometry::get<Dimension>(rhs);
0055
0056         // If both lhs and rhs equal m_value_to_be_first,
0057         // we should handle conform if lh < rh = FALSE
0058         // The first condition meets that, keep it first
0059         if (geometry::math::equals(rh, m_value_to_be_first))
0060         {
0061             // Handle conform lh < -INF, which is always false
0062             return false;
0063         }
0064
0065         if (geometry::math::equals(lh, m_value_to_be_first))
0066         {
0067             // Handle conform -INF < rh, which is always true
0068             return true;
0069         }
0070
0071         return lh < rh;
0072     }
0073 };
0074
0075 template <int Dimension, typename Collection, typename Value, typename Predicate>
0076 inline bool max_value(Collection const& collection, Value& the_max, Predicate const& predicate)
0077 {
0078     bool first = true;
0079     for (auto const& item : collection)
0080     {
0081         if (! item.empty())
0082         {
0083             Value the_value = geometry::get<Dimension>(*std::max_element(item.begin(), item.end(), predicate));
0084             if (first || the_value > the_max)
0085             {
0086                 the_max = the_value;
0087                 first = false;
0088             }
0089         }
0090     }
0091     return ! first;
0092 }
0093
0094
0095 template <int Dimension, typename Value>
0096 struct select_below
0097 {
0098     Value m_value;
0099     inline select_below(Value const& v)
0100         : m_value(v)
0101     {}
0102
0103     template <typename Intruder>
0104     inline bool operator()(Intruder const& intruder) const
0105     {
0106         if (intruder.empty())
0107         {
0108             return true;
0109         }
0110         Value max = geometry::get<Dimension>(*std::max_element(intruder.begin(), intruder.end(), detail::extreme_points::compare<Dimension>()));
0111         return geometry::math::equals(max, m_value) || max < m_value;
0112     }
0113 };
0114
0115 template <int Dimension, typename Value>
0116 struct adapt_base
0117 {
0118     Value m_value;
0119     inline adapt_base(Value const& v)
0120         : m_value(v)
0121     {}
0122
0123     template <typename Intruder>
0124     inline void operator()(Intruder& intruder) const
0125     {
0126         if (intruder.size() >= 3)
0127         {
0128             detail::extreme_points::move_along_vector<Dimension>(intruder, m_value);
0129         }
0130     }
0131 };
0132
0133 template <int Dimension, typename Value>
0134 struct min_of_intruder
0135 {
0136     template <typename Intruder>
0137     inline bool operator()(Intruder const& lhs, Intruder const& rhs) const
0138     {
0139         Value lhs_min = geometry::get<Dimension>(*std::min_element(lhs.begin(), lhs.end(), detail::extreme_points::compare<Dimension>()));
0140         Value rhs_min = geometry::get<Dimension>(*std::min_element(rhs.begin(), rhs.end(), detail::extreme_points::compare<Dimension>()));
0141         return lhs_min < rhs_min;
0142     }
0143 };
0144
0145
0146 template <typename Point, typename P>
0147 inline void calculate_average(Point& point, std::vector<P> const& points)
0148 {
0149     typedef typename geometry::coordinate_type<Point>::type coordinate_type;
0150
0151     coordinate_type x = 0;
0152     coordinate_type y = 0;
0153
0154     for (auto const& p : points)
0155     {
0156         x += geometry::get<0>(p);
0157         y += geometry::get<1>(p);
0158     }
0159
0160     signed_size_type const count = points.size();
0161     geometry::set<0>(point, x / count);
0162     geometry::set<1>(point, y / count);
0163 }
0164
0165
0166 template <int Dimension, typename Extremes, typename Intruders, typename CoordinateType>
0167 inline void replace_extremes_for_self_tangencies(Extremes& extremes, Intruders& intruders, CoordinateType const& max_intruder)
0168 {
0169     // This function handles self-tangencies.
0170     // Self-tangencies use, as usual, the major part of code...
0171
0172     //        ___ e
0173     //       /|\ \                                                            .
0174     //      / | \ \                                                           .
0175     //     /  |  \ \                                                          .
0176     //    /   |   \ \                                                         .
0177     //   / /\ |    \ \                                                        .
0178     //     i2    i1
0179
0180     // The picture above shows the extreme (outside, "e") and two intruders ("i1","i2")
0181     // Assume that "i1" is self-tangent with the extreme, in one point at the top
0182     // Now the "penultimate" value is searched, this is is the top of i2
0183     // Then everything including and below (this is "i2" here) is removed
0184     // Then the base of "i1" and of "e" is adapted to this penultimate value
0185     // It then looks like:
0186
0187     //      b ___ e
0188     //       /|\ \                                                            .
0189     //      / | \ \                                                           .
0190     //     /  |  \ \                                                          .
0191     //    /   |   \ \                                                         .
0192     //   a    c i1
0193
0194     // Then intruders (here "i1" but there may be more) are sorted from left to right
0195     // Finally points "a","b" and "c" (in this order) are selected as a new triangle.
0196     // This triangle will have a centroid which is inside (assumed that intruders left segment
0197     // is not equal to extremes left segment, but that polygon would be invalid)
0198
0199     // Find highest non-self tangent intrusion, if any
0200     CoordinateType penultimate_value;
0201     specific_coordinate_first<CoordinateType, Dimension> pu_compare(max_intruder);
0202     if (max_value<Dimension>(intruders, penultimate_value, pu_compare))
0203     {
0204         // Throw away all intrusions <= this value, and of the kept one set this as base.
0205         select_below<Dimension, CoordinateType> predicate(penultimate_value);
0206         intruders.erase
0207             (
0208                 std::remove_if(boost::begin(intruders), boost::end(intruders), predicate),
0209                 boost::end(intruders)
0210             );
0211         adapt_base<Dimension, CoordinateType> fe_predicate(penultimate_value);
0212         // Sort from left to right (or bottom to top if Dimension=0)
0213         std::for_each(boost::begin(intruders), boost::end(intruders), fe_predicate);
0214
0215         // Also adapt base of extremes
0216         detail::extreme_points::move_along_vector<Dimension>(extremes, penultimate_value);
0217     }
0218     // Then sort in 1-Dim. Take first to calc centroid.
0219     std::sort(boost::begin(intruders), boost::end(intruders), min_of_intruder<1 - Dimension, CoordinateType>());
0220
0221     Extremes triangle;
0222     triangle.reserve(3);
0223
0224     // 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)
0225     std::copy(extremes.begin(), extremes.begin() + 2, std::back_inserter(triangle));
0226     triangle.push_back(intruders.front().back());
0227
0228     // (alternatively we could use the last two points of extremes, and first point of last intruder...):
0229     //// ALTERNATIVE: std::copy(extremes.rbegin(), extremes.rbegin() + 2, std::back_inserter(triangle));
0230     //// ALTERNATIVE: triangle.push_back(intruders.back().front());
0231
0232     // Now replace extremes with this smaller subset, a triangle, such that centroid calculation will result in a point inside
0233     extremes = triangle;
0234 }
0235
0236 template <int Dimension, typename Geometry, typename Point, typename SideStrategy>
0237 inline bool calculate_point_on_surface(Geometry const& geometry, Point& point,
0238                                        SideStrategy const& strategy)
0239 {
0240     typedef typename geometry::point_type<Geometry>::type point_type;
0241     typedef typename geometry::coordinate_type<Geometry>::type coordinate_type;
0242     std::vector<point_type> extremes;
0243
0244     typedef std::vector<std::vector<point_type> > intruders_type;
0245     intruders_type intruders;
0246     geometry::extreme_points<Dimension>(geometry, extremes, intruders, strategy);
0247
0248     if (extremes.size() < 3)
0249     {
0250         return false;
0251     }
0252
0253     // If there are intruders, find the max.
0254     if (! intruders.empty())
0255     {
0256         coordinate_type max_intruder;
0257         detail::extreme_points::compare<Dimension> compare;
0258         if (max_value<Dimension>(intruders, max_intruder, compare))
0259         {
0260             coordinate_type max_extreme = geometry::get<Dimension>(*std::max_element(extremes.begin(), extremes.end(), detail::extreme_points::compare<Dimension>()));
0261             if (max_extreme > max_intruder)
0262             {
0263                 detail::extreme_points::move_along_vector<Dimension>(extremes, max_intruder);
0264             }
0265             else
0266             {
0267                 replace_extremes_for_self_tangencies<Dimension>(extremes, intruders, max_intruder);
0268             }
0269         }
0270     }
0271
0272     // Now calculate the average/centroid of the (possibly adapted) extremes
0273     calculate_average(point, extremes);
0274
0275     return true;
0276 }
0277
0278 }} // namespace detail::point_on_surface
0279 #endif // DOXYGEN_NO_DETAIL
0280
0281
0282 /*!
0283 \brief Assigns a Point guaranteed to lie on the surface of the Geometry
0284 \tparam Geometry geometry type. This also defines the type of the output point
0285 \param geometry Geometry to take point from
0286 \param point Point to assign
0287 \param strategy side strategy
0288  */
0289 template <typename Geometry, typename Point, typename SideStrategy>
0290 inline void point_on_surface(Geometry const& geometry, Point & point,
0291                              SideStrategy const& strategy)
0292 {
0293     concepts::check<Point>();
0294     concepts::check<Geometry const>();
0295
0296     // First try in Y-direction (which should always succeed for valid polygons)
0297     if (! detail::point_on_surface::calculate_point_on_surface<1>(geometry, point, strategy))
0298     {
0299         // For invalid polygons, we might try X-direction
0300         detail::point_on_surface::calculate_point_on_surface<0>(geometry, point, strategy);
0301     }
0302 }
0303
0304 /*!
0305 \brief Assigns a Point guaranteed to lie on the surface of the Geometry
0306 \tparam Geometry geometry type. This also defines the type of the output point
0307 \param geometry Geometry to take point from
0308 \param point Point to assign
0309  */
0310 template <typename Geometry, typename Point>
0311 inline void point_on_surface(Geometry const& geometry, Point & point)
0312 {
0313     typedef typename strategy::side::services::default_strategy
0314         <
0315             typename cs_tag<Geometry>::type
0316         >::type strategy_type;
0317
0318     point_on_surface(geometry, point, strategy_type());
0319 }
0320
0321
0322 /*!
0323 \brief Returns point guaranteed to lie on the surface of the Geometry
0324 \tparam Geometry geometry type. This also defines the type of the output point
0325 \param geometry Geometry to take point from
0326 \param strategy side strategy
0327 \return The Point guaranteed to lie on the surface of the Geometry
0328  */
0329 template<typename Geometry, typename SideStrategy>
0330 inline typename geometry::point_type<Geometry>::type
0331 return_point_on_surface(Geometry const& geometry, SideStrategy const& strategy)
0332 {
0333     typename geometry::point_type<Geometry>::type result;
0334     geometry::point_on_surface(geometry, result, strategy);
0335     return result;
0336 }
0337
0338 /*!
0339 \brief Returns point guaranteed to lie on the surface of the Geometry
0340 \tparam Geometry geometry type. This also defines the type of the output point
0341 \param geometry Geometry to take point from
0342 \return The Point guaranteed to lie on the surface of the Geometry
0343  */
0344 template<typename Geometry>
0345 inline typename geometry::point_type<Geometry>::type
0346 return_point_on_surface(Geometry const& geometry)
0347 {
0348     typename geometry::point_type<Geometry>::type result;
0349     geometry::point_on_surface(geometry, result);
0350     return result;
0351 }
0352
0353 }} // namespace boost::geometry
0354
0355
0356 #endif // BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP