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_LINEAR_AREAL_HPP
15 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP
16 
17 #include <boost/core/ignore_unused.hpp>
18 #include <boost/range/size.hpp>
19 
20 #include <boost/geometry/core/assert.hpp>
21 #include <boost/geometry/core/topological_dimension.hpp>
22 
23 #include <boost/geometry/util/condition.hpp>
24 #include <boost/geometry/util/range.hpp>
25 
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>
30 
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>
35 
36 #include <boost/geometry/views/detail/normalized_view.hpp>
37 
38 namespace boost { namespace geometry
39 {
40 
41 #ifndef DOXYGEN_NO_DETAIL
42 namespace detail { namespace relate {
43 
44 // WARNING!
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!
47 
48 // may be used to set IE and BE for a Linear geometry for which no turns were generated
49 template
50 <
51     typename Geometry2,
52     typename Result,
53     typename PointInArealStrategy,
54     typename BoundaryChecker,
55     bool TransposeResult
56 >
57 class no_turns_la_linestring_pred
58 {
59 public:
no_turns_la_linestring_pred(Geometry2 const & geometry2,Result & res,PointInArealStrategy const & point_in_areal_strategy,BoundaryChecker const & boundary_checker)60     no_turns_la_linestring_pred(Geometry2 const& geometry2,
61                                 Result & res,
62                                 PointInArealStrategy const& point_in_areal_strategy,
63                                 BoundaryChecker const& boundary_checker)
64         : m_geometry2(geometry2)
65         , m_result(res)
66         , m_point_in_areal_strategy(point_in_areal_strategy)
67         , m_boundary_checker(boundary_checker)
68         , m_interrupt_flags(0)
69     {
70         if ( ! may_update<interior, interior, '1', TransposeResult>(m_result) )
71         {
72             m_interrupt_flags |= 1;
73         }
74 
75         if ( ! may_update<interior, exterior, '1', TransposeResult>(m_result) )
76         {
77             m_interrupt_flags |= 2;
78         }
79 
80         if ( ! may_update<boundary, interior, '0', TransposeResult>(m_result) )
81         {
82             m_interrupt_flags |= 4;
83         }
84 
85         if ( ! may_update<boundary, exterior, '0', TransposeResult>(m_result) )
86         {
87             m_interrupt_flags |= 8;
88         }
89     }
90 
91     template <typename Linestring>
operator ()(Linestring const & linestring)92     bool operator()(Linestring const& linestring)
93     {
94         std::size_t const count = boost::size(linestring);
95 
96         // invalid input
97         if ( count < 2 )
98         {
99             // ignore
100             // TODO: throw an exception?
101             return true;
102         }
103 
104         // if those flags are set nothing will change
105         if ( m_interrupt_flags == 0xF )
106         {
107             return false;
108         }
109 
110         int const pig = detail::within::point_in_geometry(range::front(linestring),
111                                                           m_geometry2,
112                                                           m_point_in_areal_strategy);
113         //BOOST_GEOMETRY_ASSERT_MSG(pig != 0, "There should be no IPs");
114 
115         if ( pig > 0 )
116         {
117             update<interior, interior, '1', TransposeResult>(m_result);
118             m_interrupt_flags |= 1;
119         }
120         else
121         {
122             update<interior, exterior, '1', TransposeResult>(m_result);
123             m_interrupt_flags |= 2;
124         }
125 
126         // check if there is a boundary
127         if ( ( m_interrupt_flags & 0xC ) != 0xC // if wasn't already set
128           && ( m_boundary_checker.template
129                 is_endpoint_boundary<boundary_front>(range::front(linestring))
130             || m_boundary_checker.template
131                 is_endpoint_boundary<boundary_back>(range::back(linestring)) ) )
132         {
133             if ( pig > 0 )
134             {
135                 update<boundary, interior, '0', TransposeResult>(m_result);
136                 m_interrupt_flags |= 4;
137             }
138             else
139             {
140                 update<boundary, exterior, '0', TransposeResult>(m_result);
141                 m_interrupt_flags |= 8;
142             }
143         }
144 
145         return m_interrupt_flags != 0xF
146             && ! m_result.interrupt;
147     }
148 
149 private:
150     Geometry2 const& m_geometry2;
151     Result & m_result;
152     PointInArealStrategy const& m_point_in_areal_strategy;
153     BoundaryChecker const& m_boundary_checker;
154     unsigned m_interrupt_flags;
155 };
156 
157 // may be used to set EI and EB for an Areal geometry for which no turns were generated
158 template <typename Result, bool TransposeResult>
159 class no_turns_la_areal_pred
160 {
161 public:
no_turns_la_areal_pred(Result & res)162     no_turns_la_areal_pred(Result & res)
163         : m_result(res)
164         , interrupt(! may_update<interior, exterior, '2', TransposeResult>(m_result)
165                  && ! may_update<boundary, exterior, '1', TransposeResult>(m_result) )
166     {}
167 
168     template <typename Areal>
operator ()(Areal const & areal)169     bool operator()(Areal const& areal)
170     {
171         if ( interrupt )
172         {
173             return false;
174         }
175 
176         // TODO:
177         // handle empty/invalid geometries in a different way than below?
178 
179         typedef typename geometry::point_type<Areal>::type point_type;
180         point_type dummy;
181         bool const ok = boost::geometry::point_on_border(dummy, areal);
182 
183         // TODO: for now ignore, later throw an exception?
184         if ( !ok )
185         {
186             return true;
187         }
188 
189         update<interior, exterior, '2', TransposeResult>(m_result);
190         update<boundary, exterior, '1', TransposeResult>(m_result);
191 
192         return false;
193     }
194 
195 private:
196     Result & m_result;
197     bool const interrupt;
198 };
199 
200 // The implementation of an algorithm calculating relate() for L/A
201 template <typename Geometry1, typename Geometry2, bool TransposeResult = false>
202 struct linear_areal
203 {
204     // check Linear / Areal
205     BOOST_STATIC_ASSERT(topological_dimension<Geometry1>::value == 1
206                      && topological_dimension<Geometry2>::value == 2);
207 
208     static const bool interruption_enabled = true;
209 
210     typedef typename geometry::point_type<Geometry1>::type point1_type;
211     typedef typename geometry::point_type<Geometry2>::type point2_type;
212 
213     template <typename Geometry>
214         struct is_multi
215             : boost::is_base_of
216                 <
217                     multi_tag,
218                     typename tag<Geometry>::type
219                 >
220         {};
221 
222     template <typename Geom1, typename Geom2>
223     struct multi_turn_info
224         : turns::get_turns<Geom1, Geom2>::turn_info
225     {
multi_turn_infoboost::geometry::detail::relate::linear_areal::multi_turn_info226         multi_turn_info() : priority(0) {}
227         int priority; // single-geometry sorting priority
228     };
229 
230     template <typename Geom1, typename Geom2>
231     struct turn_info_type
232         : boost::mpl::if_c
233             <
234                 is_multi<Geometry2>::value,
235                 multi_turn_info<Geom1, Geom2>,
236                 typename turns::get_turns<Geom1, Geom2>::turn_info
237             >
238     {};
239 
240     template <typename Result, typename IntersectionStrategy>
applyboost::geometry::detail::relate::linear_areal241     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2,
242                              Result & result,
243                              IntersectionStrategy const& intersection_strategy)
244     {
245 // TODO: If Areal geometry may have infinite size, change the following line:
246 
247         // The result should be FFFFFFFFF
248         relate::set<exterior, exterior, result_dimension<Geometry2>::value, TransposeResult>(result);// FFFFFFFFd, d in [1,9] or T
249 
250         if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
251             return;
252 
253         // get and analyse turns
254         typedef typename turn_info_type<Geometry1, Geometry2>::type turn_type;
255         std::vector<turn_type> turns;
256 
257         interrupt_policy_linear_areal<Geometry2, Result> interrupt_policy(geometry2, result);
258 
259         turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy, intersection_strategy);
260         if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
261             return;
262 
263         boundary_checker<Geometry1> boundary_checker1(geometry1);
264         no_turns_la_linestring_pred
265             <
266                 Geometry2,
267                 Result,
268                 typename IntersectionStrategy::template point_in_geometry_strategy<Geometry1, Geometry2>::type,
269                 boundary_checker<Geometry1>,
270                 TransposeResult
271             > pred1(geometry2,
272                     result,
273                     intersection_strategy.template get_point_in_geometry_strategy<Geometry1, Geometry2>(),
274                     boundary_checker1);
275         for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
276         if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
277             return;
278 
279         no_turns_la_areal_pred<Result, !TransposeResult> pred2(result);
280         for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
281         if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
282             return;
283 
284         if ( turns.empty() )
285             return;
286 
287         // This is set here because in the case if empty Areal geometry were passed
288         // those shouldn't be set
289         relate::set<exterior, interior, '2', TransposeResult>(result);// FFFFFF2Fd
290         if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
291             return;
292 
293         {
294             sort_dispatch(turns.begin(), turns.end(), is_multi<Geometry2>());
295 
296             turns_analyser<turn_type> analyser;
297             analyse_each_turn(result, analyser,
298                               turns.begin(), turns.end(),
299                               geometry1, geometry2,
300                               boundary_checker1,
301                               intersection_strategy.get_side_strategy());
302 
303             if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) )
304                 return;
305         }
306 
307         // If 'c' (insersection_boundary) was not found we know that any Ls isn't equal to one of the Rings
308         if ( !interrupt_policy.is_boundary_found )
309         {
310             relate::set<exterior, boundary, '1', TransposeResult>(result);
311         }
312         // Don't calculate it if it's required
313         else if ( may_update<exterior, boundary, '1', TransposeResult>(result) )
314         {
315 // TODO: REVISE THIS CODE AND PROBABLY REWRITE SOME PARTS TO BE MORE HUMAN-READABLE
316 //       IN GENERAL IT ANALYSES THE RINGS OF AREAL GEOMETRY AND DETECTS THE ONES THAT
317 //       MAY OVERLAP THE INTERIOR OF LINEAR GEOMETRY (NO IPs OR NON-FAKE 'u' OPERATION)
318 // NOTE: For one case std::sort may be called again to sort data by operations for data already sorted by ring index
319 //       In the worst case scenario the complexity will be O( NlogN + R*(N/R)log(N/R) )
320 //       So always should remain O(NlogN) -> for R==1 <-> 1(N/1)log(N/1), for R==N <-> N(N/N)log(N/N)
321 //       Some benchmarking should probably be done to check if only one std::sort should be used
322 
323             // sort by multi_index and rind_index
324             std::sort(turns.begin(), turns.end(), less_ring());
325 
326             typedef typename std::vector<turn_type>::iterator turn_iterator;
327 
328             turn_iterator it = turns.begin();
329             segment_identifier * prev_seg_id_ptr = NULL;
330             // for each ring
331             for ( ; it != turns.end() ; )
332             {
333                 // it's the next single geometry
334                 if ( prev_seg_id_ptr == NULL
335                   || prev_seg_id_ptr->multi_index != it->operations[1].seg_id.multi_index )
336                 {
337                     // if the first ring has no IPs
338                     if ( it->operations[1].seg_id.ring_index > -1 )
339                     {
340                         // we can be sure that the exterior overlaps the boundary
341                         relate::set<exterior, boundary, '1', TransposeResult>(result);
342                         break;
343                     }
344                     // if there was some previous ring
345                     if ( prev_seg_id_ptr != NULL )
346                     {
347                         signed_size_type const next_ring_index = prev_seg_id_ptr->ring_index + 1;
348                         BOOST_GEOMETRY_ASSERT(next_ring_index >= 0);
349 
350                         // if one of the last rings of previous single geometry was ommited
351                         if ( static_cast<std::size_t>(next_ring_index)
352                                 < geometry::num_interior_rings(
353                                     single_geometry(geometry2, *prev_seg_id_ptr)) )
354                         {
355                             // we can be sure that the exterior overlaps the boundary
356                             relate::set<exterior, boundary, '1', TransposeResult>(result);
357                             break;
358                         }
359                     }
360                 }
361                 // if it's the same single geometry
362                 else /*if ( previous_multi_index == it->operations[1].seg_id.multi_index )*/
363                 {
364                     // and we jumped over one of the rings
365                     if ( prev_seg_id_ptr != NULL // just in case
366                       && prev_seg_id_ptr->ring_index + 1 < it->operations[1].seg_id.ring_index )
367                     {
368                         // we can be sure that the exterior overlaps the boundary
369                         relate::set<exterior, boundary, '1', TransposeResult>(result);
370                         break;
371                     }
372                 }
373 
374                 prev_seg_id_ptr = boost::addressof(it->operations[1].seg_id);
375 
376                 // find the next ring first iterator and check if the analysis should be performed
377                 has_boundary_intersection has_boundary_inters;
378                 turn_iterator next = find_next_ring(it, turns.end(), has_boundary_inters);
379 
380                 // if there is no 1d overlap with the boundary
381                 if ( !has_boundary_inters.result )
382                 {
383                     // we can be sure that the exterior overlaps the boundary
384                     relate::set<exterior, boundary, '1', TransposeResult>(result);
385                     break;
386                 }
387                 // else there is 1d overlap with the boundary so we must analyse the boundary
388                 else
389                 {
390                     // u, c
391                     typedef turns::less<1, turns::less_op_areal_linear<1> > less;
392                     std::sort(it, next, less());
393 
394                     // analyse
395                     areal_boundary_analyser<turn_type> analyser;
396                     for ( turn_iterator rit = it ; rit != next ; ++rit )
397                     {
398                         // if the analyser requests, break the search
399                         if ( !analyser.apply(it, rit, next) )
400                             break;
401                     }
402 
403                     // if the boundary of Areal goes out of the Linear
404                     if ( analyser.is_union_detected )
405                     {
406                         // we can be sure that the boundary of Areal overlaps the exterior of Linear
407                         relate::set<exterior, boundary, '1', TransposeResult>(result);
408                         break;
409                     }
410                 }
411 
412                 it = next;
413             }
414 
415             // if there was some previous ring
416             if ( prev_seg_id_ptr != NULL )
417             {
418                 signed_size_type const next_ring_index = prev_seg_id_ptr->ring_index + 1;
419                 BOOST_GEOMETRY_ASSERT(next_ring_index >= 0);
420 
421                 // if one of the last rings of previous single geometry was ommited
422                 if ( static_cast<std::size_t>(next_ring_index)
423                         < geometry::num_interior_rings(
424                             single_geometry(geometry2, *prev_seg_id_ptr)) )
425                 {
426                     // we can be sure that the exterior overlaps the boundary
427                     relate::set<exterior, boundary, '1', TransposeResult>(result);
428                 }
429             }
430         }
431     }
432 
433     template <typename It, typename Pred, typename Comp>
for_each_equal_rangeboost::geometry::detail::relate::linear_areal434     static void for_each_equal_range(It first, It last, Pred pred, Comp comp)
435     {
436         if ( first == last )
437             return;
438 
439         It first_equal = first;
440         It prev = first;
441         for ( ++first ; ; ++first, ++prev )
442         {
443             if ( first == last || !comp(*prev, *first) )
444             {
445                 pred(first_equal, first);
446                 first_equal = first;
447             }
448 
449             if ( first == last )
450                 break;
451         }
452     }
453 
454     struct same_ip
455     {
456         template <typename Turn>
operator ()boost::geometry::detail::relate::linear_areal::same_ip457         bool operator()(Turn const& left, Turn const& right) const
458         {
459             return left.operations[0].seg_id == right.operations[0].seg_id
460                 && left.operations[0].fraction == right.operations[0].fraction;
461         }
462     };
463 
464     struct same_ip_and_multi_index
465     {
466         template <typename Turn>
operator ()boost::geometry::detail::relate::linear_areal::same_ip_and_multi_index467         bool operator()(Turn const& left, Turn const& right) const
468         {
469             return same_ip()(left, right)
470                 && left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index;
471         }
472     };
473 
474     template <typename OpToPriority>
475     struct set_turns_group_priority
476     {
477         template <typename TurnIt>
operator ()boost::geometry::detail::relate::linear_areal::set_turns_group_priority478         void operator()(TurnIt first, TurnIt last) const
479         {
480             BOOST_GEOMETRY_ASSERT(first != last);
481             static OpToPriority op_to_priority;
482             // find the operation with the least priority
483             int least_priority = op_to_priority(first->operations[0]);
484             for ( TurnIt it = first + 1 ; it != last ; ++it )
485             {
486                 int priority = op_to_priority(it->operations[0]);
487                 if ( priority < least_priority )
488                     least_priority = priority;
489             }
490             // set the least priority for all turns of the group
491             for ( TurnIt it = first ; it != last ; ++it )
492             {
493                 it->priority = least_priority;
494             }
495         }
496     };
497 
498     template <typename SingleLess>
499     struct sort_turns_group
500     {
501         struct less
502         {
503             template <typename Turn>
operator ()boost::geometry::detail::relate::linear_areal::sort_turns_group::less504             bool operator()(Turn const& left, Turn const& right) const
505             {
506                 return left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index ?
507                     SingleLess()(left, right) :
508                     left.priority < right.priority;
509             }
510         };
511 
512         template <typename TurnIt>
operator ()boost::geometry::detail::relate::linear_areal::sort_turns_group513         void operator()(TurnIt first, TurnIt last) const
514         {
515             std::sort(first, last, less());
516         }
517     };
518 
519     template <typename TurnIt>
sort_dispatchboost::geometry::detail::relate::linear_areal520     static void sort_dispatch(TurnIt first, TurnIt last, boost::true_type const& /*is_multi*/)
521     {
522         // sort turns by Linear seg_id, then by fraction, then by other multi_index
523         typedef turns::less<0, turns::less_other_multi_index<0> > less;
524         std::sort(first, last, less());
525 
526         // For the same IP and multi_index - the same other's single geometry
527         // set priorities as the least operation found for the whole single geometry
528         // so e.g. single geometries containing 'u' will always be before those only containing 'i'
529         typedef turns::op_to_int<0,2,3,1,4,0> op_to_int_xuic;
530         for_each_equal_range(first, last,
531                              set_turns_group_priority<op_to_int_xuic>(), // least operation in xuic order
532                              same_ip_and_multi_index()); // other's multi_index
533 
534         // When priorities for single geometries are set now sort turns for the same IP
535         // if multi_index is the same sort them according to the single-less
536         // else use priority of the whole single-geometry set earlier
537         typedef turns::less<0, turns::less_op_linear_areal_single<0> > single_less;
538         for_each_equal_range(first, last,
539                              sort_turns_group<single_less>(),
540                              same_ip());
541     }
542 
543     template <typename TurnIt>
sort_dispatchboost::geometry::detail::relate::linear_areal544     static void sort_dispatch(TurnIt first, TurnIt last, boost::false_type const& /*is_multi*/)
545     {
546         // sort turns by Linear seg_id, then by fraction, then
547         // for same ring id: x, u, i, c
548         // for different ring id: c, i, u, x
549         typedef turns::less<0, turns::less_op_linear_areal_single<0> > less;
550         std::sort(first, last, less());
551     }
552 
553 
554     // interrupt policy which may be passed to get_turns to interrupt the analysis
555     // based on the info in the passed result/mask
556     template <typename Areal, typename Result>
557     class interrupt_policy_linear_areal
558     {
559     public:
560         static bool const enabled = true;
561 
interrupt_policy_linear_areal(Areal const & areal,Result & result)562         interrupt_policy_linear_areal(Areal const& areal, Result & result)
563             : m_result(result), m_areal(areal)
564             , is_boundary_found(false)
565         {}
566 
567 // TODO: since we update result for some operations here, we may not do it in the analyser!
568 
569         template <typename Range>
apply(Range const & turns)570         inline bool apply(Range const& turns)
571         {
572             typedef typename boost::range_iterator<Range const>::type iterator;
573 
574             for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
575             {
576                 if ( it->operations[0].operation == overlay::operation_intersection )
577                 {
578                     bool const no_interior_rings
579                         = geometry::num_interior_rings(
580                                 single_geometry(m_areal, it->operations[1].seg_id)) == 0;
581 
582                     // WARNING! THIS IS TRUE ONLY IF THE POLYGON IS SIMPLE!
583                     // OR WITHOUT INTERIOR RINGS (AND OF COURSE VALID)
584                     if ( no_interior_rings )
585                         update<interior, interior, '1', TransposeResult>(m_result);
586                 }
587                 else if ( it->operations[0].operation == overlay::operation_continue )
588                 {
589                     update<interior, boundary, '1', TransposeResult>(m_result);
590                     is_boundary_found = true;
591                 }
592                 else if ( ( it->operations[0].operation == overlay::operation_union
593                          || it->operations[0].operation == overlay::operation_blocked )
594                        && it->operations[0].position == overlay::position_middle )
595                 {
596 // TODO: here we could also check the boundaries and set BB at this point
597                     update<interior, boundary, '0', TransposeResult>(m_result);
598                 }
599             }
600 
601             return m_result.interrupt;
602         }
603 
604     private:
605         Result & m_result;
606         Areal const& m_areal;
607 
608     public:
609         bool is_boundary_found;
610     };
611 
612     // This analyser should be used like Input or SinglePass Iterator
613     // IMPORTANT! It should be called also for the end iterator - last
614     template <typename TurnInfo>
615     class turns_analyser
616     {
617         typedef typename TurnInfo::point_type turn_point_type;
618 
619         static const std::size_t op_id = 0;
620         static const std::size_t other_op_id = 1;
621 
622     public:
turns_analyser()623         turns_analyser()
624             : m_previous_turn_ptr(NULL)
625             , m_previous_operation(overlay::operation_none)
626             , m_boundary_counter(0)
627             , m_interior_detected(false)
628             , m_first_interior_other_id_ptr(NULL)
629             , m_first_from_unknown(false)
630             , m_first_from_unknown_boundary_detected(false)
631         {}
632 
633         template <typename Result,
634                   typename TurnIt,
635                   typename Geometry,
636                   typename OtherGeometry,
637                   typename BoundaryChecker,
638                   typename SideStrategy>
apply(Result & res,TurnIt it,Geometry const & geometry,OtherGeometry const & other_geometry,BoundaryChecker const & boundary_checker,SideStrategy const & side_strategy)639         void apply(Result & res, TurnIt it,
640                    Geometry const& geometry,
641                    OtherGeometry const& other_geometry,
642                    BoundaryChecker const& boundary_checker,
643                    SideStrategy const& side_strategy)
644         {
645             overlay::operation_type op = it->operations[op_id].operation;
646 
647             if ( op != overlay::operation_union
648               && op != overlay::operation_intersection
649               && op != overlay::operation_blocked
650               && op != overlay::operation_continue ) // operation_boundary / operation_boundary_intersection
651             {
652                 return;
653             }
654 
655             segment_identifier const& seg_id = it->operations[op_id].seg_id;
656             segment_identifier const& other_id = it->operations[other_op_id].seg_id;
657 
658             const bool first_in_range = m_seg_watcher.update(seg_id);
659 
660             // TODO: should apply() for the post-last ip be called if first_in_range ?
661             // this would unify how last points in ranges are handled
662             // possibly replacing parts of the code below
663             // e.g. for is_multi and m_interior_detected
664 
665             // handle possible exit
666             bool fake_enter_detected = false;
667             if ( m_exit_watcher.get_exit_operation() == overlay::operation_union )
668             {
669                 // real exit point - may be multiple
670                 // we know that we entered and now we exit
671                 if ( ! turn_on_the_same_ip<op_id>(m_exit_watcher.get_exit_turn(), *it) )
672                 {
673                     m_exit_watcher.reset_detected_exit();
674 
675                     update<interior, exterior, '1', TransposeResult>(res);
676 
677                     // next single geometry
678                     if ( first_in_range && m_previous_turn_ptr )
679                     {
680                         // NOTE: similar code is in the post-last-ip-apply()
681                         segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
682 
683                         bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
684                                                     range::back(sub_range(geometry, prev_seg_id)),
685                                                     boundary_checker);
686 
687                         // if there is a boundary on the last point
688                         if ( prev_back_b )
689                         {
690                             update<boundary, exterior, '0', TransposeResult>(res);
691                         }
692                     }
693                 }
694                 // fake exit point, reset state
695                 else if ( op == overlay::operation_intersection
696                         || op == overlay::operation_continue ) // operation_boundary
697                 {
698                     m_exit_watcher.reset_detected_exit();
699                     fake_enter_detected = true;
700                 }
701             }
702             else if ( m_exit_watcher.get_exit_operation() == overlay::operation_blocked )
703             {
704                 // ignore multiple BLOCKs for this same single geometry1
705                 if ( op == overlay::operation_blocked
706                   && seg_id.multi_index == m_previous_turn_ptr->operations[op_id].seg_id.multi_index )
707                 {
708                     return;
709                 }
710 
711                 if ( ( op == overlay::operation_intersection
712                     || op == overlay::operation_continue )
713                   && turn_on_the_same_ip<op_id>(m_exit_watcher.get_exit_turn(), *it) )
714                 {
715                     fake_enter_detected = true;
716                 }
717 
718                 m_exit_watcher.reset_detected_exit();
719             }
720 
721             if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
722               && m_first_from_unknown )
723             {
724                 // For MultiPolygon many x/u operations may be generated as a first IP
725                 // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString
726                 // then we know that the LineString is outside
727                 // Similar with the u/u turns, if it was the first one it doesn't mean that the
728                 // Linestring came from the exterior
729                 if ( ( m_previous_operation == overlay::operation_blocked
730                     && ( op != overlay::operation_blocked // operation different than block
731                         || seg_id.multi_index != m_previous_turn_ptr->operations[op_id].seg_id.multi_index ) ) // or the next single-geometry
732                   || ( m_previous_operation == overlay::operation_union
733                     && ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
734                    )
735                 {
736                     update<interior, exterior, '1', TransposeResult>(res);
737                     if ( m_first_from_unknown_boundary_detected )
738                     {
739                         update<boundary, exterior, '0', TransposeResult>(res);
740                     }
741 
742                     m_first_from_unknown = false;
743                     m_first_from_unknown_boundary_detected = false;
744                 }
745             }
746 
747 // NOTE: THE WHOLE m_interior_detected HANDLING IS HERE BECAUSE WE CAN'T EFFICIENTLY SORT TURNS (CORRECTLY)
748 // BECAUSE THE SAME IP MAY BE REPRESENTED BY TWO SEGMENTS WITH DIFFERENT DISTANCES
749 // IT WOULD REQUIRE THE CALCULATION OF MAX DISTANCE
750 // TODO: WE COULD GET RID OF THE TEST IF THE DISTANCES WERE NORMALIZED
751 
752 // UPDATE: THEY SHOULD BE NORMALIZED NOW
753 
754 // TODO: THIS IS POTENTIALLY ERROREOUS!
755 // THIS ALGORITHM DEPENDS ON SOME SPECIFIC SEQUENCE OF OPERATIONS
756 // IT WOULD GIVE WRONG RESULTS E.G.
757 // IN THE CASE OF SELF-TOUCHING POINT WHEN 'i' WOULD BE BEFORE 'u'
758 
759             // handle the interior overlap
760             if ( m_interior_detected )
761             {
762                 BOOST_GEOMETRY_ASSERT_MSG(m_previous_turn_ptr, "non-NULL ptr expected");
763 
764                 // real interior overlap
765                 if ( ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) )
766                 {
767                     update<interior, interior, '1', TransposeResult>(res);
768                     m_interior_detected = false;
769 
770                     // new range detected - reset previous state and check the boundary
771                     if ( first_in_range )
772                     {
773                         segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
774 
775                         bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
776                                                     range::back(sub_range(geometry, prev_seg_id)),
777                                                     boundary_checker);
778 
779                         // if there is a boundary on the last point
780                         if ( prev_back_b )
781                         {
782                             update<boundary, interior, '0', TransposeResult>(res);
783                         }
784 
785                         // The exit_watcher is reset below
786                         // m_exit_watcher.reset();
787                     }
788                 }
789                 // fake interior overlap
790                 else if ( op == overlay::operation_continue )
791                 {
792                     m_interior_detected = false;
793                 }
794                 else if ( op == overlay::operation_union )
795                 {
796 // TODO: this probably is not a good way of handling the interiors/enters
797 //       the solution similar to exit_watcher would be more robust
798 //       all enters should be kept and handled.
799 //       maybe integrate it with the exit_watcher -> enter_exit_watcher
800                     if ( m_first_interior_other_id_ptr
801                       && m_first_interior_other_id_ptr->multi_index == other_id.multi_index )
802                     {
803                         m_interior_detected = false;
804                     }
805                 }
806             }
807 
808             // NOTE: If post-last-ip apply() was called this wouldn't be needed
809             if ( first_in_range )
810             {
811                 m_exit_watcher.reset();
812                 m_boundary_counter = 0;
813                 m_first_from_unknown = false;
814                 m_first_from_unknown_boundary_detected = false;
815             }
816 
817             // i/u, c/u
818             if ( op == overlay::operation_intersection
819               || op == overlay::operation_continue ) // operation_boundary/operation_boundary_intersection
820             {
821                 bool const first_point = first_in_range || m_first_from_unknown;
822                 bool no_enters_detected = m_exit_watcher.is_outside();
823                 m_exit_watcher.enter(*it);
824 
825                 if ( op == overlay::operation_intersection )
826                 {
827                     if ( m_boundary_counter > 0 && it->operations[op_id].is_collinear )
828                         --m_boundary_counter;
829 
830                     if ( m_boundary_counter == 0 )
831                     {
832                         // interiors overlaps
833                         //update<interior, interior, '1', TransposeResult>(res);
834 
835 // TODO: think about the implementation of the more robust version
836 //       this way only the first enter will be handled
837                         if ( !m_interior_detected )
838                         {
839                             // don't update now
840                             // we might enter a boundary of some other ring on the same IP
841                             m_interior_detected = true;
842                             m_first_interior_other_id_ptr = boost::addressof(other_id);
843                         }
844                     }
845                 }
846                 else // operation_boundary
847                 {
848                     // don't add to the count for all met boundaries
849                     // only if this is the "new" boundary
850                     if ( first_point || !it->operations[op_id].is_collinear )
851                         ++m_boundary_counter;
852 
853                     update<interior, boundary, '1', TransposeResult>(res);
854                 }
855 
856                 bool const this_b
857                     = is_ip_on_boundary<boundary_front>(it->point,
858                                                         it->operations[op_id],
859                                                         boundary_checker,
860                                                         seg_id);
861                 // going inside on boundary point
862                 if ( this_b )
863                 {
864                     update<boundary, boundary, '0', TransposeResult>(res);
865                 }
866                 // going inside on non-boundary point
867                 else
868                 {
869                     update<interior, boundary, '0', TransposeResult>(res);
870 
871                     // if we didn't enter in the past, we were outside
872                     if ( no_enters_detected
873                       && ! fake_enter_detected
874                       && it->operations[op_id].position != overlay::position_front )
875                     {
876 // TODO: calculate_from_inside() is only needed if the current Linestring is not closed
877                         bool const from_inside = first_point
878                                               && calculate_from_inside(geometry,
879                                                                        other_geometry,
880                                                                        *it,
881                                                                        side_strategy);
882 
883                         if ( from_inside )
884                             update<interior, interior, '1', TransposeResult>(res);
885                         else
886                             update<interior, exterior, '1', TransposeResult>(res);
887 
888                         // if it's the first IP then the first point is outside
889                         if ( first_point )
890                         {
891                             bool const front_b = is_endpoint_on_boundary<boundary_front>(
892                                                     range::front(sub_range(geometry, seg_id)),
893                                                     boundary_checker);
894 
895                             // if there is a boundary on the first point
896                             if ( front_b )
897                             {
898                                 if ( from_inside )
899                                     update<boundary, interior, '0', TransposeResult>(res);
900                                 else
901                                     update<boundary, exterior, '0', TransposeResult>(res);
902                             }
903                         }
904                     }
905                 }
906 
907                 if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) )
908                 {
909                     m_first_from_unknown = false;
910                     m_first_from_unknown_boundary_detected = false;
911                 }
912             }
913             // u/u, x/u
914             else if ( op == overlay::operation_union || op == overlay::operation_blocked )
915             {
916                 bool const op_blocked = op == overlay::operation_blocked;
917                 bool const no_enters_detected = m_exit_watcher.is_outside()
918 // TODO: is this condition ok?
919 // TODO: move it into the exit_watcher?
920                     && m_exit_watcher.get_exit_operation() == overlay::operation_none;
921 
922                 if ( op == overlay::operation_union )
923                 {
924                     if ( m_boundary_counter > 0 && it->operations[op_id].is_collinear )
925                         --m_boundary_counter;
926                 }
927                 else // overlay::operation_blocked
928                 {
929                     m_boundary_counter = 0;
930                 }
931 
932                 // we're inside, possibly going out right now
933                 if ( ! no_enters_detected )
934                 {
935                     if ( op_blocked
936                       && it->operations[op_id].position == overlay::position_back ) // ignore spikes!
937                     {
938                         // check if this is indeed the boundary point
939                         // NOTE: is_ip_on_boundary<>() should be called here but the result will be the same
940                         if ( is_endpoint_on_boundary<boundary_back>(it->point, boundary_checker) )
941                         {
942                             update<boundary, boundary, '0', TransposeResult>(res);
943                         }
944                     }
945                     // union, inside, but no exit -> collinear on self-intersection point
946                     // not needed since we're already inside the boundary
947                     /*else if ( !exit_detected )
948                     {
949                         update<interior, boundary, '0', TransposeResult>(res);
950                     }*/
951                 }
952                 // we're outside or inside and this is the first turn
953                 else
954                 {
955                     bool const this_b = is_ip_on_boundary<boundary_any>(it->point,
956                                                                         it->operations[op_id],
957                                                                         boundary_checker,
958                                                                         seg_id);
959                     // if current IP is on boundary of the geometry
960                     if ( this_b )
961                     {
962                         update<boundary, boundary, '0', TransposeResult>(res);
963                     }
964                     // if current IP is not on boundary of the geometry
965                     else
966                     {
967                         update<interior, boundary, '0', TransposeResult>(res);
968                     }
969 
970                     // TODO: very similar code is used in the handling of intersection
971                     if ( it->operations[op_id].position != overlay::position_front )
972                     {
973 // TODO: calculate_from_inside() is only needed if the current Linestring is not closed
974                         // NOTE: this is not enough for MultiPolygon and operation_blocked
975                         // For LS/MultiPolygon multiple x/u turns may be generated
976                         // the first checked Polygon may be the one which LS is outside for.
977                         bool const first_point = first_in_range || m_first_from_unknown;
978                         bool const first_from_inside = first_point
979                                                     && calculate_from_inside(geometry,
980                                                                              other_geometry,
981                                                                              *it,
982                                                                              side_strategy);
983                         if ( first_from_inside )
984                         {
985                             update<interior, interior, '1', TransposeResult>(res);
986 
987                             // notify the exit_watcher that we started inside
988                             m_exit_watcher.enter(*it);
989                             // and reset unknown flags since we know that we started inside
990                             m_first_from_unknown = false;
991                             m_first_from_unknown_boundary_detected = false;
992                         }
993                         else
994                         {
995                             if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
996                               /*&& ( op == overlay::operation_blocked
997                                 || op == overlay::operation_union )*/ ) // if we're here it's u or x
998                             {
999                                 m_first_from_unknown = true;
1000                             }
1001                             else
1002                             {
1003                                 update<interior, exterior, '1', TransposeResult>(res);
1004                             }
1005                         }
1006 
1007                         // first IP on the last segment point - this means that the first point is outside or inside
1008                         if ( first_point && ( !this_b || op_blocked ) )
1009                         {
1010                             bool const front_b = is_endpoint_on_boundary<boundary_front>(
1011                                                     range::front(sub_range(geometry, seg_id)),
1012                                                     boundary_checker);
1013 
1014                             // if there is a boundary on the first point
1015                             if ( front_b )
1016                             {
1017                                 if ( first_from_inside )
1018                                 {
1019                                     update<boundary, interior, '0', TransposeResult>(res);
1020                                 }
1021                                 else
1022                                 {
1023                                     if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
1024                                       /*&& ( op == overlay::operation_blocked
1025                                         || op == overlay::operation_union )*/ ) // if we're here it's u or x
1026                                     {
1027                                         BOOST_GEOMETRY_ASSERT(m_first_from_unknown);
1028                                         m_first_from_unknown_boundary_detected = true;
1029                                     }
1030                                     else
1031                                     {
1032                                         update<boundary, exterior, '0', TransposeResult>(res);
1033                                     }
1034                                 }
1035                             }
1036                         }
1037                     }
1038                 }
1039 
1040                 // if we're going along a boundary, we exit only if the linestring was collinear
1041                 if ( m_boundary_counter == 0
1042                   || it->operations[op_id].is_collinear )
1043                 {
1044                     // notify the exit watcher about the possible exit
1045                     m_exit_watcher.exit(*it);
1046                 }
1047             }
1048 
1049             // store ref to previously analysed (valid) turn
1050             m_previous_turn_ptr = boost::addressof(*it);
1051             // and previously analysed (valid) operation
1052             m_previous_operation = op;
1053         }
1054 
1055         // it == last
1056         template <typename Result,
1057                   typename TurnIt,
1058                   typename Geometry,
1059                   typename OtherGeometry,
1060                   typename BoundaryChecker>
apply(Result & res,TurnIt first,TurnIt last,Geometry const & geometry,OtherGeometry const &,BoundaryChecker const & boundary_checker)1061         void apply(Result & res,
1062                    TurnIt first, TurnIt last,
1063                    Geometry const& geometry,
1064                    OtherGeometry const& /*other_geometry*/,
1065                    BoundaryChecker const& boundary_checker)
1066         {
1067             boost::ignore_unused(first, last);
1068             //BOOST_GEOMETRY_ASSERT( first != last );
1069 
1070             // For MultiPolygon many x/u operations may be generated as a first IP
1071             // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString
1072             // then we know that the LineString is outside
1073             if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value )
1074               && m_first_from_unknown )
1075             {
1076                 update<interior, exterior, '1', TransposeResult>(res);
1077                 if ( m_first_from_unknown_boundary_detected )
1078                 {
1079                     update<boundary, exterior, '0', TransposeResult>(res);
1080                 }
1081 
1082                 // done below
1083                 //m_first_from_unknown = false;
1084                 //m_first_from_unknown_boundary_detected = false;
1085             }
1086 
1087             // here, the possible exit is the real one
1088             // we know that we entered and now we exit
1089             if ( /*m_exit_watcher.get_exit_operation() == overlay::operation_union // THIS CHECK IS REDUNDANT
1090                 ||*/ m_previous_operation == overlay::operation_union
1091                 && !m_interior_detected )
1092             {
1093                 // for sure
1094                 update<interior, exterior, '1', TransposeResult>(res);
1095 
1096                 BOOST_GEOMETRY_ASSERT(first != last);
1097                 BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr);
1098 
1099                 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1100 
1101                 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1102                                             range::back(sub_range(geometry, prev_seg_id)),
1103                                             boundary_checker);
1104 
1105                 // if there is a boundary on the last point
1106                 if ( prev_back_b )
1107                 {
1108                     update<boundary, exterior, '0', TransposeResult>(res);
1109                 }
1110             }
1111             // we might enter some Areal and didn't go out,
1112             else if ( m_previous_operation == overlay::operation_intersection
1113                    || m_interior_detected )
1114             {
1115                 // just in case
1116                 update<interior, interior, '1', TransposeResult>(res);
1117                 m_interior_detected = false;
1118 
1119                 BOOST_GEOMETRY_ASSERT(first != last);
1120                 BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr);
1121 
1122                 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1123 
1124                 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1125                                             range::back(sub_range(geometry, prev_seg_id)),
1126                                             boundary_checker);
1127 
1128                 // if there is a boundary on the last point
1129                 if ( prev_back_b )
1130                 {
1131                     update<boundary, interior, '0', TransposeResult>(res);
1132                 }
1133             }
1134 
1135             // This condition may be false if the Linestring is lying on the Polygon's collinear spike
1136             // if Polygon's spikes are not handled in get_turns() or relate() (they currently aren't)
1137             //BOOST_GEOMETRY_ASSERT_MSG(m_previous_operation != overlay::operation_continue,
1138             //                    "Unexpected operation! Probably the error in get_turns(L,A) or relate(L,A)");
1139             // Currently one c/c turn is generated for the exit
1140             //   when a Linestring is going out on a collinear spike
1141             // When a Linestring is going in on a collinear spike
1142             //   the turn is not generated for the entry
1143             // So assume it's the former
1144             if ( m_previous_operation == overlay::operation_continue )
1145             {
1146                 update<interior, exterior, '1', TransposeResult>(res);
1147 
1148                 segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id;
1149 
1150                 bool const prev_back_b = is_endpoint_on_boundary<boundary_back>(
1151                                             range::back(sub_range(geometry, prev_seg_id)),
1152                                             boundary_checker);
1153 
1154                 // if there is a boundary on the last point
1155                 if ( prev_back_b )
1156                 {
1157                     update<boundary, exterior, '0', TransposeResult>(res);
1158                 }
1159             }
1160 
1161             // Reset exit watcher before the analysis of the next Linestring
1162             m_exit_watcher.reset();
1163             m_boundary_counter = 0;
1164             m_first_from_unknown = false;
1165             m_first_from_unknown_boundary_detected = false;
1166         }
1167 
1168         // check if the passed turn's segment of Linear geometry arrived
1169         // from the inside or the outside of the Areal geometry
1170         template <typename Turn, typename SideStrategy>
calculate_from_inside(Geometry1 const & geometry1,Geometry2 const & geometry2,Turn const & turn,SideStrategy const & side_strategy)1171         static inline bool calculate_from_inside(Geometry1 const& geometry1,
1172                                                  Geometry2 const& geometry2,
1173                                                  Turn const& turn,
1174                                                  SideStrategy const& side_strategy)
1175         {
1176             typedef typename cs_tag<typename Turn::point_type>::type cs_tag;
1177 
1178             if ( turn.operations[op_id].position == overlay::position_front )
1179                 return false;
1180 
1181             typename sub_range_return_type<Geometry1 const>::type
1182                 range1 = sub_range(geometry1, turn.operations[op_id].seg_id);
1183 
1184             typedef detail::normalized_view<Geometry2 const> const range2_type;
1185             typedef typename boost::range_iterator<range2_type>::type range2_iterator;
1186             range2_type range2(sub_range(geometry2, turn.operations[other_op_id].seg_id));
1187 
1188             BOOST_GEOMETRY_ASSERT(boost::size(range1));
1189             std::size_t const s2 = boost::size(range2);
1190             BOOST_GEOMETRY_ASSERT(s2 > 2);
1191             std::size_t const seg_count2 = s2 - 1;
1192 
1193             std::size_t const p_seg_ij = static_cast<std::size_t>(turn.operations[op_id].seg_id.segment_index);
1194             std::size_t const q_seg_ij = static_cast<std::size_t>(turn.operations[other_op_id].seg_id.segment_index);
1195 
1196             BOOST_GEOMETRY_ASSERT(p_seg_ij + 1 < boost::size(range1));
1197             BOOST_GEOMETRY_ASSERT(q_seg_ij + 1 < s2);
1198 
1199             point1_type const& pi = range::at(range1, p_seg_ij);
1200             point2_type const& qi = range::at(range2, q_seg_ij);
1201             point2_type const& qj = range::at(range2, q_seg_ij + 1);
1202             point1_type qi_conv;
1203             geometry::convert(qi, qi_conv);
1204             bool const is_ip_qj = equals::equals_point_point(turn.point, qj);
1205 // TODO: test this!
1206 //            BOOST_GEOMETRY_ASSERT(!equals::equals_point_point(turn.point, pi));
1207 //            BOOST_GEOMETRY_ASSERT(!equals::equals_point_point(turn.point, qi));
1208             point1_type new_pj;
1209             geometry::convert(turn.point, new_pj);
1210 
1211             if ( is_ip_qj )
1212             {
1213                 std::size_t const q_seg_jk = (q_seg_ij + 1) % seg_count2;
1214 // TODO: the following function should return immediately, however the worst case complexity is O(N)
1215 // It would be good to replace it with some O(1) mechanism
1216                 range2_iterator qk_it = find_next_non_duplicated(boost::begin(range2),
1217                                                                  range::pos(range2, q_seg_jk),
1218                                                                  boost::end(range2));
1219 
1220                 // Will this sequence of points be always correct?
1221                 overlay::side_calculator<cs_tag, point1_type, point2_type, SideStrategy>
1222                     side_calc(qi_conv, new_pj, pi, qi, qj, *qk_it, side_strategy);
1223 
1224                 return calculate_from_inside_sides(side_calc);
1225             }
1226             else
1227             {
1228                 point2_type new_qj;
1229                 geometry::convert(turn.point, new_qj);
1230 
1231                 overlay::side_calculator<cs_tag, point1_type, point2_type, SideStrategy>
1232                     side_calc(qi_conv, new_pj, pi, qi, new_qj, qj, side_strategy);
1233 
1234                 return calculate_from_inside_sides(side_calc);
1235             }
1236         }
1237 
1238         template <typename It>
find_next_non_duplicated(It first,It current,It last)1239         static inline It find_next_non_duplicated(It first, It current, It last)
1240         {
1241             BOOST_GEOMETRY_ASSERT( current != last );
1242 
1243             It it = current;
1244 
1245             for ( ++it ; it != last ; ++it )
1246             {
1247                 if ( !equals::equals_point_point(*current, *it) )
1248                     return it;
1249             }
1250 
1251             // if not found start from the beginning
1252             for ( it = first ; it != current ; ++it )
1253             {
1254                 if ( !equals::equals_point_point(*current, *it) )
1255                     return it;
1256             }
1257 
1258             return current;
1259         }
1260 
1261         // calculate inside or outside based on side_calc
1262         // this is simplified version of a check from equal<>
1263         template <typename SideCalc>
calculate_from_inside_sides(SideCalc const & side_calc)1264         static inline bool calculate_from_inside_sides(SideCalc const& side_calc)
1265         {
1266             int const side_pk_p = side_calc.pk_wrt_p1();
1267             int const side_qk_p = side_calc.qk_wrt_p1();
1268             // If they turn to same side (not opposite sides)
1269             if (! overlay::base_turn_handler::opposite(side_pk_p, side_qk_p))
1270             {
1271                 int const side_pk_q2 = side_calc.pk_wrt_q2();
1272                 return side_pk_q2 == -1;
1273             }
1274             else
1275             {
1276                 return side_pk_p == -1;
1277             }
1278         }
1279 
1280     private:
1281         exit_watcher<TurnInfo, op_id> m_exit_watcher;
1282         segment_watcher<same_single> m_seg_watcher;
1283         TurnInfo * m_previous_turn_ptr;
1284         overlay::operation_type m_previous_operation;
1285         unsigned m_boundary_counter;
1286         bool m_interior_detected;
1287         const segment_identifier * m_first_interior_other_id_ptr;
1288         bool m_first_from_unknown;
1289         bool m_first_from_unknown_boundary_detected;
1290     };
1291 
1292     // call analyser.apply() for each turn in range
1293     // IMPORTANT! The analyser is also called for the end iterator - last
1294     template <typename Result,
1295               typename TurnIt,
1296               typename Analyser,
1297               typename Geometry,
1298               typename OtherGeometry,
1299               typename BoundaryChecker,
1300               typename SideStrategy>
analyse_each_turnboost::geometry::detail::relate::linear_areal1301     static inline void analyse_each_turn(Result & res,
1302                                          Analyser & analyser,
1303                                          TurnIt first, TurnIt last,
1304                                          Geometry const& geometry,
1305                                          OtherGeometry const& other_geometry,
1306                                          BoundaryChecker const& boundary_checker,
1307                                          SideStrategy const& side_strategy)
1308     {
1309         if ( first == last )
1310             return;
1311 
1312         for ( TurnIt it = first ; it != last ; ++it )
1313         {
1314             analyser.apply(res, it,
1315                            geometry, other_geometry,
1316                            boundary_checker,
1317                            side_strategy);
1318 
1319             if ( BOOST_GEOMETRY_CONDITION( res.interrupt ) )
1320                 return;
1321         }
1322 
1323         analyser.apply(res, first, last,
1324                        geometry, other_geometry,
1325                        boundary_checker);
1326     }
1327 
1328     // less comparator comparing multi_index then ring_index
1329     // may be used to sort turns by ring
1330     struct less_ring
1331     {
1332         template <typename Turn>
operator ()boost::geometry::detail::relate::linear_areal::less_ring1333         inline bool operator()(Turn const& left, Turn const& right) const
1334         {
1335             return left.operations[1].seg_id.multi_index < right.operations[1].seg_id.multi_index
1336                 || ( left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index
1337                   && left.operations[1].seg_id.ring_index < right.operations[1].seg_id.ring_index );
1338         }
1339     };
1340 
1341     // policy/functor checking if a turn's operation is operation_continue
1342     struct has_boundary_intersection
1343     {
has_boundary_intersectionboost::geometry::detail::relate::linear_areal::has_boundary_intersection1344         has_boundary_intersection()
1345             : result(false) {}
1346 
1347         template <typename Turn>
operator ()boost::geometry::detail::relate::linear_areal::has_boundary_intersection1348         inline void operator()(Turn const& turn)
1349         {
1350             if ( turn.operations[1].operation == overlay::operation_continue )
1351                 result = true;
1352         }
1353 
1354         bool result;
1355     };
1356 
1357     // iterate through the range and search for the different multi_index or ring_index
1358     // also call fun for each turn in the current range
1359     template <typename It, typename Fun>
find_next_ringboost::geometry::detail::relate::linear_areal1360     static inline It find_next_ring(It first, It last, Fun & fun)
1361     {
1362         if ( first == last )
1363             return last;
1364 
1365         signed_size_type const multi_index = first->operations[1].seg_id.multi_index;
1366         signed_size_type const ring_index = first->operations[1].seg_id.ring_index;
1367 
1368         fun(*first);
1369         ++first;
1370 
1371         for ( ; first != last ; ++first )
1372         {
1373             if ( multi_index != first->operations[1].seg_id.multi_index
1374               || ring_index != first->operations[1].seg_id.ring_index )
1375             {
1376                 return first;
1377             }
1378 
1379             fun(*first);
1380         }
1381 
1382         return last;
1383     }
1384 
1385     // analyser which called for turns sorted by seg/distance/operation
1386     // checks if the boundary of Areal geometry really got out
1387     // into the exterior of Linear geometry
1388     template <typename TurnInfo>
1389     class areal_boundary_analyser
1390     {
1391     public:
areal_boundary_analyser()1392         areal_boundary_analyser()
1393             : is_union_detected(false)
1394             , m_previous_turn_ptr(NULL)
1395         {}
1396 
1397         template <typename TurnIt>
apply(TurnIt,TurnIt it,TurnIt last)1398         bool apply(TurnIt /*first*/, TurnIt it, TurnIt last)
1399         {
1400             overlay::operation_type op = it->operations[1].operation;
1401 
1402             if ( it != last )
1403             {
1404                 if ( op != overlay::operation_union
1405                   && op != overlay::operation_continue )
1406                 {
1407                     return true;
1408                 }
1409 
1410                 if ( is_union_detected )
1411                 {
1412                     BOOST_GEOMETRY_ASSERT(m_previous_turn_ptr != NULL);
1413                     if ( !detail::equals::equals_point_point(it->point, m_previous_turn_ptr->point) )
1414                     {
1415                         // break
1416                         return false;
1417                     }
1418                     else if ( op == overlay::operation_continue ) // operation_boundary
1419                     {
1420                         is_union_detected = false;
1421                     }
1422                 }
1423 
1424                 if ( op == overlay::operation_union )
1425                 {
1426                     is_union_detected = true;
1427                     m_previous_turn_ptr = boost::addressof(*it);
1428                 }
1429 
1430                 return true;
1431             }
1432             else
1433             {
1434                 return false;
1435             }
1436         }
1437 
1438         bool is_union_detected;
1439 
1440     private:
1441         const TurnInfo * m_previous_turn_ptr;
1442     };
1443 };
1444 
1445 template <typename Geometry1, typename Geometry2>
1446 struct areal_linear
1447 {
1448     typedef linear_areal<Geometry2, Geometry1, true> linear_areal_type;
1449 
1450     static const bool interruption_enabled = linear_areal_type::interruption_enabled;
1451 
1452     template <typename Result, typename IntersectionStrategy>
applyboost::geometry::detail::relate::areal_linear1453     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2,
1454                              Result & result,
1455                              IntersectionStrategy const& intersection_strategy)
1456     {
1457         linear_areal_type::apply(geometry2, geometry1, result, intersection_strategy);
1458     }
1459 };
1460 
1461 }} // namespace detail::relate
1462 #endif // DOXYGEN_NO_DETAIL
1463 
1464 }} // namespace boost::geometry
1465 
1466 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP
1467