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-2014 Adam Wulkiewicz, Lodz, Poland.
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_DETAIL_EXTREME_POINTS_HPP
Chris@102: #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP
Chris@102: 
Chris@102: 
Chris@102: #include <cstddef>
Chris@102: 
Chris@102: #include <boost/range.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/algorithms/detail/interior_iterator.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/core/cs.hpp>
Chris@102: #include <boost/geometry/core/point_type.hpp>
Chris@102: #include <boost/geometry/core/ring_type.hpp>
Chris@102: #include <boost/geometry/core/tags.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/geometries/concepts/check.hpp>
Chris@102: #include <boost/geometry/iterators/ever_circling_iterator.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/algorithms/detail/assign_box_corners.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/strategies/side.hpp>
Chris@102: 
Chris@102: #include <boost/geometry/util/math.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 extreme_points
Chris@102: {
Chris@102: 
Chris@102: template <std::size_t Dimension>
Chris@102: struct compare
Chris@102: {
Chris@102:     template <typename Point>
Chris@102:     inline bool operator()(Point const& lhs, Point const& rhs)
Chris@102:     {
Chris@102:         return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs);
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: template <std::size_t Dimension, typename PointType, typename CoordinateType>
Chris@102: inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value)
Chris@102: {
Chris@102:     // Moves a point along the vector (point, extreme) in the direction of the extreme  point
Chris@102:     // This adapts the possibly uneven legs of the triangle (or trapezium-like shape)
Chris@102:     //      _____extreme            _____
Chris@102:     //     /     \                 /     \                                    .
Chris@102:     //    /base   \         =>    /       \ point                             .
Chris@102:     //             \ point                                                    .
Chris@102:     //
Chris@102:     // For so-called intruders, it can be used to adapt both legs to the level of "base"
Chris@102:     // For the base, it can be used to adapt both legs to the level of the max-value of the intruders
Chris@102:     // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector
Chris@102: 
Chris@102:     CoordinateType const value = geometry::get<Dimension>(point);
Chris@102:     //if (geometry::math::equals(value, base_value))
Chris@102:     if (value >= base_value)
Chris@102:     {
Chris@102:         return;
Chris@102:     }
Chris@102: 
Chris@102:     PointType vector = point;
Chris@102:     subtract_point(vector, extreme);
Chris@102: 
Chris@102:     CoordinateType const diff = geometry::get<Dimension>(vector);
Chris@102: 
Chris@102:     // diff should never be zero
Chris@102:     // because of the way our triangle/trapezium is build.
Chris@102:     // We just return if it would be the case.
Chris@102:     if (geometry::math::equals(diff, 0))
Chris@102:     {
Chris@102:         return;
Chris@102:     }
Chris@102: 
Chris@102:     CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme);
Chris@102: 
Chris@102:     multiply_value(vector, base_diff);
Chris@102:     divide_value(vector, diff);
Chris@102: 
Chris@102:     // The real move:
Chris@102:     point = extreme;
Chris@102:     add_point(point, vector);
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: template <std::size_t Dimension, typename Range, typename CoordinateType>
Chris@102: inline void move_along_vector(Range& range, CoordinateType const& base_value)
Chris@102: {
Chris@102:     if (range.size() >= 3)
Chris@102:     {
Chris@102:         move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value);
Chris@102:         move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value);
Chris@102:     }
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: template<typename Ring, std::size_t Dimension>
Chris@102: struct extreme_points_on_ring
Chris@102: {
Chris@102: 
Chris@102:     typedef typename geometry::coordinate_type<Ring>::type coordinate_type;
Chris@102:     typedef typename boost::range_iterator<Ring const>::type range_iterator;
Chris@102:     typedef typename geometry::point_type<Ring>::type point_type;
Chris@102: 
Chris@102:     typedef typename geometry::strategy::side::services::default_strategy
Chris@102:         <
Chris@102:             typename geometry::cs_tag<point_type>::type
Chris@102:         >::type side_strategy;
Chris@102: 
Chris@102: 
Chris@102:     template <typename CirclingIterator, typename Points>
Chris@102:     static inline bool extend(CirclingIterator& it,
Chris@102:             std::size_t n,
Chris@102:             coordinate_type max_coordinate_value,
Chris@102:             Points& points, int direction)
Chris@102:     {
Chris@102:         std::size_t safe_index = 0;
Chris@102:         do
Chris@102:         {
Chris@102:             it += direction;
Chris@102: 
Chris@102:             points.push_back(*it);
Chris@102: 
Chris@102:             if (safe_index++ >= n)
Chris@102:             {
Chris@102:                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
Chris@102:                 return false;
Chris@102:             }
Chris@102:         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: 
Chris@102:     // Overload without adding to poinst
Chris@102:     template <typename CirclingIterator>
Chris@102:     static inline bool extend(CirclingIterator& it,
Chris@102:             std::size_t n,
Chris@102:             coordinate_type max_coordinate_value,
Chris@102:             int direction)
Chris@102:     {
Chris@102:         std::size_t safe_index = 0;
Chris@102:         do
Chris@102:         {
Chris@102:             it += direction;
Chris@102: 
Chris@102:             if (safe_index++ >= n)
Chris@102:             {
Chris@102:                 // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop)
Chris@102:                 return false;
Chris@102:             }
Chris@102:         } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value));
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: 
Chris@102:     template <typename CirclingIterator>
Chris@102:     static inline bool extent_both_sides(Ring const& ring,
Chris@102:             point_type extreme,
Chris@102:             CirclingIterator& left,
Chris@102:             CirclingIterator& right)
Chris@102:     {
Chris@102:         std::size_t const n = boost::size(ring);
Chris@102:         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
Chris@102: 
Chris@102:         if (! extend(left, n, max_coordinate_value, -1))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102:         if (! extend(right, n, max_coordinate_value, +1))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: 
Chris@102:     template <typename Collection, typename CirclingIterator>
Chris@102:     static inline bool collect(Ring const& ring,
Chris@102:             point_type extreme,
Chris@102:             Collection& points,
Chris@102:             CirclingIterator& left,
Chris@102:             CirclingIterator& right)
Chris@102:     {
Chris@102:         std::size_t const n = boost::size(ring);
Chris@102:         coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme);
Chris@102: 
Chris@102:         // Collects first left, which is reversed (if more than one point) then adds the top itself, then right
Chris@102:         if (! extend(left, n, max_coordinate_value, points, -1))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102:         std::reverse(points.begin(), points.end());
Chris@102:         points.push_back(extreme);
Chris@102:         if (! extend(right, n, max_coordinate_value, points, +1))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: 
Chris@102:     template <typename Extremes, typename Intruders, typename CirclingIterator>
Chris@102:     static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right,
Chris@102:             Extremes const& extremes,
Chris@102:             Intruders& intruders)
Chris@102:     {
Chris@102:         if (boost::size(extremes) < 3)
Chris@102:         {
Chris@102:             return;
Chris@102:         }
Chris@102:         coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>()));
Chris@102: 
Chris@102:         // Also select left/right (if Dimension=1)
Chris@102:         coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
Chris@102:         coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>()));
Chris@102: 
Chris@102:         std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too
Chris@102:         std::size_t const n = boost::size(ring);
Chris@102:         while (left != right && defensive_check_index < n)
Chris@102:         {
Chris@102:             coordinate_type const coordinate = geometry::get<Dimension>(*right);
Chris@102:             coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right);
Chris@102:             if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max)
Chris@102:             {
Chris@102:                 int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
Chris@102:                 int const first_side = side_strategy::apply(*right, extremes.front(), *(extremes.begin() + 1)) * factor;
Chris@102:                 int const last_side = side_strategy::apply(*right, *(extremes.rbegin() + 1), extremes.back()) * factor;
Chris@102: 
Chris@102:                 // If not lying left from any of the extemes side
Chris@102:                 if (first_side != 1 && last_side != 1)
Chris@102:                 {
Chris@102:                     //std::cout << "first " << first_side << " last " << last_side << std::endl;
Chris@102: 
Chris@102:                     // we start at this intrusion until it is handled, and don't affect our initial left iterator
Chris@102:                     CirclingIterator left_intrusion_it = right;
Chris@102:                     typename boost::range_value<Intruders>::type intruder;
Chris@102:                     collect(ring, *right, intruder, left_intrusion_it, right);
Chris@102: 
Chris@102:                     // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies
Chris@102:                     // (we might postpone this action, it is often not necessary. However it is not time-consuming)
Chris@102:                     move_along_vector<Dimension>(intruder, min_value);
Chris@102:                     intruders.push_back(intruder);
Chris@102:                     --right;
Chris@102:                 }
Chris@102:             }
Chris@102:             ++right;
Chris@102:             defensive_check_index++;
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline void get_intruders(Ring const& ring,
Chris@102:             Extremes const& extremes,
Chris@102:             Intruders& intruders)
Chris@102:     {
Chris@102:         std::size_t const n = boost::size(ring);
Chris@102:         if (n >= 3)
Chris@102:         {
Chris@102:             geometry::ever_circling_range_iterator<Ring const> left(ring);
Chris@102:             geometry::ever_circling_range_iterator<Ring const> right(ring);
Chris@102:             ++right;
Chris@102: 
Chris@102:             get_intruders(ring, left, right, extremes, intruders);
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     template <typename Iterator>
Chris@102:     static inline bool right_turn(Ring const& ring, Iterator it)
Chris@102:     {
Chris@102:         typename std::iterator_traits<Iterator>::difference_type const index
Chris@102:             = std::distance(boost::begin(ring), it);
Chris@102:         geometry::ever_circling_range_iterator<Ring const> left(ring);
Chris@102:         geometry::ever_circling_range_iterator<Ring const> right(ring);
Chris@102:         left += index;
Chris@102:         right += index;
Chris@102: 
Chris@102:         if (! extent_both_sides(ring, *it, left, right))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1;
Chris@102:         int const first_side = side_strategy::apply(*(right - 1), *right, *left) * factor;
Chris@102:         int const last_side = side_strategy::apply(*left, *(left + 1), *right) * factor;
Chris@102: 
Chris@102: //std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl;
Chris@102: 
Chris@102:         // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases)
Chris@102:         return first_side != 1 && last_side != 1;
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline bool apply(Ring const& ring, Extremes& extremes, Intruders& intruders)
Chris@102:     {
Chris@102:         std::size_t const n = boost::size(ring);
Chris@102:         if (n < 3)
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         // Get all maxima, usually one. In case of self-tangencies, or self-crossings,
Chris@102:         // the max might be is not valid. A valid max should make a right turn
Chris@102:         range_iterator max_it = boost::begin(ring);
Chris@102:         compare<Dimension> smaller;
Chris@102:         for (range_iterator it = max_it + 1; it != boost::end(ring); ++it)
Chris@102:         {
Chris@102:             if (smaller(*max_it, *it) && right_turn(ring, it))
Chris@102:             {
Chris@102:                 max_it = it;
Chris@102:             }
Chris@102:         }
Chris@102: 
Chris@102:         if (max_it == boost::end(ring))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         typename std::iterator_traits<range_iterator>::difference_type const
Chris@102:             index = std::distance(boost::begin(ring), max_it);
Chris@102: //std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl;
Chris@102: 
Chris@102:         geometry::ever_circling_range_iterator<Ring const> left(ring);
Chris@102:         geometry::ever_circling_range_iterator<Ring const> right(ring);
Chris@102:         left += index;
Chris@102:         right += index;
Chris@102: 
Chris@102:         // Collect all points (often 3) in a temporary vector
Chris@102:         std::vector<point_type> points;
Chris@102:         points.reserve(3);
Chris@102:         if (! collect(ring, *max_it, points, left, right))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102: //std::cout << "Built vector of " << points.size() << std::endl;
Chris@102: 
Chris@102:         coordinate_type const front_value = geometry::get<Dimension>(points.front());
Chris@102:         coordinate_type const back_value = geometry::get<Dimension>(points.back());
Chris@102:         coordinate_type const base_value = (std::max)(front_value, back_value);
Chris@102:         if (front_value < back_value)
Chris@102:         {
Chris@102:             move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value);
Chris@102:         }
Chris@102:         else
Chris@102:         {
Chris@102:             move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value);
Chris@102:         }
Chris@102: 
Chris@102:         std::copy(points.begin(), points.end(), std::back_inserter(extremes));
Chris@102: 
Chris@102:         get_intruders(ring, left, right, extremes, intruders);
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: }} // namespace detail::extreme_points
Chris@102: #endif // DOXYGEN_NO_DETAIL
Chris@102: 
Chris@102: 
Chris@102: #ifndef DOXYGEN_NO_DISPATCH
Chris@102: namespace dispatch
Chris@102: {
Chris@102: 
Chris@102: 
Chris@102: template
Chris@102: <
Chris@102:     typename Geometry,
Chris@102:     std::size_t Dimension,
Chris@102:     typename GeometryTag = typename tag<Geometry>::type
Chris@102: >
Chris@102: struct extreme_points
Chris@102: {};
Chris@102: 
Chris@102: 
Chris@102: template<typename Ring, std::size_t Dimension>
Chris@102: struct extreme_points<Ring, Dimension, ring_tag>
Chris@102:     : detail::extreme_points::extreme_points_on_ring<Ring, Dimension>
Chris@102: {};
Chris@102: 
Chris@102: 
Chris@102: template<typename Polygon, std::size_t Dimension>
Chris@102: struct extreme_points<Polygon, Dimension, polygon_tag>
Chris@102: {
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders)
Chris@102:     {
Chris@102:         typedef typename geometry::ring_type<Polygon>::type ring_type;
Chris@102:         typedef detail::extreme_points::extreme_points_on_ring
Chris@102:             <
Chris@102:                 ring_type, Dimension
Chris@102:             > ring_implementation;
Chris@102: 
Chris@102:         if (! ring_implementation::apply(geometry::exterior_ring(polygon), extremes, intruders))
Chris@102:         {
Chris@102:             return false;
Chris@102:         }
Chris@102: 
Chris@102:         // For a polygon, its interior rings can contain intruders
Chris@102:         typename interior_return_type<Polygon const>::type
Chris@102:             rings = interior_rings(polygon);
Chris@102:         for (typename detail::interior_iterator<Polygon const>::type
Chris@102:                 it = boost::begin(rings); it != boost::end(rings); ++it)
Chris@102:         {
Chris@102:             ring_implementation::get_intruders(*it, extremes,  intruders);
Chris@102:         }
Chris@102: 
Chris@102:         return true;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template<typename Box>
Chris@102: struct extreme_points<Box, 1, box_tag>
Chris@102: {
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline bool apply(Box const& box, Extremes& extremes, Intruders& )
Chris@102:     {
Chris@102:         extremes.resize(4);
Chris@102:         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
Chris@102:         // ll,ul,ur,lr, contains too exactly the right info
Chris@102:         return true;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template<typename Box>
Chris@102: struct extreme_points<Box, 0, box_tag>
Chris@102: {
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline bool apply(Box const& box, Extremes& extremes, Intruders& )
Chris@102:     {
Chris@102:         extremes.resize(4);
Chris@102:         geometry::detail::assign_box_corners_oriented<false>(box, extremes);
Chris@102:         // ll,ul,ur,lr, rotate one to start with UL and end with LL
Chris@102:         std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end());
Chris@102:         return true;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: template<typename MultiPolygon, std::size_t Dimension>
Chris@102: struct extreme_points<MultiPolygon, Dimension, multi_polygon_tag>
Chris@102: {
Chris@102:     template <typename Extremes, typename Intruders>
Chris@102:     static inline bool apply(MultiPolygon const& multi, Extremes& extremes, Intruders& intruders)
Chris@102:     {
Chris@102:         // Get one for the very first polygon, that is (for the moment) enough.
Chris@102:         // It is not guaranteed the "extreme" then, but for the current purpose
Chris@102:         // (point_on_surface) it can just be this point.
Chris@102:         if (boost::size(multi) >= 1)
Chris@102:         {
Chris@102:             return extreme_points
Chris@102:                 <
Chris@102:                     typename boost::range_value<MultiPolygon const>::type,
Chris@102:                     Dimension,
Chris@102:                     polygon_tag
Chris@102:                 >::apply(*boost::begin(multi), extremes, intruders);
Chris@102:         }
Chris@102: 
Chris@102:         return false;
Chris@102:     }
Chris@102: };
Chris@102: 
Chris@102: } // namespace dispatch
Chris@102: #endif // DOXYGEN_NO_DISPATCH
Chris@102: 
Chris@102: 
Chris@102: /*!
Chris@102: \brief Returns extreme points (for Edge=1 in dimension 1, so the top,
Chris@102:        for Edge=0 in dimension 0, the right side)
Chris@102: \note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that.
Chris@102:  */
Chris@102: template <std::size_t Edge, typename Geometry, typename Extremes, typename Intruders>
Chris@102: inline bool extreme_points(Geometry const& geometry, Extremes& extremes, Intruders& intruders)
Chris@102: {
Chris@102:     concept::check<Geometry const>();
Chris@102: 
Chris@102:     // Extremes is not required to follow a geometry concept (but it should support an output iterator),
Chris@102:     // but its elements should fulfil the point-concept
Chris@102:     concept::check<typename boost::range_value<Extremes>::type>();
Chris@102: 
Chris@102:     // Intruders should contain collections which value type is point-concept
Chris@102:     // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept
Chris@102:     concept::check
Chris@102:         <
Chris@102:             typename boost::range_value
Chris@102:                 <
Chris@102:                     typename boost::range_value<Intruders>::type
Chris@102:                 >::type
Chris@102:             const
Chris@102:         >();
Chris@102: 
Chris@102:     return dispatch::extreme_points<Geometry, Edge>::apply(geometry, extremes, intruders);
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: 
Chris@102: }} // namespace boost::geometry
Chris@102: 
Chris@102: 
Chris@102: #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP