1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
4 
5 // This file was modified by Oracle on 2013, 2014, 2015, 2017.
6 // Modifications copyright (c) 2013-2017 Oracle and/or its affiliates.
7 
8 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
9 
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)
13 
14 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP
15 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP
16 
17 #include <boost/geometry/core/topological_dimension.hpp>
18 
19 #include <boost/geometry/util/condition.hpp>
20 #include <boost/geometry/util/range.hpp>
21 
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>
26 
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>
31 
32 namespace boost { namespace geometry
33 {
34 
35 #ifndef DOXYGEN_NO_DETAIL
36 namespace detail { namespace relate {
37 
38 // WARNING!
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!
41 
42 // may be used to set EI and EB for an Areal geometry for which no turns were generated
43 template
44 <
45     typename OtherAreal,
46     typename Result,
47     typename PointInArealStrategy,
48     bool TransposeResult
49 >
50 class no_turns_aa_pred
51 {
52 public:
no_turns_aa_pred(OtherAreal const & other_areal,Result & res,PointInArealStrategy const & point_in_areal_strategy)53     no_turns_aa_pred(OtherAreal const& other_areal,
54                      Result & res,
55                      PointInArealStrategy const& point_in_areal_strategy)
56         : m_result(res)
57         , m_point_in_areal_strategy(point_in_areal_strategy)
58         , m_other_areal(other_areal)
59         , m_flags(0)
60     {
61         // check which relations must be analysed
62 
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) )
66         {
67             m_flags |= 1;
68         }
69 
70         if ( ! may_update<interior, exterior, '2', TransposeResult>(m_result)
71           && ! may_update<boundary, exterior, '1', TransposeResult>(m_result) )
72         {
73             m_flags |= 2;
74         }
75     }
76 
77     template <typename Areal>
operator ()(Areal const & areal)78     bool operator()(Areal const& areal)
79     {
80         using detail::within::point_in_geometry;
81 
82         // if those flags are set nothing will change
83         if ( m_flags == 3 )
84         {
85             return false;
86         }
87 
88         typedef typename geometry::point_type<Areal>::type point_type;
89         point_type pt;
90         bool const ok = boost::geometry::point_on_border(pt, areal);
91 
92         // TODO: for now ignore, later throw an exception?
93         if ( !ok )
94         {
95             return true;
96         }
97 
98         // check if the areal is inside the other_areal
99         // TODO: This is O(N)
100         // Run in a loop O(NM) - optimize!
101         int const pig = point_in_geometry(pt,
102                                           m_other_areal,
103                                           m_point_in_areal_strategy);
104         //BOOST_GEOMETRY_ASSERT( pig != 0 );
105 
106         // inside
107         if ( pig > 0 )
108         {
109             update<interior, interior, '2', TransposeResult>(m_result);
110             update<boundary, interior, '1', TransposeResult>(m_result);
111             update<exterior, interior, '2', TransposeResult>(m_result);
112             m_flags |= 1;
113 
114             // TODO: OPTIMIZE!
115             // Only the interior rings of other ONE single geometry must be checked
116             // NOT all geometries
117 
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 )
123             {
124                 typename detail::sub_range_return_type<Areal const>::type
125                     range_ref = detail::sub_range(areal, ring_id);
126 
127                 if ( boost::empty(range_ref) )
128                 {
129                     // TODO: throw exception?
130                     continue; // ignore
131                 }
132 
133                 // TODO: O(N)
134                 // Optimize!
135                 int const hpig = point_in_geometry(range::front(range_ref),
136                                                    m_other_areal,
137                                                    m_point_in_areal_strategy);
138 
139                 // hole outside
140                 if ( hpig < 0 )
141                 {
142                     update<interior, exterior, '2', TransposeResult>(m_result);
143                     update<boundary, exterior, '1', TransposeResult>(m_result);
144                     m_flags |= 2;
145                     break;
146                 }
147             }
148         }
149         // outside
150         else
151         {
152             update<interior, exterior, '2', TransposeResult>(m_result);
153             update<boundary, exterior, '1', TransposeResult>(m_result);
154             m_flags |= 2;
155 
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 )
161             {
162                 typename detail::sub_range_return_type<Areal const>::type
163                     range_ref = detail::sub_range(areal, ring_id);
164 
165                 if ( boost::empty(range_ref) )
166                 {
167                     // TODO: throw exception?
168                     continue; // ignore
169                 }
170 
171                 // TODO: O(N)
172                 // Optimize!
173                 int const hpig = point_in_geometry(range::front(range_ref),
174                                                    m_other_areal,
175                                                    m_point_in_areal_strategy);
176 
177                 // hole inside
178                 if ( hpig > 0 )
179                 {
180                     update<interior, interior, '2', TransposeResult>(m_result);
181                     update<boundary, interior, '1', TransposeResult>(m_result);
182                     update<exterior, interior, '2', TransposeResult>(m_result);
183                     m_flags |= 1;
184                     break;
185                 }
186             }
187         }
188 
189         return m_flags != 3 && !m_result.interrupt;
190     }
191 
192 private:
193     Result & m_result;
194     PointInArealStrategy const& m_point_in_areal_strategy;
195     OtherAreal const& m_other_areal;
196     int m_flags;
197 };
198 
199 // The implementation of an algorithm calculating relate() for A/A
200 template <typename Geometry1, typename Geometry2>
201 struct areal_areal
202 {
203     // check Linear / Areal
204     BOOST_STATIC_ASSERT(topological_dimension<Geometry1>::value == 2
205                      && topological_dimension<Geometry2>::value == 2);
206 
207     static const bool interruption_enabled = true;
208 
209     typedef typename geometry::point_type<Geometry1>::type point1_type;
210     typedef typename geometry::point_type<Geometry2>::type point2_type;
211 
212     template <typename Result, typename IntersectionStrategy>
applyboost::geometry::detail::relate::areal_areal213     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2,
214                              Result & result,
215                              IntersectionStrategy const& intersection_strategy)
216     {
217 // TODO: If Areal geometry may have infinite size, change the following line:
218 
219         // The result should be FFFFFFFFF
220         relate::set<exterior, exterior, result_dimension<Geometry2>::value>(result);// FFFFFFFFd, d in [1,9] or T
221 
222         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
223             return;
224 
225         // get and analyse turns
226         typedef typename turns::get_turns<Geometry1, Geometry2>::turn_info turn_type;
227         std::vector<turn_type> turns;
228 
229         interrupt_policy_areal_areal<Result> interrupt_policy(geometry1, geometry2, result);
230 
231         turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy, intersection_strategy);
232         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
233             return;
234 
235         typedef typename IntersectionStrategy::template point_in_geometry_strategy
236             <
237                 Geometry1, Geometry2
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
242             <
243                 Geometry2, Geometry1
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>();
247 
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) )
252             return;
253 
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) )
258             return;
259 
260         if ( turns.empty() )
261             return;
262 
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) )
268         {
269             // sort turns
270             typedef turns::less<0, turns::less_op_areal_areal<0> > less;
271             std::sort(turns.begin(), turns.end(), less());
272 
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) )*/
277             {
278                 // analyse sorted turns
279                 turns_analyser<turn_type, 0> analyser;
280                 analyse_each_turn(result, analyser, turns.begin(), turns.end());
281 
282                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
283                     return;
284             }
285 
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) )
291             {
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());
297 
298                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
299                     return;
300             }
301         }
302 
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) )
308         {
309             // sort turns
310             typedef turns::less<1, turns::less_op_areal_areal<1> > less;
311             std::sort(turns.begin(), turns.end(), less());
312 
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) )*/
317             {
318                 // analyse sorted turns
319                 turns_analyser<turn_type, 1> analyser;
320                 analyse_each_turn(result, analyser, turns.begin(), turns.end());
321 
322                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
323                     return;
324             }
325 
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) )
331             {
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());
337 
338                 //if ( result.interrupt )
339                 //    return;
340             }
341         }
342     }
343 
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
348     {
349     public:
350         static bool const enabled = true;
351 
interrupt_policy_areal_areal(Geometry1 const & geometry1,Geometry2 const & geometry2,Result & result)352         interrupt_policy_areal_areal(Geometry1 const& geometry1,
353                                      Geometry2 const& geometry2,
354                                      Result & result)
355             : m_result(result)
356             , m_geometry1(geometry1)
357             , m_geometry2(geometry2)
358         {}
359 
360         template <typename Range>
apply(Range const & turns)361         inline bool apply(Range const& turns)
362         {
363             typedef typename boost::range_iterator<Range const>::type iterator;
364 
365             for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
366             {
367                 per_turn<0>(*it);
368                 per_turn<1>(*it);
369             }
370 
371             return m_result.interrupt;
372         }
373 
374     private:
375         template <std::size_t OpId, typename Turn>
per_turn(Turn const & turn)376         inline void per_turn(Turn const& turn)
377         {
378             //static const std::size_t other_op_id = (OpId + 1) % 2;
379             static const bool transpose_result = OpId != 0;
380 
381             overlay::operation_type const op = turn.operations[OpId].operation;
382 
383             if ( op == overlay::operation_union )
384             {
385                 // ignore u/u
386                 /*if ( turn.operations[other_op_id].operation != overlay::operation_union )
387                 {
388                     update<interior, exterior, '2', transpose_result>(m_result);
389                     update<boundary, exterior, '1', transpose_result>(m_result);
390                 }*/
391 
392                 update<boundary, boundary, '0', transpose_result>(m_result);
393             }
394             else if ( op == overlay::operation_intersection )
395             {
396                 // ignore i/i
397                 /*if ( turn.operations[other_op_id].operation != overlay::operation_intersection )
398                 {
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);
402 
403                     //update<boundary, interior, '1', transpose_result>(m_result);
404                 }*/
405 
406                 update<boundary, boundary, '0', transpose_result>(m_result);
407             }
408             else if ( op == overlay::operation_continue )
409             {
410                 update<boundary, boundary, '1', transpose_result>(m_result);
411                 update<interior, interior, '2', transpose_result>(m_result);
412             }
413             else if ( op == overlay::operation_blocked )
414             {
415                 update<boundary, boundary, '1', transpose_result>(m_result);
416                 update<interior, exterior, '2', transpose_result>(m_result);
417             }
418         }
419 
420         Result & m_result;
421         Geometry1 const& m_geometry1;
422         Geometry2 const& m_geometry2;
423     };
424 
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>
428     class turns_analyser
429     {
430         typedef typename TurnInfo::point_type turn_point_type;
431 
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;
435 
436     public:
turns_analyser()437         turns_analyser()
438             : m_previous_turn_ptr(0)
439             , m_previous_operation(overlay::operation_none)
440             , m_enter_detected(false)
441             , m_exit_detected(false)
442         {}
443 
444         template <typename Result,
445                   typename TurnIt>
apply(Result & result,TurnIt it)446         void apply(Result & result, TurnIt it)
447         {
448             //BOOST_GEOMETRY_ASSERT( it != last );
449 
450             overlay::operation_type const op = it->operations[op_id].operation;
451 
452             if ( op != overlay::operation_union
453               && op != overlay::operation_intersection
454               && op != overlay::operation_blocked
455               && op != overlay::operation_continue )
456             {
457                 return;
458             }
459 
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;
462 
463             const bool first_in_range = m_seg_watcher.update(seg_id);
464 
465             if ( m_previous_turn_ptr )
466             {
467                 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
468                 {
469                     // real exit point - may be multiple
470                     if ( first_in_range
471                       || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
472                     {
473                         update_exit(result);
474                         m_exit_detected = false;
475                     }
476                     // fake exit point, reset state
477                     else if ( op != overlay::operation_union )
478                     {
479                         m_exit_detected = false;
480                     }
481                 }
482                 /*else*/
483                 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
484                 {
485                     // real entry point
486                     if ( first_in_range
487                       || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
488                     {
489                         update_enter(result);
490                         m_enter_detected = false;
491                     }
492                     // fake entry point, reset state
493                     else if ( op != overlay::operation_intersection )
494                     {
495                         m_enter_detected = false;
496                     }
497                 }
498             }
499 
500             if ( op == overlay::operation_union )
501             {
502                 // already set in interrupt policy
503                 //update<boundary, boundary, '0', transpose_result>(m_result);
504 
505                 // ignore u/u
506                 //if ( it->operations[other_op_id].operation != overlay::operation_union )
507                 {
508                     m_exit_detected = true;
509                 }
510             }
511             else if ( op == overlay::operation_intersection )
512             {
513                 // ignore i/i
514                 if ( it->operations[other_op_id].operation != overlay::operation_intersection )
515                 {
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);
519 
520                     // already set in interrupt policy
521                     //update<boundary, boundary, '0', transpose_result>(result);
522                     m_enter_detected = true;
523                 }
524             }
525             else if ( op == overlay::operation_blocked )
526             {
527                 // already set in interrupt policy
528             }
529             else // if ( op == overlay::operation_continue )
530             {
531                 // already set in interrupt policy
532             }
533 
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;
538         }
539 
540         // it == last
541         template <typename Result>
apply(Result & result)542         void apply(Result & result)
543         {
544             //BOOST_GEOMETRY_ASSERT( first != last );
545 
546             if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
547             {
548                 update_exit(result);
549                 m_exit_detected = false;
550             }
551 
552             if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
553             {
554                 update_enter(result);
555                 m_enter_detected = false;
556             }
557         }
558 
559         template <typename Result>
update_exit(Result & result)560         static inline void update_exit(Result & result)
561         {
562             update<interior, exterior, '2', transpose_result>(result);
563             update<boundary, exterior, '1', transpose_result>(result);
564         }
565 
566         template <typename Result>
update_enter(Result & result)567         static inline void update_enter(Result & result)
568         {
569             update<interior, interior, '2', transpose_result>(result);
570             update<boundary, interior, '1', transpose_result>(result);
571             update<exterior, interior, '2', transpose_result>(result);
572         }
573 
574     private:
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;
580     };
581 
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,
585               typename Analyser,
586               typename TurnIt>
analyse_each_turnboost::geometry::detail::relate::areal_areal587     static inline void analyse_each_turn(Result & res,
588                                          Analyser & analyser,
589                                          TurnIt first, TurnIt last)
590     {
591         if ( first == last )
592             return;
593 
594         for ( TurnIt it = first ; it != last ; ++it )
595         {
596             analyser.apply(res, it);
597 
598             if ( BOOST_GEOMETRY_CONDITION(res.interrupt) )
599                 return;
600         }
601 
602         analyser.apply(res);
603     }
604 
605     template
606     <
607         std::size_t OpId,
608         typename Result,
609         typename Geometry,
610         typename OtherGeometry,
611         typename PointInArealStrategy
612     >
613     class uncertain_rings_analyser
614     {
615         static const bool transpose_result = OpId != 0;
616         static const int other_id = (OpId + 1) % 2;
617 
618     public:
uncertain_rings_analyser(Result & result,Geometry const & geom,OtherGeometry const & other_geom,PointInArealStrategy const & point_in_areal_strategy)619         inline uncertain_rings_analyser(Result & result,
620                                         Geometry const& geom,
621                                         OtherGeometry const& other_geom,
622                                         PointInArealStrategy const& point_in_areal_strategy)
623             : geometry(geom)
624             , other_geometry(other_geom)
625             , interrupt(result.interrupt) // just in case, could be false as well
626             , m_result(result)
627             , m_point_in_areal_strategy(point_in_areal_strategy)
628             , m_flags(0)
629         {
630             // check which relations must be analysed
631             // NOTE: 1 and 4 could probably be connected
632 
633             if ( ! may_update<interior, interior, '2', transpose_result>(m_result) )
634             {
635                 m_flags |= 1;
636             }
637 
638             if ( ! may_update<interior, exterior, '2', transpose_result>(m_result)
639               && ! may_update<boundary, exterior, '1', transpose_result>(m_result) )
640             {
641                 m_flags |= 2;
642             }
643 
644             if ( ! may_update<boundary, interior, '1', transpose_result>(m_result)
645               && ! may_update<exterior, interior, '2', transpose_result>(m_result) )
646             {
647                 m_flags |= 4;
648             }
649         }
650 
no_turns(segment_identifier const & seg_id)651         inline void no_turns(segment_identifier const& seg_id)
652         {
653             // if those flags are set nothing will change
654             if ( m_flags == 7 )
655             {
656                 return;
657             }
658 
659             typename detail::sub_range_return_type<Geometry const>::type
660                 range_ref = detail::sub_range(geometry, seg_id);
661 
662             if ( boost::empty(range_ref) )
663             {
664                 // TODO: throw an exception?
665                 return; // ignore
666             }
667 
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
671 
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),
676                                               other_geometry,
677                                               m_point_in_areal_strategy);
678 
679             //BOOST_GEOMETRY_ASSERT(pig != 0);
680             if ( pig > 0 )
681             {
682                 update<interior, interior, '2', transpose_result>(m_result);
683                 m_flags |= 1;
684 
685                 update<boundary, interior, '1', transpose_result>(m_result);
686                 update<exterior, interior, '2', transpose_result>(m_result);
687                 m_flags |= 4;
688             }
689             else
690             {
691                 update<boundary, exterior, '1', transpose_result>(m_result);
692                 update<interior, exterior, '2', transpose_result>(m_result);
693                 m_flags |= 2;
694             }
695 
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
702 
703             interrupt = m_flags == 7 || m_result.interrupt;
704         }
705 
706         template <typename TurnIt>
turns(TurnIt first,TurnIt last)707         inline void turns(TurnIt first, TurnIt last)
708         {
709             // if those flags are set nothing will change
710             if ( (m_flags & 6) == 6 )
711             {
712                 return;
713             }
714 
715             bool found_ii = false;
716             bool found_uu = false;
717 
718             for ( TurnIt it = first ; it != last ; ++it )
719             {
720                 if ( it->operations[0].operation == overlay::operation_intersection
721                   && it->operations[1].operation == overlay::operation_intersection )
722                 {
723                     found_ii = true;
724                 }
725                 else if ( it->operations[0].operation == overlay::operation_union
726                        && it->operations[1].operation == overlay::operation_union )
727                 {
728                     found_uu = true;
729                 }
730                 else // ignore
731                 {
732                     return; // don't interrupt
733                 }
734             }
735 
736             // only i/i was generated for this ring
737             if ( found_ii )
738             {
739                 update<interior, interior, '2', transpose_result>(m_result);
740                 m_flags |= 1;
741 
742                 //update<boundary, boundary, '0', transpose_result>(m_result);
743 
744                 update<boundary, interior, '1', transpose_result>(m_result);
745                 update<exterior, interior, '2', transpose_result>(m_result);
746                 m_flags |= 4;
747             }
748 
749             // only u/u was generated for this ring
750             if ( found_uu )
751             {
752                 update<boundary, exterior, '1', transpose_result>(m_result);
753                 update<interior, exterior, '2', transpose_result>(m_result);
754                 m_flags |= 2;
755             }
756 
757             interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future
758         }
759 
760         Geometry const& geometry;
761         OtherGeometry const& other_geometry;
762         bool interrupt;
763 
764     private:
765         Result & m_result;
766         PointInArealStrategy const& m_point_in_areal_strategy;
767         int m_flags;
768     };
769 
770     template <std::size_t OpId>
771     class analyse_uncertain_rings
772     {
773     public:
774         template <typename Analyser, typename TurnIt>
apply(Analyser & analyser,TurnIt first,TurnIt last)775         static inline void apply(Analyser & analyser, TurnIt first, TurnIt last)
776         {
777             if ( first == last )
778                 return;
779 
780             for_preceding_rings(analyser, *first);
781             //analyser.per_turn(*first);
782 
783             TurnIt prev = first;
784             for ( ++first ; first != last ; ++first, ++prev )
785             {
786                 // same multi
787                 if ( prev->operations[OpId].seg_id.multi_index
788                   == first->operations[OpId].seg_id.multi_index )
789                 {
790                     // same ring
791                     if ( prev->operations[OpId].seg_id.ring_index
792                       == first->operations[OpId].seg_id.ring_index )
793                     {
794                         //analyser.per_turn(*first);
795                     }
796                     // same multi, next ring
797                     else
798                     {
799                         //analyser.end_ring(*prev);
800                         analyser.turns(prev, first);
801 
802                         //if ( prev->operations[OpId].seg_id.ring_index + 1
803                         //   < first->operations[OpId].seg_id.ring_index)
804                         {
805                             for_no_turns_rings(analyser,
806                                                *first,
807                                                prev->operations[OpId].seg_id.ring_index + 1,
808                                                first->operations[OpId].seg_id.ring_index);
809                         }
810 
811                         //analyser.per_turn(*first);
812                     }
813                 }
814                 // next multi
815                 else
816                 {
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);
822                 }
823 
824                 if ( analyser.interrupt )
825                 {
826                     return;
827                 }
828             }
829 
830             //analyser.end_ring(*prev);
831             analyser.turns(prev, first); // first == last
832             for_following_rings(analyser, *prev);
833         }
834 
835     private:
836         template <typename Analyser, typename Turn>
for_preceding_rings(Analyser & analyser,Turn const & turn)837         static inline void for_preceding_rings(Analyser & analyser, Turn const& turn)
838         {
839             segment_identifier const& seg_id = turn.operations[OpId].seg_id;
840 
841             for_no_turns_rings(analyser, turn, -1, seg_id.ring_index);
842         }
843 
844         template <typename Analyser, typename Turn>
for_following_rings(Analyser & analyser,Turn const & turn)845         static inline void for_following_rings(Analyser & analyser, Turn const& turn)
846         {
847             segment_identifier const& seg_id = turn.operations[OpId].seg_id;
848 
849             signed_size_type
850                 count = boost::numeric_cast<signed_size_type>(
851                             geometry::num_interior_rings(
852                                 detail::single_geometry(analyser.geometry, seg_id)));
853 
854             for_no_turns_rings(analyser, turn, seg_id.ring_index + 1, count);
855         }
856 
857         template <typename Analyser, typename Turn>
for_no_turns_rings(Analyser & analyser,Turn const & turn,signed_size_type first,signed_size_type last)858         static inline void for_no_turns_rings(Analyser & analyser,
859                                               Turn const& turn,
860                                               signed_size_type first,
861                                               signed_size_type last)
862         {
863             segment_identifier seg_id = turn.operations[OpId].seg_id;
864 
865             for ( seg_id.ring_index = first ; seg_id.ring_index < last ; ++seg_id.ring_index )
866             {
867                 analyser.no_turns(seg_id);
868             }
869         }
870     };
871 };
872 
873 }} // namespace detail::relate
874 #endif // DOXYGEN_NO_DETAIL
875 
876 }} // namespace boost::geometry
877 
878 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP
879