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