Chris@102: // Boost.Geometry (aka GGL, Generic Geometry Library)
Chris@102: 
Chris@102: // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
Chris@102: // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
Chris@102: // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
Chris@102: // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.
Chris@102: 
Chris@102: // This file was modified by Oracle on 2014.
Chris@102: // Modifications copyright (c) 2014 Oracle and/or its affiliates.
Chris@102: 
Chris@102: // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
Chris@102: 
Chris@102: // Use, modification and distribution is subject to the Boost Software License,
Chris@102: // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
Chris@102: // http://www.boost.org/LICENSE_1_0.txt)
Chris@102: 
Chris@102: #ifndef BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
Chris@102: #define BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
Chris@102: 
Chris@102: 
Chris@102: #include <cstddef>
Chris@102: 
Chris@102: #include <numeric>
Chris@102: 
Chris@102: #include <boost/concept_check.hpp>
Chris@102: #include <boost/range.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/core/point_type.hpp>
Chris@102: #include <boost/geometry/core/ring_type.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/geometries/concepts/check.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/algorithms/detail/extreme_points.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp>
Chris@102: 
Chris@102: 
Chris@102: namespace boost { namespace geometry
Chris@102: {
Chris@102: 
Chris@102: 
Chris@102: #ifndef DOXYGEN_NO_DETAIL
Chris@102: namespace detail { namespace point_on_surface
Chris@102: {
Chris@102: 
Chris@102: template <typename CoordinateType, int Dimension>
Chris@102: struct specific_coordinate_first
Chris@102: {
Chris@102:     CoordinateType const m_value_to_be_first;
Chris@102: 
Chris@102:     inline specific_coordinate_first(CoordinateType value_to_be_skipped)
Chris@102:         : m_value_to_be_first(value_to_be_skipped)
Chris@102:     {}
Chris@102: 
Chris@102:     template <typename Point>
Chris@102:     inline bool operator()(Point const& lhs, Point const& rhs)
Chris@102:     {
Chris@102:         CoordinateType const lh = geometry::get<Dimension>(lhs);
Chris@102:         CoordinateType const rh = geometry::get<Dimension>(rhs);
Chris@102: 
Chris@102:         // If both lhs and rhs equal m_value_to_be_first,
Chris@102:         // we should handle conform if lh < rh = FALSE
Chris@102:         // The first condition meets that, keep it first
Chris@102:         if (geometry::math::equals(rh, m_value_to_be_first))
Chris@102:         {
Chris@102:             // Handle conform lh < -INF, which is always false
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         if (geometry::math::equals(lh, m_value_to_be_first))
Chris@102:         {
Chris@102:             // Handle conform -INF < rh, which is always true
Chris@102:             return true;
Chris@102:         }
Chris@102: 
Chris@102:         return lh < rh;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <int Dimension, typename Collection, typename Value, typename Predicate>
Chris@102: inline bool max_value(Collection const& collection, Value& the_max, Predicate const& predicate)
Chris@102: {
Chris@102:     bool first = true;
Chris@102:     for (typename Collection::const_iterator it = collection.begin(); it != collection.end(); ++it)
Chris@102:     {
Chris@102:         if (! it->empty())
Chris@102:         {
Chris@102:             Value the_value = geometry::get<Dimension>(*std::max_element(it->begin(), it->end(), predicate));
Chris@102:             if (first || the_value > the_max)
Chris@102:             {
Chris@102:                 the_max = the_value;
Chris@102:                 first = false;
Chris@102:             }
Chris@102:         }
Chris@102:     }
Chris@102:     return ! first;
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: template <int Dimension, typename Value>
Chris@102: struct select_below
Chris@102: {
Chris@102:     Value m_value;
Chris@102:     inline select_below(Value const& v)
Chris@102:         : m_value(v)
Chris@102:     {}
Chris@102: 
Chris@102:     template <typename Intruder>
Chris@102:     inline bool operator()(Intruder const& intruder) const
Chris@102:     {
Chris@102:         if (intruder.empty())
Chris@102:         {
Chris@102:             return true;
Chris@102:         }
Chris@102:         Value max = geometry::get<Dimension>(*std::max_element(intruder.begin(), intruder.end(), detail::extreme_points::compare<Dimension>()));
Chris@102:         return geometry::math::equals(max, m_value) || max < m_value;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <int Dimension, typename Value>
Chris@102: struct adapt_base
Chris@102: {
Chris@102:     Value m_value;
Chris@102:     inline adapt_base(Value const& v)
Chris@102:         : m_value(v)
Chris@102:     {}
Chris@102: 
Chris@102:     template <typename Intruder>
Chris@102:     inline void operator()(Intruder& intruder) const
Chris@102:     {
Chris@102:         if (intruder.size() >= 3)
Chris@102:         {
Chris@102:             detail::extreme_points::move_along_vector<Dimension>(intruder, m_value);
Chris@102:         }
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template <int Dimension, typename Value>
Chris@102: struct min_of_intruder
Chris@102: {
Chris@102:     template <typename Intruder>
Chris@102:     inline bool operator()(Intruder const& lhs, Intruder const& rhs) const
Chris@102:     {
Chris@102:         Value lhs_min = geometry::get<Dimension>(*std::min_element(lhs.begin(), lhs.end(), detail::extreme_points::compare<Dimension>()));
Chris@102:         Value rhs_min = geometry::get<Dimension>(*std::min_element(rhs.begin(), rhs.end(), detail::extreme_points::compare<Dimension>()));
Chris@102:         return lhs_min < rhs_min;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: template <typename Point, typename P>
Chris@102: inline void calculate_average(Point& point, std::vector<P> const& points)
Chris@102: {
Chris@102:     typedef typename geometry::coordinate_type<Point>::type coordinate_type;
Chris@102:     typedef typename std::vector<P>::const_iterator iterator_type;
Chris@102:     typedef typename std::vector<P>::size_type size_type;
Chris@102: 
Chris@102:     coordinate_type x = 0;
Chris@102:     coordinate_type y = 0;
Chris@102: 
Chris@102:     iterator_type end = points.end();
Chris@102:     for ( iterator_type it = points.begin() ; it != end ; ++it)
Chris@102:     {
Chris@102:         x += geometry::get<0>(*it);
Chris@102:         y += geometry::get<1>(*it);
Chris@102:     }
Chris@102: 
Chris@102:     size_type const count = points.size();
Chris@102:     geometry::set<0>(point, x / count);
Chris@102:     geometry::set<1>(point, y / count);
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: template <int Dimension, typename Extremes, typename Intruders, typename CoordinateType>
Chris@102: inline void replace_extremes_for_self_tangencies(Extremes& extremes, Intruders& intruders, CoordinateType const& max_intruder)
Chris@102: {
Chris@102:     // This function handles self-tangencies.
Chris@102:     // Self-tangencies use, as usual, the major part of code...
Chris@102: 
Chris@102:     //        ___ e
Chris@102:     //       /|\ \                                                            .
Chris@102:     //      / | \ \                                                           .
Chris@102:     //     /  |  \ \                                                          .
Chris@102:     //    /   |   \ \                                                         .
Chris@102:     //   / /\ |    \ \                                                        .
Chris@102:     //     i2    i1
Chris@102: 
Chris@102:     // The picture above shows the extreme (outside, "e") and two intruders ("i1","i2")
Chris@102:     // Assume that "i1" is self-tangent with the extreme, in one point at the top
Chris@102:     // Now the "penultimate" value is searched, this is is the top of i2
Chris@102:     // Then everything including and below (this is "i2" here) is removed
Chris@102:     // Then the base of "i1" and of "e" is adapted to this penultimate value
Chris@102:     // It then looks like:
Chris@102: 
Chris@102:     //      b ___ e
Chris@102:     //       /|\ \                                                            .
Chris@102:     //      / | \ \                                                           .
Chris@102:     //     /  |  \ \                                                          .
Chris@102:     //    /   |   \ \                                                         .
Chris@102:     //   a    c i1
Chris@102: 
Chris@102:     // Then intruders (here "i1" but there may be more) are sorted from left to right
Chris@102:     // Finally points "a","b" and "c" (in this order) are selected as a new triangle.
Chris@102:     // This triangle will have a centroid which is inside (assumed that intruders left segment
Chris@102:     // is not equal to extremes left segment, but that polygon would be invalid)
Chris@102: 
Chris@102:     // Find highest non-self tangent intrusion, if any
Chris@102:     CoordinateType penultimate_value;
Chris@102:     specific_coordinate_first<CoordinateType, Dimension> pu_compare(max_intruder);
Chris@102:     if (max_value<Dimension>(intruders, penultimate_value, pu_compare))
Chris@102:     {
Chris@102:         // Throw away all intrusions <= this value, and of the kept one set this as base.
Chris@102:         select_below<Dimension, CoordinateType> predicate(penultimate_value);
Chris@102:         intruders.erase
Chris@102:             (
Chris@102:                 std::remove_if(boost::begin(intruders), boost::end(intruders), predicate),
Chris@102:                 boost::end(intruders)
Chris@102:             );
Chris@102:         adapt_base<Dimension, CoordinateType> fe_predicate(penultimate_value);
Chris@102:         // Sort from left to right (or bottom to top if Dimension=0)
Chris@102:         std::for_each(boost::begin(intruders), boost::end(intruders), fe_predicate);
Chris@102: 
Chris@102:         // Also adapt base of extremes
Chris@102:         detail::extreme_points::move_along_vector<Dimension>(extremes, penultimate_value);
Chris@102:     }
Chris@102:     // Then sort in 1-Dim. Take first to calc centroid.
Chris@102:     std::sort(boost::begin(intruders), boost::end(intruders), min_of_intruder<1 - Dimension, CoordinateType>());
Chris@102: 
Chris@102:     Extremes triangle;
Chris@102:     triangle.reserve(3);
Chris@102: 
Chris@102:     // 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)
Chris@102:     std::copy(extremes.begin(), extremes.begin() + 2, std::back_inserter(triangle));
Chris@102:     triangle.push_back(intruders.front().back());
Chris@102: 
Chris@102:     // (alternatively we could use the last two points of extremes, and first point of last intruder...):
Chris@102:     //// ALTERNATIVE: std::copy(extremes.rbegin(), extremes.rbegin() + 2, std::back_inserter(triangle));
Chris@102:     //// ALTERNATIVE: triangle.push_back(intruders.back().front());
Chris@102: 
Chris@102:     // Now replace extremes with this smaller subset, a triangle, such that centroid calculation will result in a point inside
Chris@102:     extremes = triangle;
Chris@102: }
Chris@102: 
Chris@102: template <int Dimension, typename Geometry, typename Point>
Chris@102: inline bool calculate_point_on_surface(Geometry const& geometry, Point& point)
Chris@102: {
Chris@102:     typedef typename geometry::point_type<Geometry>::type point_type;
Chris@102:     typedef typename geometry::coordinate_type<Geometry>::type coordinate_type;
Chris@102:     std::vector<point_type> extremes;
Chris@102: 
Chris@102:     typedef std::vector<std::vector<point_type> > intruders_type;
Chris@102:     intruders_type intruders;
Chris@102:     geometry::extreme_points<Dimension>(geometry, extremes, intruders);
Chris@102: 
Chris@102:     if (extremes.size() < 3)
Chris@102:     {
Chris@102:         return false;
Chris@102:     }
Chris@102: 
Chris@102:     // If there are intruders, find the max.
Chris@102:     if (! intruders.empty())
Chris@102:     {
Chris@102:         coordinate_type max_intruder;
Chris@102:         detail::extreme_points::compare<Dimension> compare;
Chris@102:         if (max_value<Dimension>(intruders, max_intruder, compare))
Chris@102:         {
Chris@102:             coordinate_type max_extreme = geometry::get<Dimension>(*std::max_element(extremes.begin(), extremes.end(), detail::extreme_points::compare<Dimension>()));
Chris@102:             if (max_extreme > max_intruder)
Chris@102:             {
Chris@102:                 detail::extreme_points::move_along_vector<Dimension>(extremes, max_intruder);
Chris@102:             }
Chris@102:             else
Chris@102:             {
Chris@102:                 replace_extremes_for_self_tangencies<Dimension>(extremes, intruders, max_intruder);
Chris@102:             }
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     // Now calculate the average/centroid of the (possibly adapted) extremes
Chris@102:     calculate_average(point, extremes);
Chris@102: 
Chris@102:     return true;
Chris@102: }
Chris@102: 
Chris@102: }} // namespace detail::point_on_surface
Chris@102: #endif // DOXYGEN_NO_DETAIL
Chris@102: 
Chris@102: 
Chris@102: /*!
Chris@102: \brief Assigns a Point guaranteed to lie on the surface of the Geometry
Chris@102: \tparam Geometry geometry type. This also defines the type of the output point
Chris@102: \param geometry Geometry to take point from
Chris@102: \param point Point to assign
Chris@102:  */
Chris@102: template <typename Geometry, typename Point>
Chris@102: inline void point_on_surface(Geometry const& geometry, Point & point)
Chris@102: {
Chris@102:     concept::check<Point>();
Chris@102:     concept::check<Geometry const>();
Chris@102: 
Chris@102:     // First try in Y-direction (which should always succeed for valid polygons)
Chris@102:     if (! detail::point_on_surface::calculate_point_on_surface<1>(geometry, point))
Chris@102:     {
Chris@102:         // For invalid polygons, we might try X-direction
Chris@102:         detail::point_on_surface::calculate_point_on_surface<0>(geometry, point);
Chris@102:     }
Chris@102: }
Chris@102: 
Chris@102: /*!
Chris@102: \brief Returns point guaranteed to lie on the surface of the Geometry
Chris@102: \tparam Geometry geometry type. This also defines the type of the output point
Chris@102: \param geometry Geometry to take point from
Chris@102: \return The Point guaranteed to lie on the surface of the Geometry
Chris@102:  */
Chris@102: template<typename Geometry>
Chris@102: inline typename geometry::point_type<Geometry>::type
Chris@102: return_point_on_surface(Geometry const& geometry)
Chris@102: {
Chris@102:     typename geometry::point_type<Geometry>::type result;
Chris@102:     geometry::point_on_surface(geometry, result);
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: }} // namespace boost::geometry
Chris@102: 
Chris@102: 
Chris@102: #endif // BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP