boost/geometry/algorithms/detail/relate/linear_areal.hpp

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