1 // Boost.Geometry (aka GGL, Generic Geometry Library)
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 Adam Wulkiewicz, Lodz, Poland.
8 // This file was modified by Oracle on 2014.
9 // Modifications copyright (c) 2014 Oracle and/or its affiliates.
11 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
13 // Use, modification and distribution is subject to the Boost Software License,
14 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
15 // http://www.boost.org/LICENSE_1_0.txt)
17 #ifndef BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
18 #define BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP
25 #include <boost/concept_check.hpp>
26 #include <boost/range.hpp>
28 #include <boost/geometry/core/point_type.hpp>
29 #include <boost/geometry/core/ring_type.hpp>
31 #include <boost/geometry/geometries/concepts/check.hpp>
33 #include <boost/geometry/algorithms/detail/extreme_points.hpp>
35 #include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp>
38 namespace boost { namespace geometry
42 #ifndef DOXYGEN_NO_DETAIL
43 namespace detail { namespace point_on_surface
46 template <typename CoordinateType, int Dimension>
47 struct specific_coordinate_first
49 CoordinateType const m_value_to_be_first;
51 inline specific_coordinate_first(CoordinateType value_to_be_skipped)
52 : m_value_to_be_first(value_to_be_skipped)
55 template <typename Point>
56 inline bool operator()(Point const& lhs, Point const& rhs)
58 CoordinateType const lh = geometry::get<Dimension>(lhs);
59 CoordinateType const rh = geometry::get<Dimension>(rhs);
61 // If both lhs and rhs equal m_value_to_be_first,
62 // we should handle conform if lh < rh = FALSE
63 // The first condition meets that, keep it first
64 if (geometry::math::equals(rh, m_value_to_be_first))
66 // Handle conform lh < -INF, which is always false
70 if (geometry::math::equals(lh, m_value_to_be_first))
72 // Handle conform -INF < rh, which is always true
80 template <int Dimension, typename Collection, typename Value, typename Predicate>
81 inline bool max_value(Collection const& collection, Value& the_max, Predicate const& predicate)
84 for (typename Collection::const_iterator it = collection.begin(); it != collection.end(); ++it)
88 Value the_value = geometry::get<Dimension>(*std::max_element(it->begin(), it->end(), predicate));
89 if (first || the_value > the_max)
100 template <int Dimension, typename Value>
104 inline select_below(Value const& v)
108 template <typename Intruder>
109 inline bool operator()(Intruder const& intruder) const
111 if (intruder.empty())
115 Value max = geometry::get<Dimension>(*std::max_element(intruder.begin(), intruder.end(), detail::extreme_points::compare<Dimension>()));
116 return geometry::math::equals(max, m_value) || max < m_value;
120 template <int Dimension, typename Value>
124 inline adapt_base(Value const& v)
128 template <typename Intruder>
129 inline void operator()(Intruder& intruder) const
131 if (intruder.size() >= 3)
133 detail::extreme_points::move_along_vector<Dimension>(intruder, m_value);
138 template <int Dimension, typename Value>
139 struct min_of_intruder
141 template <typename Intruder>
142 inline bool operator()(Intruder const& lhs, Intruder const& rhs) const
144 Value lhs_min = geometry::get<Dimension>(*std::min_element(lhs.begin(), lhs.end(), detail::extreme_points::compare<Dimension>()));
145 Value rhs_min = geometry::get<Dimension>(*std::min_element(rhs.begin(), rhs.end(), detail::extreme_points::compare<Dimension>()));
146 return lhs_min < rhs_min;
151 template <typename Point, typename P>
152 inline void calculate_average(Point& point, std::vector<P> const& points)
154 typedef typename geometry::coordinate_type<Point>::type coordinate_type;
155 typedef typename std::vector<P>::const_iterator iterator_type;
156 typedef typename std::vector<P>::size_type size_type;
158 coordinate_type x = 0;
159 coordinate_type y = 0;
161 iterator_type end = points.end();
162 for ( iterator_type it = points.begin() ; it != end ; ++it)
164 x += geometry::get<0>(*it);
165 y += geometry::get<1>(*it);
168 size_type const count = points.size();
169 geometry::set<0>(point, x / count);
170 geometry::set<1>(point, y / count);
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)
177 // This function handles self-tangencies.
178 // Self-tangencies use, as usual, the major part of code...
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:
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)
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))
212 // Throw away all intrusions <= this value, and of the kept one set this as base.
213 select_below<Dimension, CoordinateType> predicate(penultimate_value);
216 std::remove_if(boost::begin(intruders), boost::end(intruders), predicate),
217 boost::end(intruders)
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);
223 // Also adapt base of extremes
224 detail::extreme_points::move_along_vector<Dimension>(extremes, penultimate_value);
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>());
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());
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());
240 // Now replace extremes with this smaller subset, a triangle, such that centroid calculation will result in a point inside
244 template <int Dimension, typename Geometry, typename Point>
245 inline bool calculate_point_on_surface(Geometry const& geometry, Point& point)
247 typedef typename geometry::point_type<Geometry>::type point_type;
248 typedef typename geometry::coordinate_type<Geometry>::type coordinate_type;
249 std::vector<point_type> extremes;
251 typedef std::vector<std::vector<point_type> > intruders_type;
252 intruders_type intruders;
253 geometry::extreme_points<Dimension>(geometry, extremes, intruders);
255 if (extremes.size() < 3)
260 // If there are intruders, find the max.
261 if (! intruders.empty())
263 coordinate_type max_intruder;
264 detail::extreme_points::compare<Dimension> compare;
265 if (max_value<Dimension>(intruders, max_intruder, compare))
267 coordinate_type max_extreme = geometry::get<Dimension>(*std::max_element(extremes.begin(), extremes.end(), detail::extreme_points::compare<Dimension>()));
268 if (max_extreme > max_intruder)
270 detail::extreme_points::move_along_vector<Dimension>(extremes, max_intruder);
274 replace_extremes_for_self_tangencies<Dimension>(extremes, intruders, max_intruder);
279 // Now calculate the average/centroid of the (possibly adapted) extremes
280 calculate_average(point, extremes);
285 }} // namespace detail::point_on_surface
286 #endif // DOXYGEN_NO_DETAIL
290 \brief Assigns a Point guaranteed to lie on the surface of the Geometry
291 \tparam Geometry geometry type. This also defines the type of the output point
292 \param geometry Geometry to take point from
293 \param point Point to assign
295 template <typename Geometry, typename Point>
296 inline void point_on_surface(Geometry const& geometry, Point & point)
298 concepts::check<Point>();
299 concepts::check<Geometry const>();
301 // First try in Y-direction (which should always succeed for valid polygons)
302 if (! detail::point_on_surface::calculate_point_on_surface<1>(geometry, point))
304 // For invalid polygons, we might try X-direction
305 detail::point_on_surface::calculate_point_on_surface<0>(geometry, point);
310 \brief Returns point guaranteed to lie on the surface of the Geometry
311 \tparam Geometry geometry type. This also defines the type of the output point
312 \param geometry Geometry to take point from
313 \return The Point guaranteed to lie on the surface of the Geometry
315 template<typename Geometry>
316 inline typename geometry::point_type<Geometry>::type
317 return_point_on_surface(Geometry const& geometry)
319 typename geometry::point_type<Geometry>::type result;
320 geometry::point_on_surface(geometry, result);
324 }} // namespace boost::geometry
327 #endif // BOOST_GEOMETRY_ALGORITHMS_POINT_ON_SURFACE_HPP