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, 2017, 2018, 2019.
6 // Modifications copyright (c) 2013-2019 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_AREAL_AREAL_HPP
15 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP
17 #include <boost/geometry/core/topological_dimension.hpp>
19 #include <boost/geometry/util/condition.hpp>
20 #include <boost/geometry/util/range.hpp>
22 #include <boost/geometry/algorithms/num_interior_rings.hpp>
23 #include <boost/geometry/algorithms/detail/point_on_border.hpp>
24 #include <boost/geometry/algorithms/detail/sub_range.hpp>
25 #include <boost/geometry/algorithms/detail/single_geometry.hpp>
27 #include <boost/geometry/algorithms/detail/relate/point_geometry.hpp>
28 #include <boost/geometry/algorithms/detail/relate/turns.hpp>
29 #include <boost/geometry/algorithms/detail/relate/boundary_checker.hpp>
30 #include <boost/geometry/algorithms/detail/relate/follow_helpers.hpp>
32 namespace boost { namespace geometry
35 #ifndef DOXYGEN_NO_DETAIL
36 namespace detail { namespace relate {
39 // TODO: In the worst case calling this Pred in a loop for MultiPolygon/MultiPolygon may take O(NM)
40 // Use the rtree in this case!
42 // may be used to set EI and EB for an Areal geometry for which no turns were generated
47 typename PointInArealStrategy,
50 class no_turns_aa_pred
53 no_turns_aa_pred(OtherAreal const& other_areal,
55 PointInArealStrategy const& point_in_areal_strategy)
57 , m_point_in_areal_strategy(point_in_areal_strategy)
58 , m_other_areal(other_areal)
61 // check which relations must be analysed
63 if ( ! may_update<interior, interior, '2', TransposeResult>(m_result)
64 && ! may_update<boundary, interior, '1', TransposeResult>(m_result)
65 && ! may_update<exterior, interior, '2', TransposeResult>(m_result) )
70 if ( ! may_update<interior, exterior, '2', TransposeResult>(m_result)
71 && ! may_update<boundary, exterior, '1', TransposeResult>(m_result) )
77 template <typename Areal>
78 bool operator()(Areal const& areal)
80 using detail::within::point_in_geometry;
82 // if those flags are set nothing will change
88 typedef typename geometry::point_type<Areal>::type point_type;
90 bool const ok = boost::geometry::point_on_border(pt, areal);
92 // TODO: for now ignore, later throw an exception?
98 // check if the areal is inside the other_areal
100 // Run in a loop O(NM) - optimize!
101 int const pig = point_in_geometry(pt,
103 m_point_in_areal_strategy);
104 //BOOST_GEOMETRY_ASSERT( pig != 0 );
109 update<interior, interior, '2', TransposeResult>(m_result);
110 update<boundary, interior, '1', TransposeResult>(m_result);
111 update<exterior, interior, '2', TransposeResult>(m_result);
115 // Only the interior rings of other ONE single geometry must be checked
116 // NOT all geometries
118 // Check if any interior ring is outside
119 ring_identifier ring_id(0, -1, 0);
120 std::size_t const irings_count = geometry::num_interior_rings(areal);
121 for ( ; static_cast<std::size_t>(ring_id.ring_index) < irings_count ;
122 ++ring_id.ring_index )
124 typename detail::sub_range_return_type<Areal const>::type
125 range_ref = detail::sub_range(areal, ring_id);
127 if ( boost::empty(range_ref) )
129 // TODO: throw exception?
135 int const hpig = point_in_geometry(range::front(range_ref),
137 m_point_in_areal_strategy);
142 update<interior, exterior, '2', TransposeResult>(m_result);
143 update<boundary, exterior, '1', TransposeResult>(m_result);
152 update<interior, exterior, '2', TransposeResult>(m_result);
153 update<boundary, exterior, '1', TransposeResult>(m_result);
156 // Check if any interior ring is inside
157 ring_identifier ring_id(0, -1, 0);
158 std::size_t const irings_count = geometry::num_interior_rings(areal);
159 for ( ; static_cast<std::size_t>(ring_id.ring_index) < irings_count ;
160 ++ring_id.ring_index )
162 typename detail::sub_range_return_type<Areal const>::type
163 range_ref = detail::sub_range(areal, ring_id);
165 if ( boost::empty(range_ref) )
167 // TODO: throw exception?
173 int const hpig = point_in_geometry(range::front(range_ref),
175 m_point_in_areal_strategy);
180 update<interior, interior, '2', TransposeResult>(m_result);
181 update<boundary, interior, '1', TransposeResult>(m_result);
182 update<exterior, interior, '2', TransposeResult>(m_result);
189 return m_flags != 3 && !m_result.interrupt;
194 PointInArealStrategy const& m_point_in_areal_strategy;
195 OtherAreal const& m_other_areal;
199 // The implementation of an algorithm calculating relate() for A/A
200 template <typename Geometry1, typename Geometry2>
203 // check Linear / Areal
204 BOOST_STATIC_ASSERT(topological_dimension<Geometry1>::value == 2
205 && topological_dimension<Geometry2>::value == 2);
207 static const bool interruption_enabled = true;
209 typedef typename geometry::point_type<Geometry1>::type point1_type;
210 typedef typename geometry::point_type<Geometry2>::type point2_type;
212 template <typename Result, typename IntersectionStrategy>
213 static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2,
215 IntersectionStrategy const& intersection_strategy)
217 // TODO: If Areal geometry may have infinite size, change the following line:
219 // The result should be FFFFFFFFF
220 relate::set<exterior, exterior, result_dimension<Geometry2>::value>(result);// FFFFFFFFd, d in [1,9] or T
222 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
225 // get and analyse turns
226 typedef typename turns::get_turns
229 >::template turn_info_type<IntersectionStrategy>::type turn_type;
230 std::vector<turn_type> turns;
232 interrupt_policy_areal_areal<Result> interrupt_policy(geometry1, geometry2, result);
234 turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy, intersection_strategy);
235 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
238 typedef typename IntersectionStrategy::cs_tag cs_tag;
240 typedef typename IntersectionStrategy::template point_in_geometry_strategy
243 >::type point_in_areal_strategy12_type;
244 point_in_areal_strategy12_type point_in_areal_strategy12
245 = intersection_strategy.template get_point_in_geometry_strategy<Geometry1, Geometry2>();
246 typedef typename IntersectionStrategy::template point_in_geometry_strategy
249 >::type point_in_areal_strategy21_type;
250 point_in_areal_strategy21_type point_in_areal_strategy21
251 = intersection_strategy.template get_point_in_geometry_strategy<Geometry2, Geometry1>();
253 no_turns_aa_pred<Geometry2, Result, point_in_areal_strategy12_type, false>
254 pred1(geometry2, result, point_in_areal_strategy12);
255 for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
256 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
259 no_turns_aa_pred<Geometry1, Result, point_in_areal_strategy21_type, true>
260 pred2(geometry1, result, point_in_areal_strategy21);
261 for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
262 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
268 if ( may_update<interior, interior, '2'>(result)
269 || may_update<interior, exterior, '2'>(result)
270 || may_update<boundary, interior, '1'>(result)
271 || may_update<boundary, exterior, '1'>(result)
272 || may_update<exterior, interior, '2'>(result) )
275 typedef turns::less<0, turns::less_op_areal_areal<0>, cs_tag> less;
276 std::sort(turns.begin(), turns.end(), less());
278 /*if ( may_update<interior, exterior, '2'>(result)
279 || may_update<boundary, exterior, '1'>(result)
280 || may_update<boundary, interior, '1'>(result)
281 || may_update<exterior, interior, '2'>(result) )*/
283 // analyse sorted turns
284 turns_analyser<turn_type, 0> analyser;
285 analyse_each_turn(result, analyser, turns.begin(), turns.end(),
286 point_in_areal_strategy12.get_equals_point_point_strategy());
288 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
292 if ( may_update<interior, interior, '2'>(result)
293 || may_update<interior, exterior, '2'>(result)
294 || may_update<boundary, interior, '1'>(result)
295 || may_update<boundary, exterior, '1'>(result)
296 || may_update<exterior, interior, '2'>(result) )
298 // analyse rings for which turns were not generated
299 // or only i/i or u/u was generated
300 uncertain_rings_analyser<0, Result, Geometry1, Geometry2, point_in_areal_strategy12_type>
301 rings_analyser(result, geometry1, geometry2, point_in_areal_strategy12);
302 analyse_uncertain_rings<0>::apply(rings_analyser, turns.begin(), turns.end());
304 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
309 if ( may_update<interior, interior, '2', true>(result)
310 || may_update<interior, exterior, '2', true>(result)
311 || may_update<boundary, interior, '1', true>(result)
312 || may_update<boundary, exterior, '1', true>(result)
313 || may_update<exterior, interior, '2', true>(result) )
316 typedef turns::less<1, turns::less_op_areal_areal<1>, cs_tag> less;
317 std::sort(turns.begin(), turns.end(), less());
319 /*if ( may_update<interior, exterior, '2', true>(result)
320 || may_update<boundary, exterior, '1', true>(result)
321 || may_update<boundary, interior, '1', true>(result)
322 || may_update<exterior, interior, '2', true>(result) )*/
324 // analyse sorted turns
325 turns_analyser<turn_type, 1> analyser;
326 analyse_each_turn(result, analyser, turns.begin(), turns.end(),
327 point_in_areal_strategy21.get_equals_point_point_strategy());
329 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
333 if ( may_update<interior, interior, '2', true>(result)
334 || may_update<interior, exterior, '2', true>(result)
335 || may_update<boundary, interior, '1', true>(result)
336 || may_update<boundary, exterior, '1', true>(result)
337 || may_update<exterior, interior, '2', true>(result) )
339 // analyse rings for which turns were not generated
340 // or only i/i or u/u was generated
341 uncertain_rings_analyser<1, Result, Geometry2, Geometry1, point_in_areal_strategy21_type>
342 rings_analyser(result, geometry2, geometry1, point_in_areal_strategy21);
343 analyse_uncertain_rings<1>::apply(rings_analyser, turns.begin(), turns.end());
345 //if ( result.interrupt )
351 // interrupt policy which may be passed to get_turns to interrupt the analysis
352 // based on the info in the passed result/mask
353 template <typename Result>
354 class interrupt_policy_areal_areal
357 static bool const enabled = true;
359 interrupt_policy_areal_areal(Geometry1 const& geometry1,
360 Geometry2 const& geometry2,
363 , m_geometry1(geometry1)
364 , m_geometry2(geometry2)
367 template <typename Range>
368 inline bool apply(Range const& turns)
370 typedef typename boost::range_iterator<Range const>::type iterator;
372 for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
378 return m_result.interrupt;
382 template <std::size_t OpId, typename Turn>
383 inline void per_turn(Turn const& turn)
385 //static const std::size_t other_op_id = (OpId + 1) % 2;
386 static const bool transpose_result = OpId != 0;
388 overlay::operation_type const op = turn.operations[OpId].operation;
390 if ( op == overlay::operation_union )
393 /*if ( turn.operations[other_op_id].operation != overlay::operation_union )
395 update<interior, exterior, '2', transpose_result>(m_result);
396 update<boundary, exterior, '1', transpose_result>(m_result);
399 update<boundary, boundary, '0', transpose_result>(m_result);
401 else if ( op == overlay::operation_intersection )
404 /*if ( turn.operations[other_op_id].operation != overlay::operation_intersection )
406 // not correct e.g. for G1 touching G2 in a point where a hole is touching the exterior ring
407 // in this case 2 turns i/... and u/u will be generated for this IP
408 //update<interior, interior, '2', transpose_result>(m_result);
410 //update<boundary, interior, '1', transpose_result>(m_result);
413 update<boundary, boundary, '0', transpose_result>(m_result);
415 else if ( op == overlay::operation_continue )
417 update<boundary, boundary, '1', transpose_result>(m_result);
418 update<interior, interior, '2', transpose_result>(m_result);
420 else if ( op == overlay::operation_blocked )
422 update<boundary, boundary, '1', transpose_result>(m_result);
423 update<interior, exterior, '2', transpose_result>(m_result);
428 Geometry1 const& m_geometry1;
429 Geometry2 const& m_geometry2;
432 // This analyser should be used like Input or SinglePass Iterator
433 // IMPORTANT! It should be called also for the end iterator - last
434 template <typename TurnInfo, std::size_t OpId>
437 typedef typename TurnInfo::point_type turn_point_type;
439 static const std::size_t op_id = OpId;
440 static const std::size_t other_op_id = (OpId + 1) % 2;
441 static const bool transpose_result = OpId != 0;
445 : m_previous_turn_ptr(0)
446 , m_previous_operation(overlay::operation_none)
447 , m_enter_detected(false)
448 , m_exit_detected(false)
451 template <typename Result, typename TurnIt, typename EqPPStrategy>
452 void apply(Result & result, TurnIt it, EqPPStrategy const& strategy)
454 //BOOST_GEOMETRY_ASSERT( it != last );
456 overlay::operation_type const op = it->operations[op_id].operation;
458 if ( op != overlay::operation_union
459 && op != overlay::operation_intersection
460 && op != overlay::operation_blocked
461 && op != overlay::operation_continue )
466 segment_identifier const& seg_id = it->operations[op_id].seg_id;
467 //segment_identifier const& other_id = it->operations[other_op_id].seg_id;
469 const bool first_in_range = m_seg_watcher.update(seg_id);
471 if ( m_previous_turn_ptr )
473 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
475 // real exit point - may be multiple
477 || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it, strategy) )
480 m_exit_detected = false;
482 // fake exit point, reset state
483 else if ( op != overlay::operation_union )
485 m_exit_detected = false;
489 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
493 || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it, strategy) )
495 update_enter(result);
496 m_enter_detected = false;
498 // fake entry point, reset state
499 else if ( op != overlay::operation_intersection )
501 m_enter_detected = false;
506 if ( op == overlay::operation_union )
508 // already set in interrupt policy
509 //update<boundary, boundary, '0', transpose_result>(m_result);
512 //if ( it->operations[other_op_id].operation != overlay::operation_union )
514 m_exit_detected = true;
517 else if ( op == overlay::operation_intersection )
520 if ( it->operations[other_op_id].operation != overlay::operation_intersection )
522 // this was set in the interrupt policy but it was wrong
523 // also here it's wrong since it may be a fake entry point
524 //update<interior, interior, '2', transpose_result>(result);
526 // already set in interrupt policy
527 //update<boundary, boundary, '0', transpose_result>(result);
528 m_enter_detected = true;
531 else if ( op == overlay::operation_blocked )
533 // already set in interrupt policy
535 else // if ( op == overlay::operation_continue )
537 // already set in interrupt policy
540 // store ref to previously analysed (valid) turn
541 m_previous_turn_ptr = boost::addressof(*it);
542 // and previously analysed (valid) operation
543 m_previous_operation = op;
547 template <typename Result>
548 void apply(Result & result)
550 //BOOST_GEOMETRY_ASSERT( first != last );
552 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
555 m_exit_detected = false;
558 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
560 update_enter(result);
561 m_enter_detected = false;
565 template <typename Result>
566 static inline void update_exit(Result & result)
568 update<interior, exterior, '2', transpose_result>(result);
569 update<boundary, exterior, '1', transpose_result>(result);
572 template <typename Result>
573 static inline void update_enter(Result & result)
575 update<interior, interior, '2', transpose_result>(result);
576 update<boundary, interior, '1', transpose_result>(result);
577 update<exterior, interior, '2', transpose_result>(result);
581 segment_watcher<same_ring> m_seg_watcher;
582 TurnInfo * m_previous_turn_ptr;
583 overlay::operation_type m_previous_operation;
584 bool m_enter_detected;
585 bool m_exit_detected;
588 // call analyser.apply() for each turn in range
589 // IMPORTANT! The analyser is also called for the end iterator - last
595 typename EqPPStrategy
597 static inline void analyse_each_turn(Result & res,
599 TurnIt first, TurnIt last,
600 EqPPStrategy const& strategy)
605 for ( TurnIt it = first ; it != last ; ++it )
607 analyser.apply(res, it, strategy);
609 if ( BOOST_GEOMETRY_CONDITION(res.interrupt) )
621 typename OtherGeometry,
622 typename PointInArealStrategy
624 class uncertain_rings_analyser
626 static const bool transpose_result = OpId != 0;
627 static const int other_id = (OpId + 1) % 2;
630 inline uncertain_rings_analyser(Result & result,
631 Geometry const& geom,
632 OtherGeometry const& other_geom,
633 PointInArealStrategy const& point_in_areal_strategy)
635 , other_geometry(other_geom)
636 , interrupt(result.interrupt) // just in case, could be false as well
638 , m_point_in_areal_strategy(point_in_areal_strategy)
641 // check which relations must be analysed
642 // NOTE: 1 and 4 could probably be connected
644 if ( ! may_update<interior, interior, '2', transpose_result>(m_result) )
649 if ( ! may_update<interior, exterior, '2', transpose_result>(m_result)
650 && ! may_update<boundary, exterior, '1', transpose_result>(m_result) )
655 if ( ! may_update<boundary, interior, '1', transpose_result>(m_result)
656 && ! may_update<exterior, interior, '2', transpose_result>(m_result) )
662 inline void no_turns(segment_identifier const& seg_id)
664 // if those flags are set nothing will change
670 typename detail::sub_range_return_type<Geometry const>::type
671 range_ref = detail::sub_range(geometry, seg_id);
673 if ( boost::empty(range_ref) )
675 // TODO: throw an exception?
679 // TODO: possible optimization
680 // if the range is an interior ring we may use other IPs generated for this single geometry
681 // to know which other single geometries should be checked
683 // TODO: optimize! e.g. use spatial index
684 // O(N) - running it in a loop gives O(NM)
685 using detail::within::point_in_geometry;
686 int const pig = point_in_geometry(range::front(range_ref),
688 m_point_in_areal_strategy);
690 //BOOST_GEOMETRY_ASSERT(pig != 0);
693 update<interior, interior, '2', transpose_result>(m_result);
696 update<boundary, interior, '1', transpose_result>(m_result);
697 update<exterior, interior, '2', transpose_result>(m_result);
702 update<boundary, exterior, '1', transpose_result>(m_result);
703 update<interior, exterior, '2', transpose_result>(m_result);
707 // TODO: break if all things are set
708 // also some of them could be checked outside, before the analysis
709 // In this case we shouldn't relay just on dummy flags
710 // Flags should be initialized with proper values
711 // or the result should be checked directly
712 // THIS IS ALSO TRUE FOR OTHER ANALYSERS! in L/L and L/A
714 interrupt = m_flags == 7 || m_result.interrupt;
717 template <typename TurnIt>
718 inline void turns(TurnIt first, TurnIt last)
720 // if those flags are set nothing will change
721 if ( (m_flags & 6) == 6 )
726 bool found_ii = false;
727 bool found_uu = false;
729 for ( TurnIt it = first ; it != last ; ++it )
731 if ( it->operations[0].operation == overlay::operation_intersection
732 && it->operations[1].operation == overlay::operation_intersection )
736 else if ( it->operations[0].operation == overlay::operation_union
737 && it->operations[1].operation == overlay::operation_union )
743 return; // don't interrupt
747 // only i/i was generated for this ring
750 update<interior, interior, '2', transpose_result>(m_result);
753 //update<boundary, boundary, '0', transpose_result>(m_result);
755 update<boundary, interior, '1', transpose_result>(m_result);
756 update<exterior, interior, '2', transpose_result>(m_result);
760 // only u/u was generated for this ring
763 update<boundary, exterior, '1', transpose_result>(m_result);
764 update<interior, exterior, '2', transpose_result>(m_result);
768 interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future
771 Geometry const& geometry;
772 OtherGeometry const& other_geometry;
777 PointInArealStrategy const& m_point_in_areal_strategy;
781 template <std::size_t OpId>
782 class analyse_uncertain_rings
785 template <typename Analyser, typename TurnIt>
786 static inline void apply(Analyser & analyser, TurnIt first, TurnIt last)
791 for_preceding_rings(analyser, *first);
792 //analyser.per_turn(*first);
795 for ( ++first ; first != last ; ++first, ++prev )
798 if ( prev->operations[OpId].seg_id.multi_index
799 == first->operations[OpId].seg_id.multi_index )
802 if ( prev->operations[OpId].seg_id.ring_index
803 == first->operations[OpId].seg_id.ring_index )
805 //analyser.per_turn(*first);
807 // same multi, next ring
810 //analyser.end_ring(*prev);
811 analyser.turns(prev, first);
813 //if ( prev->operations[OpId].seg_id.ring_index + 1
814 // < first->operations[OpId].seg_id.ring_index)
816 for_no_turns_rings(analyser,
818 prev->operations[OpId].seg_id.ring_index + 1,
819 first->operations[OpId].seg_id.ring_index);
822 //analyser.per_turn(*first);
828 //analyser.end_ring(*prev);
829 analyser.turns(prev, first);
830 for_following_rings(analyser, *prev);
831 for_preceding_rings(analyser, *first);
832 //analyser.per_turn(*first);
835 if ( analyser.interrupt )
841 //analyser.end_ring(*prev);
842 analyser.turns(prev, first); // first == last
843 for_following_rings(analyser, *prev);
847 template <typename Analyser, typename Turn>
848 static inline void for_preceding_rings(Analyser & analyser, Turn const& turn)
850 segment_identifier const& seg_id = turn.operations[OpId].seg_id;
852 for_no_turns_rings(analyser, turn, -1, seg_id.ring_index);
855 template <typename Analyser, typename Turn>
856 static inline void for_following_rings(Analyser & analyser, Turn const& turn)
858 segment_identifier const& seg_id = turn.operations[OpId].seg_id;
861 count = boost::numeric_cast<signed_size_type>(
862 geometry::num_interior_rings(
863 detail::single_geometry(analyser.geometry, seg_id)));
865 for_no_turns_rings(analyser, turn, seg_id.ring_index + 1, count);
868 template <typename Analyser, typename Turn>
869 static inline void for_no_turns_rings(Analyser & analyser,
871 signed_size_type first,
872 signed_size_type last)
874 segment_identifier seg_id = turn.operations[OpId].seg_id;
876 for ( seg_id.ring_index = first ; seg_id.ring_index < last ; ++seg_id.ring_index )
878 analyser.no_turns(seg_id);
884 }} // namespace detail::relate
885 #endif // DOXYGEN_NO_DETAIL
887 }} // namespace boost::geometry
889 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP