• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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