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.
6 // Modifications copyright (c) 2013-2017 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<Geometry1, Geometry2>::turn_info turn_type;
227 std::vector<turn_type> turns;
229 interrupt_policy_areal_areal<Result> interrupt_policy(geometry1, geometry2, result);
231 turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy, intersection_strategy);
232 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
235 typedef typename IntersectionStrategy::template point_in_geometry_strategy
238 >::type point_in_areal_strategy12_type;
239 point_in_areal_strategy12_type point_in_areal_strategy12
240 = intersection_strategy.template get_point_in_geometry_strategy<Geometry1, Geometry2>();
241 typedef typename IntersectionStrategy::template point_in_geometry_strategy
244 >::type point_in_areal_strategy21_type;
245 point_in_areal_strategy21_type point_in_areal_strategy21
246 = intersection_strategy.template get_point_in_geometry_strategy<Geometry2, Geometry1>();
248 no_turns_aa_pred<Geometry2, Result, point_in_areal_strategy12_type, false>
249 pred1(geometry2, result, point_in_areal_strategy12);
250 for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
251 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
254 no_turns_aa_pred<Geometry1, Result, point_in_areal_strategy21_type, true>
255 pred2(geometry1, result, point_in_areal_strategy21);
256 for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
257 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
263 if ( may_update<interior, interior, '2'>(result)
264 || may_update<interior, exterior, '2'>(result)
265 || may_update<boundary, interior, '1'>(result)
266 || may_update<boundary, exterior, '1'>(result)
267 || may_update<exterior, interior, '2'>(result) )
270 typedef turns::less<0, turns::less_op_areal_areal<0> > less;
271 std::sort(turns.begin(), turns.end(), less());
273 /*if ( may_update<interior, exterior, '2'>(result)
274 || may_update<boundary, exterior, '1'>(result)
275 || may_update<boundary, interior, '1'>(result)
276 || may_update<exterior, interior, '2'>(result) )*/
278 // analyse sorted turns
279 turns_analyser<turn_type, 0> analyser;
280 analyse_each_turn(result, analyser, turns.begin(), turns.end());
282 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
286 if ( may_update<interior, interior, '2'>(result)
287 || may_update<interior, exterior, '2'>(result)
288 || may_update<boundary, interior, '1'>(result)
289 || may_update<boundary, exterior, '1'>(result)
290 || may_update<exterior, interior, '2'>(result) )
292 // analyse rings for which turns were not generated
293 // or only i/i or u/u was generated
294 uncertain_rings_analyser<0, Result, Geometry1, Geometry2, point_in_areal_strategy12_type>
295 rings_analyser(result, geometry1, geometry2, point_in_areal_strategy12);
296 analyse_uncertain_rings<0>::apply(rings_analyser, turns.begin(), turns.end());
298 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
303 if ( may_update<interior, interior, '2', true>(result)
304 || may_update<interior, exterior, '2', true>(result)
305 || may_update<boundary, interior, '1', true>(result)
306 || may_update<boundary, exterior, '1', true>(result)
307 || may_update<exterior, interior, '2', true>(result) )
310 typedef turns::less<1, turns::less_op_areal_areal<1> > less;
311 std::sort(turns.begin(), turns.end(), less());
313 /*if ( may_update<interior, exterior, '2', true>(result)
314 || may_update<boundary, exterior, '1', true>(result)
315 || may_update<boundary, interior, '1', true>(result)
316 || may_update<exterior, interior, '2', true>(result) )*/
318 // analyse sorted turns
319 turns_analyser<turn_type, 1> analyser;
320 analyse_each_turn(result, analyser, turns.begin(), turns.end());
322 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
326 if ( may_update<interior, interior, '2', true>(result)
327 || may_update<interior, exterior, '2', true>(result)
328 || may_update<boundary, interior, '1', true>(result)
329 || may_update<boundary, exterior, '1', true>(result)
330 || may_update<exterior, interior, '2', true>(result) )
332 // analyse rings for which turns were not generated
333 // or only i/i or u/u was generated
334 uncertain_rings_analyser<1, Result, Geometry2, Geometry1, point_in_areal_strategy21_type>
335 rings_analyser(result, geometry2, geometry1, point_in_areal_strategy21);
336 analyse_uncertain_rings<1>::apply(rings_analyser, turns.begin(), turns.end());
338 //if ( result.interrupt )
344 // interrupt policy which may be passed to get_turns to interrupt the analysis
345 // based on the info in the passed result/mask
346 template <typename Result>
347 class interrupt_policy_areal_areal
350 static bool const enabled = true;
352 interrupt_policy_areal_areal(Geometry1 const& geometry1,
353 Geometry2 const& geometry2,
356 , m_geometry1(geometry1)
357 , m_geometry2(geometry2)
360 template <typename Range>
361 inline bool apply(Range const& turns)
363 typedef typename boost::range_iterator<Range const>::type iterator;
365 for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
371 return m_result.interrupt;
375 template <std::size_t OpId, typename Turn>
376 inline void per_turn(Turn const& turn)
378 //static const std::size_t other_op_id = (OpId + 1) % 2;
379 static const bool transpose_result = OpId != 0;
381 overlay::operation_type const op = turn.operations[OpId].operation;
383 if ( op == overlay::operation_union )
386 /*if ( turn.operations[other_op_id].operation != overlay::operation_union )
388 update<interior, exterior, '2', transpose_result>(m_result);
389 update<boundary, exterior, '1', transpose_result>(m_result);
392 update<boundary, boundary, '0', transpose_result>(m_result);
394 else if ( op == overlay::operation_intersection )
397 /*if ( turn.operations[other_op_id].operation != overlay::operation_intersection )
399 // not correct e.g. for G1 touching G2 in a point where a hole is touching the exterior ring
400 // in this case 2 turns i/... and u/u will be generated for this IP
401 //update<interior, interior, '2', transpose_result>(m_result);
403 //update<boundary, interior, '1', transpose_result>(m_result);
406 update<boundary, boundary, '0', transpose_result>(m_result);
408 else if ( op == overlay::operation_continue )
410 update<boundary, boundary, '1', transpose_result>(m_result);
411 update<interior, interior, '2', transpose_result>(m_result);
413 else if ( op == overlay::operation_blocked )
415 update<boundary, boundary, '1', transpose_result>(m_result);
416 update<interior, exterior, '2', transpose_result>(m_result);
421 Geometry1 const& m_geometry1;
422 Geometry2 const& m_geometry2;
425 // This analyser should be used like Input or SinglePass Iterator
426 // IMPORTANT! It should be called also for the end iterator - last
427 template <typename TurnInfo, std::size_t OpId>
430 typedef typename TurnInfo::point_type turn_point_type;
432 static const std::size_t op_id = OpId;
433 static const std::size_t other_op_id = (OpId + 1) % 2;
434 static const bool transpose_result = OpId != 0;
438 : m_previous_turn_ptr(0)
439 , m_previous_operation(overlay::operation_none)
440 , m_enter_detected(false)
441 , m_exit_detected(false)
444 template <typename Result,
446 void apply(Result & result, TurnIt it)
448 //BOOST_GEOMETRY_ASSERT( it != last );
450 overlay::operation_type const op = it->operations[op_id].operation;
452 if ( op != overlay::operation_union
453 && op != overlay::operation_intersection
454 && op != overlay::operation_blocked
455 && op != overlay::operation_continue )
460 segment_identifier const& seg_id = it->operations[op_id].seg_id;
461 //segment_identifier const& other_id = it->operations[other_op_id].seg_id;
463 const bool first_in_range = m_seg_watcher.update(seg_id);
465 if ( m_previous_turn_ptr )
467 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
469 // real exit point - may be multiple
471 || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
474 m_exit_detected = false;
476 // fake exit point, reset state
477 else if ( op != overlay::operation_union )
479 m_exit_detected = false;
483 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
487 || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
489 update_enter(result);
490 m_enter_detected = false;
492 // fake entry point, reset state
493 else if ( op != overlay::operation_intersection )
495 m_enter_detected = false;
500 if ( op == overlay::operation_union )
502 // already set in interrupt policy
503 //update<boundary, boundary, '0', transpose_result>(m_result);
506 //if ( it->operations[other_op_id].operation != overlay::operation_union )
508 m_exit_detected = true;
511 else if ( op == overlay::operation_intersection )
514 if ( it->operations[other_op_id].operation != overlay::operation_intersection )
516 // this was set in the interrupt policy but it was wrong
517 // also here it's wrong since it may be a fake entry point
518 //update<interior, interior, '2', transpose_result>(result);
520 // already set in interrupt policy
521 //update<boundary, boundary, '0', transpose_result>(result);
522 m_enter_detected = true;
525 else if ( op == overlay::operation_blocked )
527 // already set in interrupt policy
529 else // if ( op == overlay::operation_continue )
531 // already set in interrupt policy
534 // store ref to previously analysed (valid) turn
535 m_previous_turn_ptr = boost::addressof(*it);
536 // and previously analysed (valid) operation
537 m_previous_operation = op;
541 template <typename Result>
542 void apply(Result & result)
544 //BOOST_GEOMETRY_ASSERT( first != last );
546 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
549 m_exit_detected = false;
552 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
554 update_enter(result);
555 m_enter_detected = false;
559 template <typename Result>
560 static inline void update_exit(Result & result)
562 update<interior, exterior, '2', transpose_result>(result);
563 update<boundary, exterior, '1', transpose_result>(result);
566 template <typename Result>
567 static inline void update_enter(Result & result)
569 update<interior, interior, '2', transpose_result>(result);
570 update<boundary, interior, '1', transpose_result>(result);
571 update<exterior, interior, '2', transpose_result>(result);
575 segment_watcher<same_ring> m_seg_watcher;
576 TurnInfo * m_previous_turn_ptr;
577 overlay::operation_type m_previous_operation;
578 bool m_enter_detected;
579 bool m_exit_detected;
582 // call analyser.apply() for each turn in range
583 // IMPORTANT! The analyser is also called for the end iterator - last
584 template <typename Result,
587 static inline void analyse_each_turn(Result & res,
589 TurnIt first, TurnIt last)
594 for ( TurnIt it = first ; it != last ; ++it )
596 analyser.apply(res, it);
598 if ( BOOST_GEOMETRY_CONDITION(res.interrupt) )
610 typename OtherGeometry,
611 typename PointInArealStrategy
613 class uncertain_rings_analyser
615 static const bool transpose_result = OpId != 0;
616 static const int other_id = (OpId + 1) % 2;
619 inline uncertain_rings_analyser(Result & result,
620 Geometry const& geom,
621 OtherGeometry const& other_geom,
622 PointInArealStrategy const& point_in_areal_strategy)
624 , other_geometry(other_geom)
625 , interrupt(result.interrupt) // just in case, could be false as well
627 , m_point_in_areal_strategy(point_in_areal_strategy)
630 // check which relations must be analysed
631 // NOTE: 1 and 4 could probably be connected
633 if ( ! may_update<interior, interior, '2', transpose_result>(m_result) )
638 if ( ! may_update<interior, exterior, '2', transpose_result>(m_result)
639 && ! may_update<boundary, exterior, '1', transpose_result>(m_result) )
644 if ( ! may_update<boundary, interior, '1', transpose_result>(m_result)
645 && ! may_update<exterior, interior, '2', transpose_result>(m_result) )
651 inline void no_turns(segment_identifier const& seg_id)
653 // if those flags are set nothing will change
659 typename detail::sub_range_return_type<Geometry const>::type
660 range_ref = detail::sub_range(geometry, seg_id);
662 if ( boost::empty(range_ref) )
664 // TODO: throw an exception?
668 // TODO: possible optimization
669 // if the range is an interior ring we may use other IPs generated for this single geometry
670 // to know which other single geometries should be checked
672 // TODO: optimize! e.g. use spatial index
673 // O(N) - running it in a loop gives O(NM)
674 using detail::within::point_in_geometry;
675 int const pig = point_in_geometry(range::front(range_ref),
677 m_point_in_areal_strategy);
679 //BOOST_GEOMETRY_ASSERT(pig != 0);
682 update<interior, interior, '2', transpose_result>(m_result);
685 update<boundary, interior, '1', transpose_result>(m_result);
686 update<exterior, interior, '2', transpose_result>(m_result);
691 update<boundary, exterior, '1', transpose_result>(m_result);
692 update<interior, exterior, '2', transpose_result>(m_result);
696 // TODO: break if all things are set
697 // also some of them could be checked outside, before the analysis
698 // In this case we shouldn't relay just on dummy flags
699 // Flags should be initialized with proper values
700 // or the result should be checked directly
701 // THIS IS ALSO TRUE FOR OTHER ANALYSERS! in L/L and L/A
703 interrupt = m_flags == 7 || m_result.interrupt;
706 template <typename TurnIt>
707 inline void turns(TurnIt first, TurnIt last)
709 // if those flags are set nothing will change
710 if ( (m_flags & 6) == 6 )
715 bool found_ii = false;
716 bool found_uu = false;
718 for ( TurnIt it = first ; it != last ; ++it )
720 if ( it->operations[0].operation == overlay::operation_intersection
721 && it->operations[1].operation == overlay::operation_intersection )
725 else if ( it->operations[0].operation == overlay::operation_union
726 && it->operations[1].operation == overlay::operation_union )
732 return; // don't interrupt
736 // only i/i was generated for this ring
739 update<interior, interior, '2', transpose_result>(m_result);
742 //update<boundary, boundary, '0', transpose_result>(m_result);
744 update<boundary, interior, '1', transpose_result>(m_result);
745 update<exterior, interior, '2', transpose_result>(m_result);
749 // only u/u was generated for this ring
752 update<boundary, exterior, '1', transpose_result>(m_result);
753 update<interior, exterior, '2', transpose_result>(m_result);
757 interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future
760 Geometry const& geometry;
761 OtherGeometry const& other_geometry;
766 PointInArealStrategy const& m_point_in_areal_strategy;
770 template <std::size_t OpId>
771 class analyse_uncertain_rings
774 template <typename Analyser, typename TurnIt>
775 static inline void apply(Analyser & analyser, TurnIt first, TurnIt last)
780 for_preceding_rings(analyser, *first);
781 //analyser.per_turn(*first);
784 for ( ++first ; first != last ; ++first, ++prev )
787 if ( prev->operations[OpId].seg_id.multi_index
788 == first->operations[OpId].seg_id.multi_index )
791 if ( prev->operations[OpId].seg_id.ring_index
792 == first->operations[OpId].seg_id.ring_index )
794 //analyser.per_turn(*first);
796 // same multi, next ring
799 //analyser.end_ring(*prev);
800 analyser.turns(prev, first);
802 //if ( prev->operations[OpId].seg_id.ring_index + 1
803 // < first->operations[OpId].seg_id.ring_index)
805 for_no_turns_rings(analyser,
807 prev->operations[OpId].seg_id.ring_index + 1,
808 first->operations[OpId].seg_id.ring_index);
811 //analyser.per_turn(*first);
817 //analyser.end_ring(*prev);
818 analyser.turns(prev, first);
819 for_following_rings(analyser, *prev);
820 for_preceding_rings(analyser, *first);
821 //analyser.per_turn(*first);
824 if ( analyser.interrupt )
830 //analyser.end_ring(*prev);
831 analyser.turns(prev, first); // first == last
832 for_following_rings(analyser, *prev);
836 template <typename Analyser, typename Turn>
837 static inline void for_preceding_rings(Analyser & analyser, Turn const& turn)
839 segment_identifier const& seg_id = turn.operations[OpId].seg_id;
841 for_no_turns_rings(analyser, turn, -1, seg_id.ring_index);
844 template <typename Analyser, typename Turn>
845 static inline void for_following_rings(Analyser & analyser, Turn const& turn)
847 segment_identifier const& seg_id = turn.operations[OpId].seg_id;
850 count = boost::numeric_cast<signed_size_type>(
851 geometry::num_interior_rings(
852 detail::single_geometry(analyser.geometry, seg_id)));
854 for_no_turns_rings(analyser, turn, seg_id.ring_index + 1, count);
857 template <typename Analyser, typename Turn>
858 static inline void for_no_turns_rings(Analyser & analyser,
860 signed_size_type first,
861 signed_size_type last)
863 segment_identifier seg_id = turn.operations[OpId].seg_id;
865 for ( seg_id.ring_index = first ; seg_id.ring_index < last ; ++seg_id.ring_index )
867 analyser.no_turns(seg_id);
873 }} // namespace detail::relate
874 #endif // DOXYGEN_NO_DETAIL
876 }} // namespace boost::geometry
878 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP