]>
Commit | Line | Data |
---|---|---|
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-2014 Adam Wulkiewicz, Lodz, Poland. | |
7 | ||
8 | // Use, modification and distribution is subject to the Boost Software License, | |
9 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
10 | // http://www.boost.org/LICENSE_1_0.txt) | |
11 | ||
12 | #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP | |
13 | #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP | |
14 | ||
15 | ||
16 | #include <cstddef> | |
17 | ||
18 | #include <boost/range.hpp> | |
19 | ||
20 | #include <boost/geometry/algorithms/detail/interior_iterator.hpp> | |
21 | ||
22 | #include <boost/geometry/core/cs.hpp> | |
23 | #include <boost/geometry/core/point_type.hpp> | |
24 | #include <boost/geometry/core/ring_type.hpp> | |
25 | #include <boost/geometry/core/tags.hpp> | |
26 | ||
27 | #include <boost/geometry/geometries/concepts/check.hpp> | |
28 | #include <boost/geometry/iterators/ever_circling_iterator.hpp> | |
29 | ||
30 | #include <boost/geometry/algorithms/detail/assign_box_corners.hpp> | |
31 | ||
32 | #include <boost/geometry/strategies/side.hpp> | |
33 | ||
34 | #include <boost/geometry/util/math.hpp> | |
35 | ||
36 | ||
37 | namespace boost { namespace geometry | |
38 | { | |
39 | ||
40 | ||
41 | #ifndef DOXYGEN_NO_DETAIL | |
42 | namespace detail { namespace extreme_points | |
43 | { | |
44 | ||
45 | template <std::size_t Dimension> | |
46 | struct compare | |
47 | { | |
48 | template <typename Point> | |
49 | inline bool operator()(Point const& lhs, Point const& rhs) | |
50 | { | |
51 | return geometry::get<Dimension>(lhs) < geometry::get<Dimension>(rhs); | |
52 | } | |
53 | }; | |
54 | ||
55 | ||
56 | template <std::size_t Dimension, typename PointType, typename CoordinateType> | |
57 | inline void move_along_vector(PointType& point, PointType const& extreme, CoordinateType const& base_value) | |
58 | { | |
59 | // Moves a point along the vector (point, extreme) in the direction of the extreme point | |
60 | // This adapts the possibly uneven legs of the triangle (or trapezium-like shape) | |
61 | // _____extreme _____ | |
62 | // / \ / \ . | |
63 | // /base \ => / \ point . | |
64 | // \ point . | |
65 | // | |
66 | // For so-called intruders, it can be used to adapt both legs to the level of "base" | |
67 | // For the base, it can be used to adapt both legs to the level of the max-value of the intruders | |
68 | // If there are 2 or more extreme values, use the one close to 'point' to have a correct vector | |
69 | ||
70 | CoordinateType const value = geometry::get<Dimension>(point); | |
71 | //if (geometry::math::equals(value, base_value)) | |
72 | if (value >= base_value) | |
73 | { | |
74 | return; | |
75 | } | |
76 | ||
77 | PointType vector = point; | |
78 | subtract_point(vector, extreme); | |
79 | ||
80 | CoordinateType const diff = geometry::get<Dimension>(vector); | |
81 | ||
82 | // diff should never be zero | |
83 | // because of the way our triangle/trapezium is build. | |
84 | // We just return if it would be the case. | |
85 | if (geometry::math::equals(diff, 0)) | |
86 | { | |
87 | return; | |
88 | } | |
89 | ||
90 | CoordinateType const base_diff = base_value - geometry::get<Dimension>(extreme); | |
91 | ||
92 | multiply_value(vector, base_diff); | |
93 | divide_value(vector, diff); | |
94 | ||
95 | // The real move: | |
96 | point = extreme; | |
97 | add_point(point, vector); | |
98 | } | |
99 | ||
100 | ||
101 | template <std::size_t Dimension, typename Range, typename CoordinateType> | |
102 | inline void move_along_vector(Range& range, CoordinateType const& base_value) | |
103 | { | |
104 | if (range.size() >= 3) | |
105 | { | |
106 | move_along_vector<Dimension>(range.front(), *(range.begin() + 1), base_value); | |
107 | move_along_vector<Dimension>(range.back(), *(range.rbegin() + 1), base_value); | |
108 | } | |
109 | } | |
110 | ||
111 | ||
112 | template<typename Ring, std::size_t Dimension> | |
113 | struct extreme_points_on_ring | |
114 | { | |
115 | ||
116 | typedef typename geometry::coordinate_type<Ring>::type coordinate_type; | |
117 | typedef typename boost::range_iterator<Ring const>::type range_iterator; | |
118 | typedef typename geometry::point_type<Ring>::type point_type; | |
119 | ||
120 | typedef typename geometry::strategy::side::services::default_strategy | |
121 | < | |
122 | typename geometry::cs_tag<point_type>::type | |
123 | >::type side_strategy; | |
124 | ||
125 | ||
126 | template <typename CirclingIterator, typename Points> | |
127 | static inline bool extend(CirclingIterator& it, | |
128 | std::size_t n, | |
129 | coordinate_type max_coordinate_value, | |
130 | Points& points, int direction) | |
131 | { | |
132 | std::size_t safe_index = 0; | |
133 | do | |
134 | { | |
135 | it += direction; | |
136 | ||
137 | points.push_back(*it); | |
138 | ||
139 | if (safe_index++ >= n) | |
140 | { | |
141 | // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop) | |
142 | return false; | |
143 | } | |
144 | } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value)); | |
145 | ||
146 | return true; | |
147 | } | |
148 | ||
149 | // Overload without adding to poinst | |
150 | template <typename CirclingIterator> | |
151 | static inline bool extend(CirclingIterator& it, | |
152 | std::size_t n, | |
153 | coordinate_type max_coordinate_value, | |
154 | int direction) | |
155 | { | |
156 | std::size_t safe_index = 0; | |
157 | do | |
158 | { | |
159 | it += direction; | |
160 | ||
161 | if (safe_index++ >= n) | |
162 | { | |
163 | // E.g.: ring is completely horizontal or vertical (= invalid, but we don't want to have an infinite loop) | |
164 | return false; | |
165 | } | |
166 | } while (geometry::math::equals(geometry::get<Dimension>(*it), max_coordinate_value)); | |
167 | ||
168 | return true; | |
169 | } | |
170 | ||
171 | template <typename CirclingIterator> | |
172 | static inline bool extent_both_sides(Ring const& ring, | |
173 | point_type extreme, | |
174 | CirclingIterator& left, | |
175 | CirclingIterator& right) | |
176 | { | |
177 | std::size_t const n = boost::size(ring); | |
178 | coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme); | |
179 | ||
180 | if (! extend(left, n, max_coordinate_value, -1)) | |
181 | { | |
182 | return false; | |
183 | } | |
184 | if (! extend(right, n, max_coordinate_value, +1)) | |
185 | { | |
186 | return false; | |
187 | } | |
188 | ||
189 | return true; | |
190 | } | |
191 | ||
192 | template <typename Collection, typename CirclingIterator> | |
193 | static inline bool collect(Ring const& ring, | |
194 | point_type extreme, | |
195 | Collection& points, | |
196 | CirclingIterator& left, | |
197 | CirclingIterator& right) | |
198 | { | |
199 | std::size_t const n = boost::size(ring); | |
200 | coordinate_type const max_coordinate_value = geometry::get<Dimension>(extreme); | |
201 | ||
202 | // Collects first left, which is reversed (if more than one point) then adds the top itself, then right | |
203 | if (! extend(left, n, max_coordinate_value, points, -1)) | |
204 | { | |
205 | return false; | |
206 | } | |
207 | std::reverse(points.begin(), points.end()); | |
208 | points.push_back(extreme); | |
209 | if (! extend(right, n, max_coordinate_value, points, +1)) | |
210 | { | |
211 | return false; | |
212 | } | |
213 | ||
214 | return true; | |
215 | } | |
216 | ||
217 | template <typename Extremes, typename Intruders, typename CirclingIterator> | |
218 | static inline void get_intruders(Ring const& ring, CirclingIterator left, CirclingIterator right, | |
219 | Extremes const& extremes, | |
220 | Intruders& intruders) | |
221 | { | |
222 | if (boost::size(extremes) < 3) | |
223 | { | |
224 | return; | |
225 | } | |
226 | coordinate_type const min_value = geometry::get<Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<Dimension>())); | |
227 | ||
228 | // Also select left/right (if Dimension=1) | |
229 | coordinate_type const other_min = geometry::get<1 - Dimension>(*std::min_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>())); | |
230 | coordinate_type const other_max = geometry::get<1 - Dimension>(*std::max_element(boost::begin(extremes), boost::end(extremes), compare<1 - Dimension>())); | |
231 | ||
232 | std::size_t defensive_check_index = 0; // in case we skip over left/right check, collect modifies right too | |
233 | std::size_t const n = boost::size(ring); | |
234 | while (left != right && defensive_check_index < n) | |
235 | { | |
236 | coordinate_type const coordinate = geometry::get<Dimension>(*right); | |
237 | coordinate_type const other_coordinate = geometry::get<1 - Dimension>(*right); | |
238 | if (coordinate > min_value && other_coordinate > other_min && other_coordinate < other_max) | |
239 | { | |
240 | int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1; | |
241 | int const first_side = side_strategy::apply(*right, extremes.front(), *(extremes.begin() + 1)) * factor; | |
242 | int const last_side = side_strategy::apply(*right, *(extremes.rbegin() + 1), extremes.back()) * factor; | |
243 | ||
244 | // If not lying left from any of the extemes side | |
245 | if (first_side != 1 && last_side != 1) | |
246 | { | |
247 | //std::cout << "first " << first_side << " last " << last_side << std::endl; | |
248 | ||
249 | // we start at this intrusion until it is handled, and don't affect our initial left iterator | |
250 | CirclingIterator left_intrusion_it = right; | |
251 | typename boost::range_value<Intruders>::type intruder; | |
252 | collect(ring, *right, intruder, left_intrusion_it, right); | |
253 | ||
254 | // Also moves these to base-level, makes sorting possible which can be done in case of self-tangencies | |
255 | // (we might postpone this action, it is often not necessary. However it is not time-consuming) | |
256 | move_along_vector<Dimension>(intruder, min_value); | |
257 | intruders.push_back(intruder); | |
258 | --right; | |
259 | } | |
260 | } | |
261 | ++right; | |
262 | defensive_check_index++; | |
263 | } | |
264 | } | |
265 | ||
266 | template <typename Extremes, typename Intruders> | |
267 | static inline void get_intruders(Ring const& ring, | |
268 | Extremes const& extremes, | |
269 | Intruders& intruders) | |
270 | { | |
271 | std::size_t const n = boost::size(ring); | |
272 | if (n >= 3) | |
273 | { | |
274 | geometry::ever_circling_range_iterator<Ring const> left(ring); | |
275 | geometry::ever_circling_range_iterator<Ring const> right(ring); | |
276 | ++right; | |
277 | ||
278 | get_intruders(ring, left, right, extremes, intruders); | |
279 | } | |
280 | } | |
281 | ||
282 | template <typename Iterator> | |
283 | static inline bool right_turn(Ring const& ring, Iterator it) | |
284 | { | |
285 | typename std::iterator_traits<Iterator>::difference_type const index | |
286 | = std::distance(boost::begin(ring), it); | |
287 | geometry::ever_circling_range_iterator<Ring const> left(ring); | |
288 | geometry::ever_circling_range_iterator<Ring const> right(ring); | |
289 | left += index; | |
290 | right += index; | |
291 | ||
292 | if (! extent_both_sides(ring, *it, left, right)) | |
293 | { | |
294 | return false; | |
295 | } | |
296 | ||
297 | int const factor = geometry::point_order<Ring>::value == geometry::clockwise ? 1 : -1; | |
298 | int const first_side = side_strategy::apply(*(right - 1), *right, *left) * factor; | |
299 | int const last_side = side_strategy::apply(*left, *(left + 1), *right) * factor; | |
300 | ||
301 | //std::cout << "Candidate at " << geometry::wkt(*it) << " first=" << first_side << " last=" << last_side << std::endl; | |
302 | ||
303 | // Turn should not be left (actually, it should be right because extent removes horizontal/collinear cases) | |
304 | return first_side != 1 && last_side != 1; | |
305 | } | |
306 | ||
307 | ||
308 | // Gets the extreme segments (top point plus neighbouring points), plus intruders, if any, on the same ring | |
309 | template <typename Extremes, typename Intruders> | |
310 | static inline bool apply(Ring const& ring, Extremes& extremes, Intruders& intruders) | |
311 | { | |
312 | std::size_t const n = boost::size(ring); | |
313 | if (n < 3) | |
314 | { | |
315 | return false; | |
316 | } | |
317 | ||
318 | // Get all maxima, usually one. In case of self-tangencies, or self-crossings, | |
319 | // the max might be is not valid. A valid max should make a right turn | |
320 | range_iterator max_it = boost::begin(ring); | |
321 | compare<Dimension> smaller; | |
322 | for (range_iterator it = max_it + 1; it != boost::end(ring); ++it) | |
323 | { | |
324 | if (smaller(*max_it, *it) && right_turn(ring, it)) | |
325 | { | |
326 | max_it = it; | |
327 | } | |
328 | } | |
329 | ||
330 | if (max_it == boost::end(ring)) | |
331 | { | |
332 | return false; | |
333 | } | |
334 | ||
335 | typename std::iterator_traits<range_iterator>::difference_type const | |
336 | index = std::distance(boost::begin(ring), max_it); | |
337 | //std::cout << "Extreme point lies at " << index << " having " << geometry::wkt(*max_it) << std::endl; | |
338 | ||
339 | geometry::ever_circling_range_iterator<Ring const> left(ring); | |
340 | geometry::ever_circling_range_iterator<Ring const> right(ring); | |
341 | left += index; | |
342 | right += index; | |
343 | ||
344 | // Collect all points (often 3) in a temporary vector | |
345 | std::vector<point_type> points; | |
346 | points.reserve(3); | |
347 | if (! collect(ring, *max_it, points, left, right)) | |
348 | { | |
349 | return false; | |
350 | } | |
351 | ||
352 | //std::cout << "Built vector of " << points.size() << std::endl; | |
353 | ||
354 | coordinate_type const front_value = geometry::get<Dimension>(points.front()); | |
355 | coordinate_type const back_value = geometry::get<Dimension>(points.back()); | |
356 | coordinate_type const base_value = (std::max)(front_value, back_value); | |
357 | if (front_value < back_value) | |
358 | { | |
359 | move_along_vector<Dimension>(points.front(), *(points.begin() + 1), base_value); | |
360 | } | |
361 | else | |
362 | { | |
363 | move_along_vector<Dimension>(points.back(), *(points.rbegin() + 1), base_value); | |
364 | } | |
365 | ||
366 | std::copy(points.begin(), points.end(), std::back_inserter(extremes)); | |
367 | ||
368 | get_intruders(ring, left, right, extremes, intruders); | |
369 | ||
370 | return true; | |
371 | } | |
372 | }; | |
373 | ||
374 | ||
375 | ||
376 | ||
377 | ||
378 | }} // namespace detail::extreme_points | |
379 | #endif // DOXYGEN_NO_DETAIL | |
380 | ||
381 | ||
382 | #ifndef DOXYGEN_NO_DISPATCH | |
383 | namespace dispatch | |
384 | { | |
385 | ||
386 | ||
387 | template | |
388 | < | |
389 | typename Geometry, | |
390 | std::size_t Dimension, | |
391 | typename GeometryTag = typename tag<Geometry>::type | |
392 | > | |
393 | struct extreme_points | |
394 | {}; | |
395 | ||
396 | ||
397 | template<typename Ring, std::size_t Dimension> | |
398 | struct extreme_points<Ring, Dimension, ring_tag> | |
399 | : detail::extreme_points::extreme_points_on_ring<Ring, Dimension> | |
400 | {}; | |
401 | ||
402 | ||
403 | template<typename Polygon, std::size_t Dimension> | |
404 | struct extreme_points<Polygon, Dimension, polygon_tag> | |
405 | { | |
406 | template <typename Extremes, typename Intruders> | |
407 | static inline bool apply(Polygon const& polygon, Extremes& extremes, Intruders& intruders) | |
408 | { | |
409 | typedef typename geometry::ring_type<Polygon>::type ring_type; | |
410 | typedef detail::extreme_points::extreme_points_on_ring | |
411 | < | |
412 | ring_type, Dimension | |
413 | > ring_implementation; | |
414 | ||
415 | if (! ring_implementation::apply(geometry::exterior_ring(polygon), extremes, intruders)) | |
416 | { | |
417 | return false; | |
418 | } | |
419 | ||
420 | // For a polygon, its interior rings can contain intruders | |
421 | typename interior_return_type<Polygon const>::type | |
422 | rings = interior_rings(polygon); | |
423 | for (typename detail::interior_iterator<Polygon const>::type | |
424 | it = boost::begin(rings); it != boost::end(rings); ++it) | |
425 | { | |
426 | ring_implementation::get_intruders(*it, extremes, intruders); | |
427 | } | |
428 | ||
429 | return true; | |
430 | } | |
431 | }; | |
432 | ||
433 | template<typename Box> | |
434 | struct extreme_points<Box, 1, box_tag> | |
435 | { | |
436 | template <typename Extremes, typename Intruders> | |
437 | static inline bool apply(Box const& box, Extremes& extremes, Intruders& ) | |
438 | { | |
439 | extremes.resize(4); | |
440 | geometry::detail::assign_box_corners_oriented<false>(box, extremes); | |
441 | // ll,ul,ur,lr, contains too exactly the right info | |
442 | return true; | |
443 | } | |
444 | }; | |
445 | ||
446 | template<typename Box> | |
447 | struct extreme_points<Box, 0, box_tag> | |
448 | { | |
449 | template <typename Extremes, typename Intruders> | |
450 | static inline bool apply(Box const& box, Extremes& extremes, Intruders& ) | |
451 | { | |
452 | extremes.resize(4); | |
453 | geometry::detail::assign_box_corners_oriented<false>(box, extremes); | |
454 | // ll,ul,ur,lr, rotate one to start with UL and end with LL | |
455 | std::rotate(extremes.begin(), extremes.begin() + 1, extremes.end()); | |
456 | return true; | |
457 | } | |
458 | }; | |
459 | ||
460 | template<typename MultiPolygon, std::size_t Dimension> | |
461 | struct extreme_points<MultiPolygon, Dimension, multi_polygon_tag> | |
462 | { | |
463 | template <typename Extremes, typename Intruders> | |
464 | static inline bool apply(MultiPolygon const& multi, Extremes& extremes, Intruders& intruders) | |
465 | { | |
466 | // Get one for the very first polygon, that is (for the moment) enough. | |
467 | // It is not guaranteed the "extreme" then, but for the current purpose | |
468 | // (point_on_surface) it can just be this point. | |
469 | if (boost::size(multi) >= 1) | |
470 | { | |
471 | return extreme_points | |
472 | < | |
473 | typename boost::range_value<MultiPolygon const>::type, | |
474 | Dimension, | |
475 | polygon_tag | |
476 | >::apply(*boost::begin(multi), extremes, intruders); | |
477 | } | |
478 | ||
479 | return false; | |
480 | } | |
481 | }; | |
482 | ||
483 | } // namespace dispatch | |
484 | #endif // DOXYGEN_NO_DISPATCH | |
485 | ||
486 | ||
487 | /*! | |
488 | \brief Returns extreme points (for Edge=1 in dimension 1, so the top, | |
489 | for Edge=0 in dimension 0, the right side) | |
490 | \note We could specify a strategy (less/greater) to get bottom/left side too. However, until now we don't need that. | |
491 | */ | |
492 | template <std::size_t Edge, typename Geometry, typename Extremes, typename Intruders> | |
493 | inline bool extreme_points(Geometry const& geometry, Extremes& extremes, Intruders& intruders) | |
494 | { | |
495 | concepts::check<Geometry const>(); | |
496 | ||
497 | // Extremes is not required to follow a geometry concept (but it should support an output iterator), | |
498 | // but its elements should fulfil the point-concept | |
499 | concepts::check<typename boost::range_value<Extremes>::type>(); | |
500 | ||
501 | // Intruders should contain collections which value type is point-concept | |
502 | // Extremes might be anything (supporting an output iterator), but its elements should fulfil the point-concept | |
503 | concepts::check | |
504 | < | |
505 | typename boost::range_value | |
506 | < | |
507 | typename boost::range_value<Intruders>::type | |
508 | >::type | |
509 | const | |
510 | >(); | |
511 | ||
512 | return dispatch::extreme_points<Geometry, Edge>::apply(geometry, extremes, intruders); | |
513 | } | |
514 | ||
515 | ||
516 | ||
517 | }} // namespace boost::geometry | |
518 | ||
519 | ||
520 | #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_EXTREME_POINTS_HPP |