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