1 // Boost.Geometry (aka GGL, Generic Geometry Library)
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
5 // This file was modified by Oracle on 2013, 2014, 2015.
6 // Modifications copyright (c) 2013-2015 Oracle and/or its affiliates.
8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
10 // Use, modification and distribution is subject to the Boost Software License,
11 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
12 // http://www.boost.org/LICENSE_1_0.txt)
14 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP
15 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP
17 #include <boost/core/ignore_unused.hpp>
18 #include <boost/range/size.hpp>
20 #include <boost/geometry/core/assert.hpp>
21 #include <boost/geometry/core/topological_dimension.hpp>
23 #include <boost/geometry/util/condition.hpp>
24 #include <boost/geometry/util/range.hpp>
26 #include <boost/geometry/algorithms/num_interior_rings.hpp>
27 #include <boost/geometry/algorithms/detail/point_on_border.hpp>
28 #include <boost/geometry/algorithms/detail/sub_range.hpp>
29 #include <boost/geometry/algorithms/detail/single_geometry.hpp>
31 #include <boost/geometry/algorithms/detail/relate/point_geometry.hpp>
32 #include <boost/geometry/algorithms/detail/relate/turns.hpp>
33 #include <boost/geometry/algorithms/detail/relate/boundary_checker.hpp>
34 #include <boost/geometry/algorithms/detail/relate/follow_helpers.hpp>
36 #include <boost/geometry/views/detail/normalized_view.hpp>
38 namespace boost { namespace geometry
41 #ifndef DOXYGEN_NO_DETAIL
42 namespace detail { namespace relate {
45 // TODO: In the worst case calling this Pred in a loop for MultiLinestring/MultiPolygon may take O(NM)
46 // Use the rtree in this case!
48 // may be used to set IE and BE for a Linear geometry for which no turns were generated
49 template <typename Geometry2, typename Result, typename BoundaryChecker, bool TransposeResult>
50 class no_turns_la_linestring_pred
53 no_turns_la_linestring_pred(Geometry2 const& geometry2,
55 BoundaryChecker const& boundary_checker)
56 : m_geometry2(geometry2)
58 , m_boundary_checker(boundary_checker)
59 , m_interrupt_flags(0)
61 if ( ! may_update<interior, interior, '1', TransposeResult>(m_result) )
63 m_interrupt_flags |= 1;
66 if ( ! may_update<interior, exterior, '1', TransposeResult>(m_result) )
68 m_interrupt_flags |= 2;
71 if ( ! may_update<boundary, interior, '0', TransposeResult>(m_result) )
73 m_interrupt_flags |= 4;
76 if ( ! may_update<boundary, exterior, '0', TransposeResult>(m_result) )
78 m_interrupt_flags |= 8;
82 template <typename Linestring>
83 bool operator()(Linestring const& linestring)
85 std::size_t const count = boost::size(linestring);
91 // TODO: throw an exception?
95 // if those flags are set nothing will change
96 if ( m_interrupt_flags == 0xF )
101 int const pig = detail::within::point_in_geometry(range::front(linestring), m_geometry2);
102 //BOOST_GEOMETRY_ASSERT_MSG(pig != 0, "There should be no IPs");
106 update<interior, interior, '1', TransposeResult>(m_result);
107 m_interrupt_flags |= 1;
111 update<interior, exterior, '1', TransposeResult>(m_result);
112 m_interrupt_flags |= 2;
115 // check if there is a boundary
116 if ( ( m_interrupt_flags & 0xC ) != 0xC // if wasn't already set
117 && ( m_boundary_checker.template
118 is_endpoint_boundary<boundary_front>(range::front(linestring))
119 || m_boundary_checker.template
120 is_endpoint_boundary<boundary_back>(range::back(linestring)) ) )
124 update<boundary, interior, '0', TransposeResult>(m_result);
125 m_interrupt_flags |= 4;
129 update<boundary, exterior, '0', TransposeResult>(m_result);
130 m_interrupt_flags |= 8;
134 return m_interrupt_flags != 0xF
135 && ! m_result.interrupt;
139 Geometry2 const& m_geometry2;
141 BoundaryChecker const& m_boundary_checker;
142 unsigned m_interrupt_flags;
145 // may be used to set EI and EB for an Areal geometry for which no turns were generated
146 template <typename Result, bool TransposeResult>
147 class no_turns_la_areal_pred
150 no_turns_la_areal_pred(Result & res)
152 , interrupt(! may_update<interior, exterior, '2', TransposeResult>(m_result)
153 && ! may_update<boundary, exterior, '1', TransposeResult>(m_result) )
156 template <typename Areal>
157 bool operator()(Areal const& areal)
165 // handle empty/invalid geometries in a different way than below?
167 typedef typename geometry::point_type<Areal>::type point_type;
169 bool const ok = boost::geometry::point_on_border(dummy, areal);
171 // TODO: for now ignore, later throw an exception?
177 update<interior, exterior, '2', TransposeResult>(m_result);
178 update<boundary, exterior, '1', TransposeResult>(m_result);
185 bool const interrupt;
188 // The implementation of an algorithm calculating relate() for L/A
189 template <typename Geometry1, typename Geometry2, bool TransposeResult = false>
192 // check Linear / Areal
193 BOOST_STATIC_ASSERT(topological_dimension<Geometry1>::value == 1
194 && topological_dimension<Geometry2>::value == 2);
196 static const bool interruption_enabled = true;
198 typedef typename geometry::point_type<Geometry1>::type point1_type;
199 typedef typename geometry::point_type<Geometry2>::type point2_type;
201 template <typename Geometry>
206 typename tag<Geometry>::type
210 template <typename Geom1, typename Geom2>
211 struct multi_turn_info
212 : turns::get_turns<Geom1, Geom2>::turn_info
214 multi_turn_info() : priority(0) {}
215 int priority; // single-geometry sorting priority
218 template <typename Geom1, typename Geom2>
219 struct turn_info_type
222 is_multi<Geometry2>::value,
223 multi_turn_info<Geom1, Geom2>,
224 typename turns::get_turns<Geom1, Geom2>::turn_info
228 template <typename Result>
229 static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result)
231 // TODO: If Areal geometry may have infinite size, change the following line:
233 // The result should be FFFFFFFFF
234 relate::set<exterior, exterior, result_dimension<Geometry2>::value, TransposeResult>(result);// FFFFFFFFd, d in [1,9] or T
236 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
239 // get and analyse turns
240 typedef typename turn_info_type<Geometry1, Geometry2>::type turn_type;
241 std::vector<turn_type> turns;
243 interrupt_policy_linear_areal<Geometry2, Result> interrupt_policy(geometry2, result);
245 turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy);
246 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
249 boundary_checker<Geometry1> boundary_checker1(geometry1);
250 no_turns_la_linestring_pred
254 boundary_checker<Geometry1>,
256 > pred1(geometry2, result, boundary_checker1);
257 for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
258 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
261 no_turns_la_areal_pred<Result, !TransposeResult> pred2(result);
262 for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
263 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
269 // This is set here because in the case if empty Areal geometry were passed
270 // those shouldn't be set
271 relate::set<exterior, interior, '2', TransposeResult>(result);// FFFFFF2Fd
272 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
276 sort_dispatch(turns.begin(), turns.end(), is_multi<Geometry2>());
278 turns_analyser<turn_type> analyser;
279 analyse_each_turn(result, analyser,
280 turns.begin(), turns.end(),
281 geometry1, geometry2,
284 if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
288 // If 'c' (insersection_boundary) was not found we know that any Ls isn't equal to one of the Rings
289 if ( !interrupt_policy.is_boundary_found )
291 relate::set<exterior, boundary, '1', TransposeResult>(result);
293 // Don't calculate it if it's required
294 else if ( may_update<exterior, boundary, '1', TransposeResult>(result) )
296 // TODO: REVISE THIS CODE AND PROBABLY REWRITE SOME PARTS TO BE MORE HUMAN-READABLE
297 // IN GENERAL IT ANALYSES THE RINGS OF AREAL GEOMETRY AND DETECTS THE ONES THAT
298 // MAY OVERLAP THE INTERIOR OF LINEAR GEOMETRY (NO IPs OR NON-FAKE 'u' OPERATION)
299 // NOTE: For one case std::sort may be called again to sort data by operations for data already sorted by ring index
300 // In the worst case scenario the complexity will be O( NlogN + R*(N/R)log(N/R) )
301 // So always should remain O(NlogN) -> for R==1 <-> 1(N/1)log(N/1), for R==N <-> N(N/N)log(N/N)
302 // Some benchmarking should probably be done to check if only one std::sort should be used
304 // sort by multi_index and rind_index
305 std::sort(turns.begin(), turns.end(), less_ring());
307 typedef typename std::vector<turn_type>::iterator turn_iterator;
309 turn_iterator it = turns.begin();
310 segment_identifier * prev_seg_id_ptr = NULL;
312 for ( ; it != turns.end() ; )
314 // it's the next single geometry
315 if ( prev_seg_id_ptr == NULL
316 || prev_seg_id_ptr->multi_index != it->operations[1].seg_id.multi_index )
318 // if the first ring has no IPs
319 if ( it->operations[1].seg_id.ring_index > -1 )
321 // we can be sure that the exterior overlaps the boundary
322 relate::set<exterior, boundary, '1', TransposeResult>(result);
325 // if there was some previous ring
326 if ( prev_seg_id_ptr != NULL )
328 signed_size_type const next_ring_index = prev_seg_id_ptr->ring_index + 1;
329 BOOST_GEOMETRY_ASSERT(next_ring_index >= 0);
331 // if one of the last rings of previous single geometry was ommited
332 if ( static_cast<std::size_t>(next_ring_index)
333 < geometry::num_interior_rings(
334 single_geometry(geometry2, *prev_seg_id_ptr)) )
336 // we can be sure that the exterior overlaps the boundary
337 relate::set<exterior, boundary, '1', TransposeResult>(result);
342 // if it's the same single geometry
343 else /*if ( previous_multi_index == it->operations[1].seg_id.multi_index )*/
345 // and we jumped over one of the rings
346 if ( prev_seg_id_ptr != NULL // just in case
347 && prev_seg_id_ptr->ring_index + 1 < it->operations[1].seg_id.ring_index )
349 // we can be sure that the exterior overlaps the boundary
350 relate::set<exterior, boundary, '1', TransposeResult>(result);
355 prev_seg_id_ptr = boost::addressof(it->operations[1].seg_id);
357 // find the next ring first iterator and check if the analysis should be performed
358 has_boundary_intersection has_boundary_inters;
359 turn_iterator next = find_next_ring(it, turns.end(), has_boundary_inters);
361 // if there is no 1d overlap with the boundary
362 if ( !has_boundary_inters.result )
364 // we can be sure that the exterior overlaps the boundary
365 relate::set<exterior, boundary, '1', TransposeResult>(result);
368 // else there is 1d overlap with the boundary so we must analyse the boundary
372 typedef turns::less<1, turns::less_op_areal_linear<1> > less;
373 std::sort(it, next, less());
376 areal_boundary_analyser<turn_type> analyser;
377 for ( turn_iterator rit = it ; rit != next ; ++rit )
379 // if the analyser requests, break the search
380 if ( !analyser.apply(it, rit, next) )
384 // if the boundary of Areal goes out of the Linear
385 if ( analyser.is_union_detected )
387 // we can be sure that the boundary of Areal overlaps the exterior of Linear
388 relate::set<exterior, boundary, '1', TransposeResult>(result);
396 // if there was some previous ring
397 if ( prev_seg_id_ptr != NULL )
399 signed_size_type const next_ring_index = prev_seg_id_ptr->ring_index + 1;
400 BOOST_GEOMETRY_ASSERT(next_ring_index >= 0);
402 // if one of the last rings of previous single geometry was ommited
403 if ( static_cast<std::size_t>(next_ring_index)
404 < geometry::num_interior_rings(
405 single_geometry(geometry2, *prev_seg_id_ptr)) )
407 // we can be sure that the exterior overlaps the boundary
408 relate::set<exterior, boundary, '1', TransposeResult>(result);
414 template <typename It, typename Pred, typename Comp>
415 static void for_each_equal_range(It first, It last, Pred pred, Comp comp)
420 It first_equal = first;
422 for ( ++first ; ; ++first, ++prev )
424 if ( first == last || !comp(*prev, *first) )
426 pred(first_equal, first);
437 template <typename Turn>
438 bool operator()(Turn const& left, Turn const& right) const
440 return left.operations[0].seg_id == right.operations[0].seg_id
441 && left.operations[0].fraction == right.operations[0].fraction;
445 struct same_ip_and_multi_index
447 template <typename Turn>
448 bool operator()(Turn const& left, Turn const& right) const
450 return same_ip()(left, right)
451 && left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index;
455 template <typename OpToPriority>
456 struct set_turns_group_priority
458 template <typename TurnIt>
459 void operator()(TurnIt first, TurnIt last) const
461 BOOST_GEOMETRY_ASSERT(first != last);
462 static OpToPriority op_to_priority;
463 // find the operation with the least priority
464 int least_priority = op_to_priority(first->operations[0]);
465 for ( TurnIt it = first + 1 ; it != last ; ++it )
467 int priority = op_to_priority(it->operations[0]);
468 if ( priority < least_priority )
469 least_priority = priority;
471 // set the least priority for all turns of the group
472 for ( TurnIt it = first ; it != last ; ++it )
474 it->priority = least_priority;
479 template <typename SingleLess>
480 struct sort_turns_group
484 template <typename Turn>
485 bool operator()(Turn const& left, Turn const& right) const
487 return left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index ?
488 SingleLess()(left, right) :
489 left.priority < right.priority;
493 template <typename TurnIt>
494 void operator()(TurnIt first, TurnIt last) const
496 std::sort(first, last, less());
500 template <typename TurnIt>
501 static void sort_dispatch(TurnIt first, TurnIt last, boost::true_type const& /*is_multi*/)
503 // sort turns by Linear seg_id, then by fraction, then by other multi_index
504 typedef turns::less<0, turns::less_other_multi_index<0> > less;
505 std::sort(first, last, less());
507 // For the same IP and multi_index - the same other's single geometry
508 // set priorities as the least operation found for the whole single geometry
509 // so e.g. single geometries containing 'u' will always be before those only containing 'i'
510 typedef turns::op_to_int<0,2,3,1,4,0> op_to_int_xuic;
511 for_each_equal_range(first, last,
512 set_turns_group_priority<op_to_int_xuic>(), // least operation in xuic order
513 same_ip_and_multi_index()); // other's multi_index
515 // When priorities for single geometries are set now sort turns for the same IP
516 // if multi_index is the same sort them according to the single-less
517 // else use priority of the whole single-geometry set earlier
518 typedef turns::less<0, turns::less_op_linear_areal_single<0> > single_less;
519 for_each_equal_range(first, last,
520 sort_turns_group<single_less>(),
524 template <typename TurnIt>
525 static void sort_dispatch(TurnIt first, TurnIt last, boost::false_type const& /*is_multi*/)
527 // sort turns by Linear seg_id, then by fraction, then
528 // for same ring id: x, u, i, c
529 // for different ring id: c, i, u, x
530 typedef turns::less<0, turns::less_op_linear_areal_single<0> > less;
531 std::sort(first, last, less());
535 // interrupt policy which may be passed to get_turns to interrupt the analysis
536 // based on the info in the passed result/mask
537 template <typename Areal, typename Result>
538 class interrupt_policy_linear_areal
541 static bool const enabled = true;
543 interrupt_policy_linear_areal(Areal const& areal, Result & result)
544 : m_result(result), m_areal(areal)
545 , is_boundary_found(false)
548 // TODO: since we update result for some operations here, we may not do it in the analyser!
550 template <typename Range>
551 inline bool apply(Range const& turns)
553 typedef typename boost::range_iterator<Range const>::type iterator;
555 for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
557 if ( it->operations[0].operation == overlay::operation_intersection )
559 bool const no_interior_rings
560 = geometry::num_interior_rings(
561 single_geometry(m_areal, it->operations[1].seg_id)) == 0;
563 // WARNING! THIS IS TRUE ONLY IF THE POLYGON IS SIMPLE!
564 // OR WITHOUT INTERIOR RINGS (AND OF COURSE VALID)
565 if ( no_interior_rings )
566 update<interior, interior, '1', TransposeResult>(m_result);
568 else if ( it->operations[0].operation == overlay::operation_continue )
570 update<interior, boundary, '1', TransposeResult>(m_result);
571 is_boundary_found = true;
573 else if ( ( it->operations[0].operation == overlay::operation_union
574 || it->operations[0].operation == overlay::operation_blocked )
575 && it->operations[0].position == overlay::position_middle )
577 // TODO: here we could also check the boundaries and set BB at this point
578 update<interior, boundary, '0', TransposeResult>(m_result);
582 return m_result.interrupt;
587 Areal const& m_areal;
590 bool is_boundary_found;
593 // This analyser should be used like Input or SinglePass Iterator
594 // IMPORTANT! It should be called also for the end iterator - last
595 template <typename TurnInfo>
598 typedef typename TurnInfo::point_type turn_point_type;
600 static const std::size_t op_id = 0;
601 static const std::size_t other_op_id = 1;
605 : m_previous_turn_ptr(NULL)
606 , m_previous_operation(overlay::operation_none)
607 , m_boundary_counter(0)
608 , m_interior_detected(false)
609 , m_first_interior_other_id_ptr(NULL)
610 , m_first_from_unknown(false)
611 , m_first_from_unknown_boundary_detected(false)
614 template <typename Result,
617 typename OtherGeometry,
618 typename BoundaryChecker>
619 void apply(Result & res, TurnIt it,
620 Geometry const& geometry,
621 OtherGeometry const& other_geometry,
622 BoundaryChecker const& boundary_checker)
624 overlay::operation_type op = it->operations[op_id].operation;
626 if ( op != overlay::operation_union
627 && op != overlay::operation_intersection
628 && op != overlay::operation_blocked
629 && op != overlay::operation_continue ) // operation_boundary / operation_boundary_intersection
634 segment_identifier const& seg_id = it->operations[op_id].seg_id;
635 segment_identifier const& other_id = it->operations[other_op_id].seg_id;
637 const bool first_in_range = m_seg_watcher.update(seg_id);
639 // TODO: should apply() for the post-last ip be called if first_in_range ?
640 // this would unify how last points in ranges are handled
641 // possibly replacing parts of the code below
642 // e.g. for is_multi and m_interior_detected
644 // handle possible exit
645 bool fake_enter_detected = false;
646 if ( m_exit_watcher.get_exit_operation() == overlay::operation_union )
648 // real exit point - may be multiple
649 // we know that we entered and now we exit
650 if ( ! turn_on_the_same_ip<op_id>(m_exit_watcher.get_exit_turn(), *it) )
652 m_exit_watcher.reset_detected_exit();
654 update<interior, exterior, '1', TransposeResult>(res);
656 // next single geometry
657 if ( first_in_range && m_previous_turn_ptr )
659 // NOTE: similar code is in the post-last-ip-apply()
660 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
662 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
663 range::back(sub_range(geometry, prev_seg_id)),
666 // if there is a boundary on the last point
669 update<boundary, exterior, '0', TransposeResult>(res);
673 // fake exit point, reset state
674 else if ( op == overlay::operation_intersection
675 || op == overlay::operation_continue ) // operation_boundary
677 m_exit_watcher.reset_detected_exit();
678 fake_enter_detected = true;
681 else if ( m_exit_watcher.get_exit_operation() == overlay::operation_blocked )
683 // ignore multiple BLOCKs for this same single geometry1
684 if ( op == overlay::operation_blocked
685 && seg_id.multi_index == m_previous_turn_ptr->operations[op_id].seg_id.multi_index )
690 if ( ( op == overlay::operation_intersection
691 || op == overlay::operation_continue )
692 && turn_on_the_same_ip<op_id>(m_exit_watcher.get_exit_turn(), *it) )
694 fake_enter_detected = true;
697 m_exit_watcher.reset_detected_exit();
700 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
701 && m_first_from_unknown )
703 // For MultiPolygon many x/u operations may be generated as a first IP
704 // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString
705 // then we know that the LineString is outside
706 // Similar with the u/u turns, if it was the first one it doesn't mean that the
707 // Linestring came from the exterior
708 if ( ( m_previous_operation == overlay::operation_blocked
709 && ( op != overlay::operation_blocked // operation different than block
710 || seg_id.multi_index != m_previous_turn_ptr->operations[op_id].seg_id.multi_index ) ) // or the next single-geometry
711 || ( m_previous_operation == overlay::operation_union
712 && ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
715 update<interior, exterior, '1', TransposeResult>(res);
716 if ( m_first_from_unknown_boundary_detected )
718 update<boundary, exterior, '0', TransposeResult>(res);
721 m_first_from_unknown = false;
722 m_first_from_unknown_boundary_detected = false;
726 // NOTE: THE WHOLE m_interior_detected HANDLING IS HERE BECAUSE WE CAN'T EFFICIENTLY SORT TURNS (CORRECTLY)
727 // BECAUSE THE SAME IP MAY BE REPRESENTED BY TWO SEGMENTS WITH DIFFERENT DISTANCES
728 // IT WOULD REQUIRE THE CALCULATION OF MAX DISTANCE
729 // TODO: WE COULD GET RID OF THE TEST IF THE DISTANCES WERE NORMALIZED
731 // UPDATE: THEY SHOULD BE NORMALIZED NOW
733 // TODO: THIS IS POTENTIALLY ERROREOUS!
734 // THIS ALGORITHM DEPENDS ON SOME SPECIFIC SEQUENCE OF OPERATIONS
735 // IT WOULD GIVE WRONG RESULTS E.G.
736 // IN THE CASE OF SELF-TOUCHING POINT WHEN 'i' WOULD BE BEFORE 'u'
738 // handle the interior overlap
739 if ( m_interior_detected )
741 BOOST_GEOMETRY_ASSERT_MSG(m_previous_turn_ptr, "non-NULL ptr expected");
743 // real interior overlap
744 if ( ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
746 update<interior, interior, '1', TransposeResult>(res);
747 m_interior_detected = false;
749 // new range detected - reset previous state and check the boundary
750 if ( first_in_range )
752 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
754 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
755 range::back(sub_range(geometry, prev_seg_id)),
758 // if there is a boundary on the last point
761 update<boundary, interior, '0', TransposeResult>(res);
764 // The exit_watcher is reset below
765 // m_exit_watcher.reset();
768 // fake interior overlap
769 else if ( op == overlay::operation_continue )
771 m_interior_detected = false;
773 else if ( op == overlay::operation_union )
775 // TODO: this probably is not a good way of handling the interiors/enters
776 // the solution similar to exit_watcher would be more robust
777 // all enters should be kept and handled.
778 // maybe integrate it with the exit_watcher -> enter_exit_watcher
779 if ( m_first_interior_other_id_ptr
780 && m_first_interior_other_id_ptr->multi_index == other_id.multi_index )
782 m_interior_detected = false;
787 // NOTE: If post-last-ip apply() was called this wouldn't be needed
788 if ( first_in_range )
790 m_exit_watcher.reset();
791 m_boundary_counter = 0;
792 m_first_from_unknown = false;
793 m_first_from_unknown_boundary_detected = false;
797 if ( op == overlay::operation_intersection
798 || op == overlay::operation_continue ) // operation_boundary/operation_boundary_intersection
800 bool const first_point = first_in_range || m_first_from_unknown;
801 bool no_enters_detected = m_exit_watcher.is_outside();
802 m_exit_watcher.enter(*it);
804 if ( op == overlay::operation_intersection )
806 if ( m_boundary_counter > 0 && it->operations[op_id].is_collinear )
807 --m_boundary_counter;
809 if ( m_boundary_counter == 0 )
811 // interiors overlaps
812 //update<interior, interior, '1', TransposeResult>(res);
814 // TODO: think about the implementation of the more robust version
815 // this way only the first enter will be handled
816 if ( !m_interior_detected )
819 // we might enter a boundary of some other ring on the same IP
820 m_interior_detected = true;
821 m_first_interior_other_id_ptr = boost::addressof(other_id);
825 else // operation_boundary
827 // don't add to the count for all met boundaries
828 // only if this is the "new" boundary
829 if ( first_point || !it->operations[op_id].is_collinear )
830 ++m_boundary_counter;
832 update<interior, boundary, '1', TransposeResult>(res);
836 = is_ip_on_boundary<boundary_front>(it->point,
837 it->operations[op_id],
840 // going inside on boundary point
843 update<boundary, boundary, '0', TransposeResult>(res);
845 // going inside on non-boundary point
848 update<interior, boundary, '0', TransposeResult>(res);
850 // if we didn't enter in the past, we were outside
851 if ( no_enters_detected
852 && ! fake_enter_detected
853 && it->operations[op_id].position != overlay::position_front )
855 // TODO: calculate_from_inside() is only needed if the current Linestring is not closed
856 bool const from_inside = first_point
857 && calculate_from_inside(geometry,
862 update<interior, interior, '1', TransposeResult>(res);
864 update<interior, exterior, '1', TransposeResult>(res);
866 // if it's the first IP then the first point is outside
869 bool const front_b = is_endpoint_on_boundary<boundary_front>(
870 range::front(sub_range(geometry, seg_id)),
873 // if there is a boundary on the first point
877 update<boundary, interior, '0', TransposeResult>(res);
879 update<boundary, exterior, '0', TransposeResult>(res);
885 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) )
887 m_first_from_unknown = false;
888 m_first_from_unknown_boundary_detected = false;
892 else if ( op == overlay::operation_union || op == overlay::operation_blocked )
894 bool const op_blocked = op == overlay::operation_blocked;
895 bool const no_enters_detected = m_exit_watcher.is_outside()
896 // TODO: is this condition ok?
897 // TODO: move it into the exit_watcher?
898 && m_exit_watcher.get_exit_operation() == overlay::operation_none;
900 if ( op == overlay::operation_union )
902 if ( m_boundary_counter > 0 && it->operations[op_id].is_collinear )
903 --m_boundary_counter;
905 else // overlay::operation_blocked
907 m_boundary_counter = 0;
910 // we're inside, possibly going out right now
911 if ( ! no_enters_detected )
914 && it->operations[op_id].position == overlay::position_back ) // ignore spikes!
916 // check if this is indeed the boundary point
917 // NOTE: is_ip_on_boundary<>() should be called here but the result will be the same
918 if ( is_endpoint_on_boundary<boundary_back>(it->point, boundary_checker) )
920 update<boundary, boundary, '0', TransposeResult>(res);
923 // union, inside, but no exit -> collinear on self-intersection point
924 // not needed since we're already inside the boundary
925 /*else if ( !exit_detected )
927 update<interior, boundary, '0', TransposeResult>(res);
930 // we're outside or inside and this is the first turn
933 bool const this_b = is_ip_on_boundary<boundary_any>(it->point,
934 it->operations[op_id],
937 // if current IP is on boundary of the geometry
940 update<boundary, boundary, '0', TransposeResult>(res);
942 // if current IP is not on boundary of the geometry
945 update<interior, boundary, '0', TransposeResult>(res);
948 // TODO: very similar code is used in the handling of intersection
949 if ( it->operations[op_id].position != overlay::position_front )
951 // TODO: calculate_from_inside() is only needed if the current Linestring is not closed
952 // NOTE: this is not enough for MultiPolygon and operation_blocked
953 // For LS/MultiPolygon multiple x/u turns may be generated
954 // the first checked Polygon may be the one which LS is outside for.
955 bool const first_point = first_in_range || m_first_from_unknown;
956 bool const first_from_inside = first_point
957 && calculate_from_inside(geometry,
960 if ( first_from_inside )
962 update<interior, interior, '1', TransposeResult>(res);
964 // notify the exit_watcher that we started inside
965 m_exit_watcher.enter(*it);
966 // and reset unknown flags since we know that we started inside
967 m_first_from_unknown = false;
968 m_first_from_unknown_boundary_detected = false;
972 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
973 /*&& ( op == overlay::operation_blocked
974 || op == overlay::operation_union )*/ ) // if we're here it's u or x
976 m_first_from_unknown = true;
980 update<interior, exterior, '1', TransposeResult>(res);
984 // first IP on the last segment point - this means that the first point is outside or inside
985 if ( first_point && ( !this_b || op_blocked ) )
987 bool const front_b = is_endpoint_on_boundary<boundary_front>(
988 range::front(sub_range(geometry, seg_id)),
991 // if there is a boundary on the first point
994 if ( first_from_inside )
996 update<boundary, interior, '0', TransposeResult>(res);
1000 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
1001 /*&& ( op == overlay::operation_blocked
1002 || op == overlay::operation_union )*/ ) // if we're here it's u or x
1004 BOOST_GEOMETRY_ASSERT(m_first_from_unknown);
1005 m_first_from_unknown_boundary_detected = true;
1009 update<boundary, exterior, '0', TransposeResult>(res);
1017 // if we're going along a boundary, we exit only if the linestring was collinear
1018 if ( m_boundary_counter == 0
1019 || it->operations[op_id].is_collinear )
1021 // notify the exit watcher about the possible exit
1022 m_exit_watcher.exit(*it);
1026 // store ref to previously analysed (valid) turn
1027 m_previous_turn_ptr = boost::addressof(*it);
1028 // and previously analysed (valid) operation
1029 m_previous_operation = op;
1033 template <typename Result,
1036 typename OtherGeometry,
1037 typename BoundaryChecker>
1038 void apply(Result & res,
1039 TurnIt first, TurnIt last,
1040 Geometry const& geometry,
1041 OtherGeometry const& /*other_geometry*/,
1042 BoundaryChecker const& boundary_checker)
1044 boost::ignore_unused(first, last);
1045 //BOOST_GEOMETRY_ASSERT( first != last );
1047 // For MultiPolygon many x/u operations may be generated as a first IP
1048 // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString
1049 // then we know that the LineString is outside
1050 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
1051 && m_first_from_unknown )
1053 update<interior, exterior, '1', TransposeResult>(res);
1054 if ( m_first_from_unknown_boundary_detected )
1056 update<boundary, exterior, '0', TransposeResult>(res);
1060 //m_first_from_unknown = false;
1061 //m_first_from_unknown_boundary_detected = false;
1064 // here, the possible exit is the real one
1065 // we know that we entered and now we exit
1066 if ( /*m_exit_watcher.get_exit_operation() == overlay::operation_union // THIS CHECK IS REDUNDANT
1067 ||*/ m_previous_operation == overlay::operation_union
1068 && !m_interior_detected )
1071 update<interior, exterior, '1', TransposeResult>(res);
1073 BOOST_GEOMETRY_ASSERT(first != last);
1074 BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr);
1076 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1078 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1079 range::back(sub_range(geometry, prev_seg_id)),
1082 // if there is a boundary on the last point
1085 update<boundary, exterior, '0', TransposeResult>(res);
1088 // we might enter some Areal and didn't go out,
1089 else if ( m_previous_operation == overlay::operation_intersection
1090 || m_interior_detected )
1093 update<interior, interior, '1', TransposeResult>(res);
1094 m_interior_detected = false;
1096 BOOST_GEOMETRY_ASSERT(first != last);
1097 BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr);
1099 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1101 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1102 range::back(sub_range(geometry, prev_seg_id)),
1105 // if there is a boundary on the last point
1108 update<boundary, interior, '0', TransposeResult>(res);
1112 // This condition may be false if the Linestring is lying on the Polygon's collinear spike
1113 // if Polygon's spikes are not handled in get_turns() or relate() (they currently aren't)
1114 //BOOST_GEOMETRY_ASSERT_MSG(m_previous_operation != overlay::operation_continue,
1115 // "Unexpected operation! Probably the error in get_turns(L,A) or relate(L,A)");
1116 // Currently one c/c turn is generated for the exit
1117 // when a Linestring is going out on a collinear spike
1118 // When a Linestring is going in on a collinear spike
1119 // the turn is not generated for the entry
1120 // So assume it's the former
1121 if ( m_previous_operation == overlay::operation_continue )
1123 update<interior, exterior, '1', TransposeResult>(res);
1125 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1127 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1128 range::back(sub_range(geometry, prev_seg_id)),
1131 // if there is a boundary on the last point
1134 update<boundary, exterior, '0', TransposeResult>(res);
1138 // Reset exit watcher before the analysis of the next Linestring
1139 m_exit_watcher.reset();
1140 m_boundary_counter = 0;
1141 m_first_from_unknown = false;
1142 m_first_from_unknown_boundary_detected = false;
1145 // check if the passed turn's segment of Linear geometry arrived
1146 // from the inside or the outside of the Areal geometry
1147 template <typename Turn>
1148 static inline bool calculate_from_inside(Geometry1 const& geometry1,
1149 Geometry2 const& geometry2,
1152 typedef typename cs_tag<typename Turn::point_type>::type cs_tag;
1154 if ( turn.operations[op_id].position == overlay::position_front )
1157 typename sub_range_return_type<Geometry1 const>::type
1158 range1 = sub_range(geometry1, turn.operations[op_id].seg_id);
1160 typedef detail::normalized_view<Geometry2 const> const range2_type;
1161 typedef typename boost::range_iterator<range2_type>::type range2_iterator;
1162 range2_type range2(sub_range(geometry2, turn.operations[other_op_id].seg_id));
1164 BOOST_GEOMETRY_ASSERT(boost::size(range1));
1165 std::size_t const s2 = boost::size(range2);
1166 BOOST_GEOMETRY_ASSERT(s2 > 2);
1167 std::size_t const seg_count2 = s2 - 1;
1169 std::size_t const p_seg_ij = static_cast<std::size_t>(turn.operations[op_id].seg_id.segment_index);
1170 std::size_t const q_seg_ij = static_cast<std::size_t>(turn.operations[other_op_id].seg_id.segment_index);
1172 BOOST_GEOMETRY_ASSERT(p_seg_ij + 1 < boost::size(range1));
1173 BOOST_GEOMETRY_ASSERT(q_seg_ij + 1 < s2);
1175 point1_type const& pi = range::at(range1, p_seg_ij);
1176 point2_type const& qi = range::at(range2, q_seg_ij);
1177 point2_type const& qj = range::at(range2, q_seg_ij + 1);
1178 point1_type qi_conv;
1179 geometry::convert(qi, qi_conv);
1180 bool const is_ip_qj = equals::equals_point_point(turn.point, qj);
1182 // BOOST_GEOMETRY_ASSERT(!equals::equals_point_point(turn.point, pi));
1183 // BOOST_GEOMETRY_ASSERT(!equals::equals_point_point(turn.point, qi));
1185 geometry::convert(turn.point, new_pj);
1189 std::size_t const q_seg_jk = (q_seg_ij + 1) % seg_count2;
1190 // TODO: the following function should return immediately, however the worst case complexity is O(N)
1191 // It would be good to replace it with some O(1) mechanism
1192 range2_iterator qk_it = find_next_non_duplicated(boost::begin(range2),
1193 range::pos(range2, q_seg_jk),
1194 boost::end(range2));
1196 // Will this sequence of points be always correct?
1197 overlay::side_calculator<cs_tag, point1_type, point2_type> side_calc(qi_conv, new_pj, pi, qi, qj, *qk_it);
1199 return calculate_from_inside_sides(side_calc);
1204 geometry::convert(turn.point, new_qj);
1206 overlay::side_calculator<cs_tag, point1_type, point2_type> side_calc(qi_conv, new_pj, pi, qi, new_qj, qj);
1208 return calculate_from_inside_sides(side_calc);
1212 template <typename It>
1213 static inline It find_next_non_duplicated(It first, It current, It last)
1215 BOOST_GEOMETRY_ASSERT( current != last );
1219 for ( ++it ; it != last ; ++it )
1221 if ( !equals::equals_point_point(*current, *it) )
1225 // if not found start from the beginning
1226 for ( it = first ; it != current ; ++it )
1228 if ( !equals::equals_point_point(*current, *it) )
1235 // calculate inside or outside based on side_calc
1236 // this is simplified version of a check from equal<>
1237 template <typename SideCalc>
1238 static inline bool calculate_from_inside_sides(SideCalc const& side_calc)
1240 int const side_pk_p = side_calc.pk_wrt_p1();
1241 int const side_qk_p = side_calc.qk_wrt_p1();
1242 // If they turn to same side (not opposite sides)
1243 if (! overlay::base_turn_handler::opposite(side_pk_p, side_qk_p))
1245 int const side_pk_q2 = side_calc.pk_wrt_q2();
1246 return side_pk_q2 == -1;
1250 return side_pk_p == -1;
1255 exit_watcher<TurnInfo, op_id> m_exit_watcher;
1256 segment_watcher<same_single> m_seg_watcher;
1257 TurnInfo * m_previous_turn_ptr;
1258 overlay::operation_type m_previous_operation;
1259 unsigned m_boundary_counter;
1260 bool m_interior_detected;
1261 const segment_identifier * m_first_interior_other_id_ptr;
1262 bool m_first_from_unknown;
1263 bool m_first_from_unknown_boundary_detected;
1266 // call analyser.apply() for each turn in range
1267 // IMPORTANT! The analyser is also called for the end iterator - last
1268 template <typename Result,
1272 typename OtherGeometry,
1273 typename BoundaryChecker>
1274 static inline void analyse_each_turn(Result & res,
1275 Analyser & analyser,
1276 TurnIt first, TurnIt last,
1277 Geometry const& geometry,
1278 OtherGeometry const& other_geometry,
1279 BoundaryChecker const& boundary_checker)
1281 if ( first == last )
1284 for ( TurnIt it = first ; it != last ; ++it )
1286 analyser.apply(res, it,
1287 geometry, other_geometry,
1290 if ( BOOST_GEOMETRY_CONDITION( res.interrupt ) )
1294 analyser.apply(res, first, last,
1295 geometry, other_geometry,
1299 // less comparator comparing multi_index then ring_index
1300 // may be used to sort turns by ring
1303 template <typename Turn>
1304 inline bool operator()(Turn const& left, Turn const& right) const
1306 return left.operations[1].seg_id.multi_index < right.operations[1].seg_id.multi_index
1307 || ( left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index
1308 && left.operations[1].seg_id.ring_index < right.operations[1].seg_id.ring_index );
1312 // policy/functor checking if a turn's operation is operation_continue
1313 struct has_boundary_intersection
1315 has_boundary_intersection()
1318 template <typename Turn>
1319 inline void operator()(Turn const& turn)
1321 if ( turn.operations[1].operation == overlay::operation_continue )
1328 // iterate through the range and search for the different multi_index or ring_index
1329 // also call fun for each turn in the current range
1330 template <typename It, typename Fun>
1331 static inline It find_next_ring(It first, It last, Fun & fun)
1333 if ( first == last )
1336 signed_size_type const multi_index = first->operations[1].seg_id.multi_index;
1337 signed_size_type const ring_index = first->operations[1].seg_id.ring_index;
1342 for ( ; first != last ; ++first )
1344 if ( multi_index != first->operations[1].seg_id.multi_index
1345 || ring_index != first->operations[1].seg_id.ring_index )
1356 // analyser which called for turns sorted by seg/distance/operation
1357 // checks if the boundary of Areal geometry really got out
1358 // into the exterior of Linear geometry
1359 template <typename TurnInfo>
1360 class areal_boundary_analyser
1363 areal_boundary_analyser()
1364 : is_union_detected(false)
1365 , m_previous_turn_ptr(NULL)
1368 template <typename TurnIt>
1369 bool apply(TurnIt /*first*/, TurnIt it, TurnIt last)
1371 overlay::operation_type op = it->operations[1].operation;
1375 if ( op != overlay::operation_union
1376 && op != overlay::operation_continue )
1381 if ( is_union_detected )
1383 BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr != NULL);
1384 if ( !detail::equals::equals_point_point(it->point, m_previous_turn_ptr->point) )
1389 else if ( op == overlay::operation_continue ) // operation_boundary
1391 is_union_detected = false;
1395 if ( op == overlay::operation_union )
1397 is_union_detected = true;
1398 m_previous_turn_ptr = boost::addressof(*it);
1409 bool is_union_detected;
1412 const TurnInfo * m_previous_turn_ptr;
1416 template <typename Geometry1, typename Geometry2>
1419 typedef linear_areal<Geometry2, Geometry1, true> linear_areal_type;
1421 static const bool interruption_enabled = linear_areal_type::interruption_enabled;
1423 template <typename Result>
1424 static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result)
1426 linear_areal_type::apply(geometry2, geometry1, result);
1430 }} // namespace detail::relate
1431 #endif // DOXYGEN_NO_DETAIL
1433 }} // namespace boost::geometry
1435 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP