1 // Boost.Geometry (aka GGL, Generic Geometry Library)
3 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
4 // Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
6 // This file was modified by Oracle on 2017, 2019.
7 // Modifications copyright (c) 2017, 2019 Oracle and/or its affiliates.
9 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
11 // Use, modification and distribution is subject to the Boost Software License,
12 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
13 // http://www.boost.org/LICENSE_1_0.txt)
15 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP
16 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP
22 #include <boost/geometry/algorithms/detail/overlay/approximately_equals.hpp>
23 #include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp>
24 #include <boost/geometry/algorithms/detail/overlay/get_ring.hpp>
25 #include <boost/geometry/algorithms/detail/direction_code.hpp>
26 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
28 #include <boost/geometry/util/condition.hpp>
29 #include <boost/geometry/util/math.hpp>
30 #include <boost/geometry/util/select_coordinate_type.hpp>
31 #include <boost/geometry/util/select_most_precise.hpp>
33 namespace boost { namespace geometry
36 #ifndef DOXYGEN_NO_DETAIL
37 namespace detail { namespace overlay { namespace sort_by_side
40 enum direction_type { dir_unknown = -1, dir_from = 0, dir_to = 1 };
42 typedef signed_size_type rank_type;
45 // Point-wrapper, adding some properties
46 template <typename Point>
53 , direction(dir_unknown)
56 , operation(operation_none)
59 ranked_point(Point const& p, signed_size_type ti, int oi,
60 direction_type d, operation_type op, segment_identifier const& si)
73 using point_type = Point;
77 signed_size_type zone; // index of closed zone, in uu turn there would be 2 zones
78 signed_size_type turn_index;
79 int operation_index; // 0,1
80 direction_type direction;
81 std::size_t count_left;
82 std::size_t count_right;
83 operation_type operation;
84 segment_identifier seg_id;
87 struct less_by_turn_index
90 inline bool operator()(const T& first, const T& second) const
92 return first.turn_index == second.turn_index
93 ? first.index < second.index
94 : first.turn_index < second.turn_index
101 template <typename T>
102 inline bool operator()(const T& first, const T& second) const
104 // Length might be considered too
105 // First order by from/to
106 if (first.direction != second.direction)
108 return first.direction < second.direction;
110 // Then by turn index
111 if (first.turn_index != second.turn_index)
113 return first.turn_index < second.turn_index;
115 // This can also be the same (for example in buffer), but seg_id is
117 return first.seg_id < second.seg_id;
123 template <typename T>
124 inline bool operator()(const T&, const T& ) const
130 template <typename Point, typename SideStrategy, typename LessOnSame, typename Compare>
133 less_by_side(const Point& p1, const Point& p2, SideStrategy const& strategy)
136 , m_strategy(strategy)
139 template <typename T>
140 inline bool operator()(const T& first, const T& second) const
142 typedef typename SideStrategy::cs_tag cs_tag;
147 int const side_first = m_strategy.apply(m_origin, m_turn_point, first.point);
148 int const side_second = m_strategy.apply(m_origin, m_turn_point, second.point);
150 if (side_first == 0 && side_second == 0)
152 // Both collinear. They might point into different directions: <------*------>
153 // If so, order the one going backwards as the very first.
155 int const first_code = direction_code<cs_tag>(m_origin, m_turn_point, first.point);
156 int const second_code = direction_code<cs_tag>(m_origin, m_turn_point, second.point);
158 // Order by code, backwards first, then forward.
159 return first_code != second_code
160 ? first_code < second_code
161 : on_same(first, second)
164 else if (side_first == 0
165 && direction_code<cs_tag>(m_origin, m_turn_point, first.point) == -1)
167 // First collinear and going backwards.
168 // Order as the very first, so return always true
171 else if (side_second == 0
172 && direction_code<cs_tag>(m_origin, m_turn_point, second.point) == -1)
174 // Second is collinear and going backwards
175 // Order as very last, so return always false
179 // They are not both collinear
181 if (side_first != side_second)
183 return compare(side_first, side_second);
186 // They are both left, both right, and/or both collinear (with each other and/or with p1,p2)
188 int const side_second_wrt_first = m_strategy.apply(m_turn_point, first.point, second.point);
190 if (side_second_wrt_first == 0)
192 return on_same(first, second);
195 int const side_first_wrt_second = m_strategy.apply(m_turn_point, second.point, first.point);
196 if (side_second_wrt_first != -side_first_wrt_second)
198 // (FP) accuracy error in side calculation, the sides are not opposite.
199 // In that case they can be handled as collinear.
200 // If not, then the sort-order might not be stable.
201 return on_same(first, second);
204 // Both are on same side, and not collinear
205 // Union: return true if second is right w.r.t. first, so -1,
206 // so other is 1. union has greater as compare functor
207 // Intersection: v.v.
208 return compare(side_first_wrt_second, side_second_wrt_first);
212 Point const& m_origin;
213 Point const& m_turn_point;
214 SideStrategy const& m_strategy;
217 // Sorts vectors in counter clockwise order (by default)
222 overlay_type OverlayType,
224 typename SideStrategy,
229 typedef ranked_point<Point> rp;
234 inline bool operator()(rp const& ranked_point) const
236 // New candidate if there are no polygons on left side,
237 // but there are on right side
238 return ranked_point.count_left == 0
239 && ranked_point.count_right > 0;
243 struct include_intersection
245 inline bool operator()(rp const& ranked_point) const
247 // New candidate if there are two polygons on right side,
248 // and less on the left side
249 return ranked_point.count_left < 2
250 && ranked_point.count_right >= 2;
255 side_sorter(SideStrategy const& strategy)
257 , m_origin_segment_distance(0)
258 , m_strategy(strategy)
261 template <typename Operation>
262 void add_segment_from(signed_size_type turn_index, int op_index,
263 Point const& point_from,
267 m_ranked_points.push_back(rp(point_from, turn_index, op_index,
268 dir_from, op.operation, op.seg_id));
271 m_origin = point_from;
276 template <typename Operation>
277 void add_segment_to(signed_size_type turn_index, int op_index,
278 Point const& point_to,
281 m_ranked_points.push_back(rp(point_to, turn_index, op_index,
282 dir_to, op.operation, op.seg_id));
285 template <typename Operation>
286 void add_segment(signed_size_type turn_index, int op_index,
287 Point const& point_from, Point const& point_to,
288 Operation const& op, bool is_origin)
290 add_segment_from(turn_index, op_index, point_from, op, is_origin);
291 add_segment_to(turn_index, op_index, point_to, op);
294 template <typename Operation, typename Geometry1, typename Geometry2>
295 static Point walk_over_ring(Operation const& op, int offset,
296 Geometry1 const& geometry1,
297 Geometry2 const& geometry2)
300 geometry::copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, op.seg_id, offset, point);
304 template <typename Turn, typename Operation, typename Geometry1, typename Geometry2>
305 Point add(Turn const& turn, Operation const& op, signed_size_type turn_index, int op_index,
306 Geometry1 const& geometry1,
307 Geometry2 const& geometry2,
310 Point point_from, point2, point3;
311 geometry::copy_segment_points<Reverse1, Reverse2>(geometry1, geometry2,
312 op.seg_id, point_from, point2, point3);
313 Point point_to = op.fraction.is_one() ? point3 : point2;
315 // If the point is in the neighbourhood (the limit is arbitrary),
316 // then take a point (or more) further back.
317 // The limit of offset avoids theoretical infinite loops.
318 // In practice it currently walks max 1 point back in all cases.
319 // Use the coordinate type, but if it is too small (e.g. std::int16), use a double
320 using ct_type = typename geometry::select_most_precise
322 typename geometry::coordinate_type<Point>::type,
326 ct_type const tolerance = 1000000000;
329 while (approximately_equals(point_from, turn.point, tolerance)
332 point_from = walk_over_ring(op, --offset, geometry1, geometry2);
335 // Similarly for the point_to, walk forward
337 while (approximately_equals(point_to, turn.point, tolerance)
340 point_to = walk_over_ring(op, ++offset, geometry1, geometry2);
343 add_segment(turn_index, op_index, point_from, point_to, op, is_origin);
348 template <typename Turn, typename Operation, typename Geometry1, typename Geometry2>
349 void add(Turn const& turn,
350 Operation const& op, signed_size_type turn_index, int op_index,
351 segment_identifier const& departure_seg_id,
352 Geometry1 const& geometry1,
353 Geometry2 const& geometry2,
356 Point potential_origin = add(turn, op, turn_index, op_index, geometry1, geometry2, false);
361 = op.seg_id.source_index == departure_seg_id.source_index
362 && op.seg_id.ring_index == departure_seg_id.ring_index
363 && op.seg_id.multi_index == departure_seg_id.multi_index;
367 signed_size_type const sd
368 = departure_seg_id.source_index == 0
369 ? segment_distance(geometry1, departure_seg_id, op.seg_id)
370 : segment_distance(geometry2, departure_seg_id, op.seg_id);
372 if (m_origin_count == 0 || sd < m_origin_segment_distance)
374 m_origin = potential_origin;
375 m_origin_segment_distance = sd;
382 void apply(Point const& turn_point)
384 // We need three compare functors:
385 // 1) to order clockwise (union) or counter clockwise (intersection)
386 // 2) to order by side, resulting in unique ranks for all points
387 // 3) to order by side, resulting in non-unique ranks
388 // to give colinear points
390 // Sort by side and assign rank
391 less_by_side<Point, SideStrategy, less_by_index, Compare> less_unique(m_origin, turn_point, m_strategy);
392 less_by_side<Point, SideStrategy, less_false, Compare> less_non_unique(m_origin, turn_point, m_strategy);
394 std::sort(m_ranked_points.begin(), m_ranked_points.end(), less_unique);
396 std::size_t colinear_rank = 0;
397 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
400 && less_non_unique(m_ranked_points[i - 1], m_ranked_points[i]))
402 // It is not collinear
406 m_ranked_points[i].rank = colinear_rank;
410 void find_open_by_piece_index()
412 // For buffers, use piece index
413 std::set<signed_size_type> handled;
415 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
417 const rp& ranked = m_ranked_points[i];
418 if (ranked.direction != dir_from)
423 signed_size_type const& index = ranked.seg_id.piece_index;
424 if (handled.count(index) > 0)
428 find_polygons_for_source<&segment_identifier::piece_index>(index, i);
429 handled.insert(index);
433 void find_open_by_source_index()
435 // Check for source index 0 and 1
436 bool handled[2] = {false, false};
437 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
439 const rp& ranked = m_ranked_points[i];
440 if (ranked.direction != dir_from)
445 signed_size_type const& index = ranked.seg_id.source_index;
446 if (index < 0 || index > 1 || handled[index])
450 find_polygons_for_source<&segment_identifier::source_index>(index, i);
451 handled[index] = true;
457 if (BOOST_GEOMETRY_CONDITION(OverlayType == overlay_buffer))
459 find_open_by_piece_index();
463 find_open_by_source_index();
469 if (m_ranked_points.empty())
474 std::size_t const last = 1 + m_ranked_points.back().rank;
476 // Move iterator after rank==0
477 bool has_first = false;
478 typename container_type::iterator it = m_ranked_points.begin() + 1;
479 for (; it != m_ranked_points.end() && it->rank == 0; ++it)
486 // Reverse first part (having rank == 0), if any,
487 // but skip the very first row
488 std::reverse(m_ranked_points.begin() + 1, it);
489 for (typename container_type::iterator fit = m_ranked_points.begin();
492 BOOST_ASSERT(fit->rank == 0);
496 // Reverse the rest (main rank > 0)
497 std::reverse(it, m_ranked_points.end());
498 for (; it != m_ranked_points.end(); ++it)
500 BOOST_ASSERT(it->rank > 0);
501 it->rank = last - it->rank;
505 bool has_origin() const
507 return m_origin_count > 0;
512 typedef std::vector<rp> container_type;
513 container_type m_ranked_points;
515 std::size_t m_origin_count;
516 signed_size_type m_origin_segment_distance;
517 SideStrategy m_strategy;
521 //! Check how many open spaces there are
522 template <typename Include>
523 inline std::size_t open_count(Include const& include_functor) const
525 std::size_t result = 0;
526 rank_type last_rank = 0;
527 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
529 rp const& ranked_point = m_ranked_points[i];
531 if (ranked_point.rank > last_rank
532 && ranked_point.direction == sort_by_side::dir_to
533 && include_functor(ranked_point))
536 last_rank = ranked_point.rank;
542 std::size_t move(std::size_t index) const
544 std::size_t const result = index + 1;
545 return result >= m_ranked_points.size() ? 0 : result;
548 //! member is pointer to member (source_index or multi_index)
549 template <signed_size_type segment_identifier::*Member>
550 std::size_t move(signed_size_type member_index, std::size_t index) const
552 std::size_t result = move(index);
553 while (m_ranked_points[result].seg_id.*Member != member_index)
555 result = move(result);
560 void assign_ranks(rank_type min_rank, rank_type max_rank, int side_index)
562 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
564 rp& ranked = m_ranked_points[i];
565 // Suppose there are 8 ranks, if min=4,max=6: assign 4,5,6
566 // if min=5,max=2: assign from 5,6,7,1,2
568 = max_rank >= min_rank
569 ? ranked.rank >= min_rank && ranked.rank <= max_rank
570 : ranked.rank >= min_rank || ranked.rank <= max_rank
579 else if (side_index == 2)
581 ranked.count_right++;
587 template <signed_size_type segment_identifier::*Member>
588 void find_polygons_for_source(signed_size_type the_index,
589 std::size_t start_index)
591 bool in_polygon = true; // Because start_index is "from", arrives at the turn
592 rp const& start_rp = m_ranked_points[start_index];
593 rank_type last_from_rank = start_rp.rank;
594 rank_type previous_rank = start_rp.rank;
596 for (std::size_t index = move<Member>(the_index, start_index);
598 index = move<Member>(the_index, index))
600 rp& ranked = m_ranked_points[index];
602 if (ranked.rank != previous_rank && ! in_polygon)
604 assign_ranks(last_from_rank, previous_rank - 1, 1);
605 assign_ranks(last_from_rank + 1, previous_rank, 2);
608 if (index == start_index)
613 if (ranked.direction == dir_from)
615 last_from_rank = ranked.rank;
618 else if (ranked.direction == dir_to)
623 previous_rank = ranked.rank;
627 //! Find closed zones and assign it
628 template <typename Include>
629 std::size_t assign_zones(Include const& include_functor)
631 // Find a starting point (the first rank after an outgoing rank
632 // with no polygons on the left side)
633 rank_type start_rank = m_ranked_points.size() + 1;
634 std::size_t start_index = 0;
635 rank_type max_rank = 0;
636 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
638 rp const& ranked_point = m_ranked_points[i];
639 if (ranked_point.rank > max_rank)
641 max_rank = ranked_point.rank;
643 if (ranked_point.direction == sort_by_side::dir_to
644 && include_functor(ranked_point))
646 start_rank = ranked_point.rank + 1;
648 if (ranked_point.rank == start_rank && start_index == 0)
655 rank_type const undefined_rank = max_rank + 1;
656 std::size_t zone_id = 0;
657 rank_type last_rank = 0;
658 rank_type rank_at_next_zone = undefined_rank;
659 std::size_t index = start_index;
660 for (std::size_t i = 0; i < m_ranked_points.size(); i++)
662 rp& ranked_point = m_ranked_points[index];
664 // Implement cyclic behavior
666 if (index == m_ranked_points.size())
671 if (ranked_point.rank != last_rank)
673 if (ranked_point.rank == rank_at_next_zone)
676 rank_at_next_zone = undefined_rank;
679 if (ranked_point.direction == sort_by_side::dir_to
680 && include_functor(ranked_point))
682 rank_at_next_zone = ranked_point.rank + 1;
683 if (rank_at_next_zone > max_rank)
685 rank_at_next_zone = 0;
689 last_rank = ranked_point.rank;
692 ranked_point.zone = zone_id;
698 inline std::size_t open_count(operation_type for_operation) const
700 return for_operation == operation_union
701 ? open_count(include_union())
702 : open_count(include_intersection())
706 inline std::size_t assign_zones(operation_type for_operation)
708 return for_operation == operation_union
709 ? assign_zones(include_union())
710 : assign_zones(include_intersection())
717 //! Metafunction to define side_order (clockwise, ccw) by operation_type
718 template <operation_type OpType>
719 struct side_compare {};
722 struct side_compare<operation_union>
724 typedef std::greater<int> type;
728 struct side_compare<operation_intersection>
730 typedef std::less<int> type;
734 }}} // namespace detail::overlay::sort_by_side
735 #endif //DOXYGEN_NO_DETAIL
738 }} // namespace boost::geometry
740 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP