boost/geometry/algorithms/detail/relate/areal_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_AREAL_AREAL_HPP
  15 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP
  16
  17 #include <boost/geometry/core/topological_dimension.hpp>
  18
  19 #include <boost/geometry/util/condition.hpp>
  20 #include <boost/geometry/util/range.hpp>
  21
  22 #include <boost/geometry/algorithms/num_interior_rings.hpp>
  23 #include <boost/geometry/algorithms/detail/point_on_border.hpp>
  24 #include <boost/geometry/algorithms/detail/sub_range.hpp>
  25 #include <boost/geometry/algorithms/detail/single_geometry.hpp>
  26
  27 #include <boost/geometry/algorithms/detail/relate/point_geometry.hpp>
  28 #include <boost/geometry/algorithms/detail/relate/turns.hpp>
  29 #include <boost/geometry/algorithms/detail/relate/boundary_checker.hpp>
  30 #include <boost/geometry/algorithms/detail/relate/follow_helpers.hpp>
  31
  32 namespace boost { namespace geometry
  33 {
  34
  35 #ifndef DOXYGEN_NO_DETAIL
  36 namespace detail { namespace relate {
  37
  38 // WARNING!
  39 // TODO: In the worst case calling this Pred in a loop for MultiPolygon/MultiPolygon may take O(NM)
  40 // Use the rtree in this case!
  41
  42 // may be used to set EI and EB for an Areal geometry for which no turns were generated
  43 template
  44 <
  45     typename OtherAreal,
  46     typename Result,
  47     typename PointInArealStrategy,
  48     bool TransposeResult
  49 >
  50 class no_turns_aa_pred
  51 {
  52 public:
  53     no_turns_aa_pred(OtherAreal const& other_areal,
  54                      Result & res,
  55                      PointInArealStrategy const& point_in_areal_strategy)
  56         : m_result(res)
  57         , m_point_in_areal_strategy(point_in_areal_strategy)
  58         , m_other_areal(other_areal)
  59         , m_flags(0)
  60     {
  61         // check which relations must be analysed
  62
  63         if ( ! may_update<interior, interior, '2', TransposeResult>(m_result)
  64           && ! may_update<boundary, interior, '1', TransposeResult>(m_result)
  65           && ! may_update<exterior, interior, '2', TransposeResult>(m_result) )
  66         {
  67             m_flags |= 1;
  68         }
  69
  70         if ( ! may_update<interior, exterior, '2', TransposeResult>(m_result)
  71           && ! may_update<boundary, exterior, '1', TransposeResult>(m_result) )
  72         {
  73             m_flags |= 2;
  74         }
  75     }
  76
  77     template <typename Areal>
  78     bool operator()(Areal const& areal)
  79     {
  80         using detail::within::point_in_geometry;
  81
  82         // if those flags are set nothing will change
  83         if ( m_flags == 3 )
  84         {
  85             return false;
  86         }
  87
  88         typedef typename geometry::point_type<Areal>::type point_type;
  89         point_type pt;
  90         bool const ok = boost::geometry::point_on_border(pt, areal);
  91
  92         // TODO: for now ignore, later throw an exception?
  93         if ( !ok )
  94         {
  95             return true;
  96         }
  97
  98         // check if the areal is inside the other_areal
  99         // TODO: This is O(N)
 100         // Run in a loop O(NM) - optimize!
 101         int const pig = point_in_geometry(pt,
 102                                           m_other_areal,
 103                                           m_point_in_areal_strategy);
 104         //BOOST_GEOMETRY_ASSERT( pig != 0 );
 105
 106         // inside
 107         if ( pig > 0 )
 108         {
 109             update<interior, interior, '2', TransposeResult>(m_result);
 110             update<boundary, interior, '1', TransposeResult>(m_result);
 111             update<exterior, interior, '2', TransposeResult>(m_result);
 112             m_flags |= 1;
 113
 114             // TODO: OPTIMIZE!
 115             // Only the interior rings of other ONE single geometry must be checked
 116             // NOT all geometries
 117
 118             // Check if any interior ring is outside
 119             ring_identifier ring_id(0, -1, 0);
 120             std::size_t const irings_count = geometry::num_interior_rings(areal);
 121             for ( ; static_cast<std::size_t>(ring_id.ring_index) < irings_count ;
 122                     ++ring_id.ring_index )
 123             {
 124                 typename detail::sub_range_return_type<Areal const>::type
 125                     range_ref = detail::sub_range(areal, ring_id);
 126
 127                 if ( boost::empty(range_ref) )
 128                 {
 129                     // TODO: throw exception?
 130                     continue; // ignore
 131                 }
 132
 133                 // TODO: O(N)
 134                 // Optimize!
 135                 int const hpig = point_in_geometry(range::front(range_ref),
 136                                                    m_other_areal,
 137                                                    m_point_in_areal_strategy);
 138
 139                 // hole outside
 140                 if ( hpig < 0 )
 141                 {
 142                     update<interior, exterior, '2', TransposeResult>(m_result);
 143                     update<boundary, exterior, '1', TransposeResult>(m_result);
 144                     m_flags |= 2;
 145                     break;
 146                 }
 147             }
 148         }
 149         // outside
 150         else
 151         {
 152             update<interior, exterior, '2', TransposeResult>(m_result);
 153             update<boundary, exterior, '1', TransposeResult>(m_result);
 154             m_flags |= 2;
 155
 156             // Check if any interior ring is inside
 157             ring_identifier ring_id(0, -1, 0);
 158             std::size_t const irings_count = geometry::num_interior_rings(areal);
 159             for ( ; static_cast<std::size_t>(ring_id.ring_index) < irings_count ;
 160                     ++ring_id.ring_index )
 161             {
 162                 typename detail::sub_range_return_type<Areal const>::type
 163                     range_ref = detail::sub_range(areal, ring_id);
 164
 165                 if ( boost::empty(range_ref) )
 166                 {
 167                     // TODO: throw exception?
 168                     continue; // ignore
 169                 }
 170
 171                 // TODO: O(N)
 172                 // Optimize!
 173                 int const hpig = point_in_geometry(range::front(range_ref),
 174                                                    m_other_areal,
 175                                                    m_point_in_areal_strategy);
 176
 177                 // hole inside
 178                 if ( hpig > 0 )
 179                 {
 180                     update<interior, interior, '2', TransposeResult>(m_result);
 181                     update<boundary, interior, '1', TransposeResult>(m_result);
 182                     update<exterior, interior, '2', TransposeResult>(m_result);
 183                     m_flags |= 1;
 184                     break;
 185                 }
 186             }
 187         }
 188
 189         return m_flags != 3 && !m_result.interrupt;
 190     }
 191
 192 private:
 193     Result & m_result;
 194     PointInArealStrategy const& m_point_in_areal_strategy;
 195     OtherAreal const& m_other_areal;
 196     int m_flags;
 197 };
 198
 199 // The implementation of an algorithm calculating relate() for A/A
 200 template <typename Geometry1, typename Geometry2>
 201 struct areal_areal
 202 {
 203     // check Linear / Areal
 204     BOOST_STATIC_ASSERT(topological_dimension<Geometry1>::value == 2
 205                      && topological_dimension<Geometry2>::value == 2);
 206
 207     static const bool interruption_enabled = true;
 208
 209     typedef typename geometry::point_type<Geometry1>::type point1_type;
 210     typedef typename geometry::point_type<Geometry2>::type point2_type;
 211
 212     template <typename Result, typename IntersectionStrategy>
 213     static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2,
 214                              Result & result,
 215                              IntersectionStrategy const& intersection_strategy)
 216     {
 217 // TODO: If Areal geometry may have infinite size, change the following line:
 218
 219         // The result should be FFFFFFFFF
 220         relate::set<exterior, exterior, result_dimension<Geometry2>::value>(result);// FFFFFFFFd, d in [1,9] or T
 221
 222         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 223             return;
 224
 225         // get and analyse turns
 226         typedef typename turns::get_turns
 227             <
 228                 Geometry1, Geometry2
 229             >::template turn_info_type<IntersectionStrategy>::type turn_type;
 230         std::vector<turn_type> turns;
 231
 232         interrupt_policy_areal_areal<Result> interrupt_policy(geometry1, geometry2, result);
 233
 234         turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy, intersection_strategy);
 235         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 236             return;
 237
 238         typedef typename IntersectionStrategy::cs_tag cs_tag;
 239
 240         typedef typename IntersectionStrategy::template point_in_geometry_strategy
 241             <
 242                 Geometry1, Geometry2
 243             >::type point_in_areal_strategy12_type;
 244         point_in_areal_strategy12_type point_in_areal_strategy12
 245             = intersection_strategy.template get_point_in_geometry_strategy<Geometry1, Geometry2>();
 246         typedef typename IntersectionStrategy::template point_in_geometry_strategy
 247             <
 248                 Geometry2, Geometry1
 249             >::type point_in_areal_strategy21_type;
 250         point_in_areal_strategy21_type point_in_areal_strategy21
 251             = intersection_strategy.template get_point_in_geometry_strategy<Geometry2, Geometry1>();
 252
 253         no_turns_aa_pred<Geometry2, Result, point_in_areal_strategy12_type, false>
 254             pred1(geometry2, result, point_in_areal_strategy12);
 255         for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1);
 256         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 257             return;
 258
 259         no_turns_aa_pred<Geometry1, Result, point_in_areal_strategy21_type, true>
 260             pred2(geometry1, result, point_in_areal_strategy21);
 261         for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2);
 262         if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 263             return;
 264
 265         if ( turns.empty() )
 266             return;
 267
 268         if ( may_update<interior, interior, '2'>(result)
 269           || may_update<interior, exterior, '2'>(result)
 270           || may_update<boundary, interior, '1'>(result)
 271           || may_update<boundary, exterior, '1'>(result)
 272           || may_update<exterior, interior, '2'>(result) )
 273         {
 274             // sort turns
 275             typedef turns::less<0, turns::less_op_areal_areal<0>, cs_tag> less;
 276             std::sort(turns.begin(), turns.end(), less());
 277
 278             /*if ( may_update<interior, exterior, '2'>(result)
 279               || may_update<boundary, exterior, '1'>(result)
 280               || may_update<boundary, interior, '1'>(result)
 281               || may_update<exterior, interior, '2'>(result) )*/
 282             {
 283                 // analyse sorted turns
 284                 turns_analyser<turn_type, 0> analyser;
 285                 analyse_each_turn(result, analyser, turns.begin(), turns.end(),
 286                                   point_in_areal_strategy12.get_equals_point_point_strategy());
 287
 288                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 289                     return;
 290             }
 291
 292             if ( may_update<interior, interior, '2'>(result)
 293               || may_update<interior, exterior, '2'>(result)
 294               || may_update<boundary, interior, '1'>(result)
 295               || may_update<boundary, exterior, '1'>(result)
 296               || may_update<exterior, interior, '2'>(result) )
 297             {
 298                 // analyse rings for which turns were not generated
 299                 // or only i/i or u/u was generated
 300                 uncertain_rings_analyser<0, Result, Geometry1, Geometry2, point_in_areal_strategy12_type>
 301                     rings_analyser(result, geometry1, geometry2, point_in_areal_strategy12);
 302                 analyse_uncertain_rings<0>::apply(rings_analyser, turns.begin(), turns.end());
 303
 304                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 305                     return;
 306             }
 307         }
 308
 309         if ( may_update<interior, interior, '2', true>(result)
 310           || may_update<interior, exterior, '2', true>(result)
 311           || may_update<boundary, interior, '1', true>(result)
 312           || may_update<boundary, exterior, '1', true>(result)
 313           || may_update<exterior, interior, '2', true>(result) )
 314         {
 315             // sort turns
 316             typedef turns::less<1, turns::less_op_areal_areal<1>, cs_tag> less;
 317             std::sort(turns.begin(), turns.end(), less());
 318
 319             /*if ( may_update<interior, exterior, '2', true>(result)
 320               || may_update<boundary, exterior, '1', true>(result)
 321               || may_update<boundary, interior, '1', true>(result)
 322               || may_update<exterior, interior, '2', true>(result) )*/
 323             {
 324                 // analyse sorted turns
 325                 turns_analyser<turn_type, 1> analyser;
 326                 analyse_each_turn(result, analyser, turns.begin(), turns.end(),
 327                                   point_in_areal_strategy21.get_equals_point_point_strategy());
 328
 329                 if ( BOOST_GEOMETRY_CONDITION(result.interrupt) )
 330                     return;
 331             }
 332
 333             if ( may_update<interior, interior, '2', true>(result)
 334               || may_update<interior, exterior, '2', true>(result)
 335               || may_update<boundary, interior, '1', true>(result)
 336               || may_update<boundary, exterior, '1', true>(result)
 337               || may_update<exterior, interior, '2', true>(result) )
 338             {
 339                 // analyse rings for which turns were not generated
 340                 // or only i/i or u/u was generated
 341                 uncertain_rings_analyser<1, Result, Geometry2, Geometry1, point_in_areal_strategy21_type>
 342                     rings_analyser(result, geometry2, geometry1, point_in_areal_strategy21);
 343                 analyse_uncertain_rings<1>::apply(rings_analyser, turns.begin(), turns.end());
 344
 345                 //if ( result.interrupt )
 346                 //    return;
 347             }
 348         }
 349     }
 350
 351     // interrupt policy which may be passed to get_turns to interrupt the analysis
 352     // based on the info in the passed result/mask
 353     template <typename Result>
 354     class interrupt_policy_areal_areal
 355     {
 356     public:
 357         static bool const enabled = true;
 358
 359         interrupt_policy_areal_areal(Geometry1 const& geometry1,
 360                                      Geometry2 const& geometry2,
 361                                      Result & result)
 362             : m_result(result)
 363             , m_geometry1(geometry1)
 364             , m_geometry2(geometry2)
 365         {}
 366
 367         template <typename Range>
 368         inline bool apply(Range const& turns)
 369         {
 370             typedef typename boost::range_iterator<Range const>::type iterator;
 371
 372             for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it)
 373             {
 374                 per_turn<0>(*it);
 375                 per_turn<1>(*it);
 376             }
 377
 378             return m_result.interrupt;
 379         }
 380
 381     private:
 382         template <std::size_t OpId, typename Turn>
 383         inline void per_turn(Turn const& turn)
 384         {
 385             //static const std::size_t other_op_id = (OpId + 1) % 2;
 386             static const bool transpose_result = OpId != 0;
 387
 388             overlay::operation_type const op = turn.operations[OpId].operation;
 389
 390             if ( op == overlay::operation_union )
 391             {
 392                 // ignore u/u
 393                 /*if ( turn.operations[other_op_id].operation != overlay::operation_union )
 394                 {
 395                     update<interior, exterior, '2', transpose_result>(m_result);
 396                     update<boundary, exterior, '1', transpose_result>(m_result);
 397                 }*/
 398
 399                 update<boundary, boundary, '0', transpose_result>(m_result);
 400             }
 401             else if ( op == overlay::operation_intersection )
 402             {
 403                 // ignore i/i
 404                 /*if ( turn.operations[other_op_id].operation != overlay::operation_intersection )
 405                 {
 406                     // not correct e.g. for G1 touching G2 in a point where a hole is touching the exterior ring
 407                     // in this case 2 turns i/... and u/u will be generated for this IP
 408                     //update<interior, interior, '2', transpose_result>(m_result);
 409
 410                     //update<boundary, interior, '1', transpose_result>(m_result);
 411                 }*/
 412
 413                 update<boundary, boundary, '0', transpose_result>(m_result);
 414             }
 415             else if ( op == overlay::operation_continue )
 416             {
 417                 update<boundary, boundary, '1', transpose_result>(m_result);
 418                 update<interior, interior, '2', transpose_result>(m_result);
 419             }
 420             else if ( op == overlay::operation_blocked )
 421             {
 422                 update<boundary, boundary, '1', transpose_result>(m_result);
 423                 update<interior, exterior, '2', transpose_result>(m_result);
 424             }
 425         }
 426
 427         Result & m_result;
 428         Geometry1 const& m_geometry1;
 429         Geometry2 const& m_geometry2;
 430     };
 431
 432     // This analyser should be used like Input or SinglePass Iterator
 433     // IMPORTANT! It should be called also for the end iterator - last
 434     template <typename TurnInfo, std::size_t OpId>
 435     class turns_analyser
 436     {
 437         typedef typename TurnInfo::point_type turn_point_type;
 438
 439         static const std::size_t op_id = OpId;
 440         static const std::size_t other_op_id = (OpId + 1) % 2;
 441         static const bool transpose_result = OpId != 0;
 442
 443     public:
 444         turns_analyser()
 445             : m_previous_turn_ptr(0)
 446             , m_previous_operation(overlay::operation_none)
 447             , m_enter_detected(false)
 448             , m_exit_detected(false)
 449         {}
 450
 451         template <typename Result, typename TurnIt, typename EqPPStrategy>
 452         void apply(Result & result, TurnIt it, EqPPStrategy const& strategy)
 453         {
 454             //BOOST_GEOMETRY_ASSERT( it != last );
 455
 456             overlay::operation_type const op = it->operations[op_id].operation;
 457
 458             if ( op != overlay::operation_union
 459               && op != overlay::operation_intersection
 460               && op != overlay::operation_blocked
 461               && op != overlay::operation_continue )
 462             {
 463                 return;
 464             }
 465
 466             segment_identifier const& seg_id = it->operations[op_id].seg_id;
 467             //segment_identifier const& other_id = it->operations[other_op_id].seg_id;
 468
 469             const bool first_in_range = m_seg_watcher.update(seg_id);
 470
 471             if ( m_previous_turn_ptr )
 472             {
 473                 if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
 474                 {
 475                     // real exit point - may be multiple
 476                     if ( first_in_range
 477                       || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it, strategy) )
 478                     {
 479                         update_exit(result);
 480                         m_exit_detected = false;
 481                     }
 482                     // fake exit point, reset state
 483                     else if ( op != overlay::operation_union )
 484                     {
 485                         m_exit_detected = false;
 486                     }
 487                 }
 488                 /*else*/
 489                 if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
 490                 {
 491                     // real entry point
 492                     if ( first_in_range
 493                       || ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it, strategy) )
 494                     {
 495                         update_enter(result);
 496                         m_enter_detected = false;
 497                     }
 498                     // fake entry point, reset state
 499                     else if ( op != overlay::operation_intersection )
 500                     {
 501                         m_enter_detected = false;
 502                     }
 503                 }
 504             }
 505
 506             if ( op == overlay::operation_union )
 507             {
 508                 // already set in interrupt policy
 509                 //update<boundary, boundary, '0', transpose_result>(m_result);
 510
 511                 // ignore u/u
 512                 //if ( it->operations[other_op_id].operation != overlay::operation_union )
 513                 {
 514                     m_exit_detected = true;
 515                 }
 516             }
 517             else if ( op == overlay::operation_intersection )
 518             {
 519                 // ignore i/i
 520                 if ( it->operations[other_op_id].operation != overlay::operation_intersection )
 521                 {
 522                     // this was set in the interrupt policy but it was wrong
 523                     // also here it's wrong since it may be a fake entry point
 524                     //update<interior, interior, '2', transpose_result>(result);
 525
 526                     // already set in interrupt policy
 527                     //update<boundary, boundary, '0', transpose_result>(result);
 528                     m_enter_detected = true;
 529                 }
 530             }
 531             else if ( op == overlay::operation_blocked )
 532             {
 533                 // already set in interrupt policy
 534             }
 535             else // if ( op == overlay::operation_continue )
 536             {
 537                 // already set in interrupt policy
 538             }
 539
 540             // store ref to previously analysed (valid) turn
 541             m_previous_turn_ptr = boost::addressof(*it);
 542             // and previously analysed (valid) operation
 543             m_previous_operation = op;
 544         }
 545
 546         // it == last
 547         template <typename Result>
 548         void apply(Result & result)
 549         {
 550             //BOOST_GEOMETRY_ASSERT( first != last );
 551
 552             if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ )
 553             {
 554                 update_exit(result);
 555                 m_exit_detected = false;
 556             }
 557
 558             if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ )
 559             {
 560                 update_enter(result);
 561                 m_enter_detected = false;
 562             }
 563         }
 564
 565         template <typename Result>
 566         static inline void update_exit(Result & result)
 567         {
 568             update<interior, exterior, '2', transpose_result>(result);
 569             update<boundary, exterior, '1', transpose_result>(result);
 570         }
 571
 572         template <typename Result>
 573         static inline void update_enter(Result & result)
 574         {
 575             update<interior, interior, '2', transpose_result>(result);
 576             update<boundary, interior, '1', transpose_result>(result);
 577             update<exterior, interior, '2', transpose_result>(result);
 578         }
 579
 580     private:
 581         segment_watcher<same_ring> m_seg_watcher;
 582         TurnInfo * m_previous_turn_ptr;
 583         overlay::operation_type m_previous_operation;
 584         bool m_enter_detected;
 585         bool m_exit_detected;
 586     };
 587
 588     // call analyser.apply() for each turn in range
 589     // IMPORTANT! The analyser is also called for the end iterator - last
 590     template
 591     <
 592         typename Result,
 593         typename Analyser,
 594         typename TurnIt,
 595         typename EqPPStrategy
 596     >
 597     static inline void analyse_each_turn(Result & res,
 598                                          Analyser & analyser,
 599                                          TurnIt first, TurnIt last,
 600                                          EqPPStrategy const& strategy)
 601     {
 602         if ( first == last )
 603             return;
 604
 605         for ( TurnIt it = first ; it != last ; ++it )
 606         {
 607             analyser.apply(res, it, strategy);
 608
 609             if ( BOOST_GEOMETRY_CONDITION(res.interrupt) )
 610                 return;
 611         }
 612
 613         analyser.apply(res);
 614     }
 615
 616     template
 617     <
 618         std::size_t OpId,
 619         typename Result,
 620         typename Geometry,
 621         typename OtherGeometry,
 622         typename PointInArealStrategy
 623     >
 624     class uncertain_rings_analyser
 625     {
 626         static const bool transpose_result = OpId != 0;
 627         static const int other_id = (OpId + 1) % 2;
 628
 629     public:
 630         inline uncertain_rings_analyser(Result & result,
 631                                         Geometry const& geom,
 632                                         OtherGeometry const& other_geom,
 633                                         PointInArealStrategy const& point_in_areal_strategy)
 634             : geometry(geom)
 635             , other_geometry(other_geom)
 636             , interrupt(result.interrupt) // just in case, could be false as well
 637             , m_result(result)
 638             , m_point_in_areal_strategy(point_in_areal_strategy)
 639             , m_flags(0)
 640         {
 641             // check which relations must be analysed
 642             // NOTE: 1 and 4 could probably be connected
 643
 644             if ( ! may_update<interior, interior, '2', transpose_result>(m_result) )
 645             {
 646                 m_flags |= 1;
 647             }
 648
 649             if ( ! may_update<interior, exterior, '2', transpose_result>(m_result)
 650               && ! may_update<boundary, exterior, '1', transpose_result>(m_result) )
 651             {
 652                 m_flags |= 2;
 653             }
 654
 655             if ( ! may_update<boundary, interior, '1', transpose_result>(m_result)
 656               && ! may_update<exterior, interior, '2', transpose_result>(m_result) )
 657             {
 658                 m_flags |= 4;
 659             }
 660         }
 661
 662         inline void no_turns(segment_identifier const& seg_id)
 663         {
 664             // if those flags are set nothing will change
 665             if ( m_flags == 7 )
 666             {
 667                 return;
 668             }
 669
 670             typename detail::sub_range_return_type<Geometry const>::type
 671                 range_ref = detail::sub_range(geometry, seg_id);
 672
 673             if ( boost::empty(range_ref) )
 674             {
 675                 // TODO: throw an exception?
 676                 return; // ignore
 677             }
 678
 679             // TODO: possible optimization
 680             // if the range is an interior ring we may use other IPs generated for this single geometry
 681             // to know which other single geometries should be checked
 682
 683             // TODO: optimize! e.g. use spatial index
 684             // O(N) - running it in a loop gives O(NM)
 685             using detail::within::point_in_geometry;
 686             int const pig = point_in_geometry(range::front(range_ref),
 687                                               other_geometry,
 688                                               m_point_in_areal_strategy);
 689
 690             //BOOST_GEOMETRY_ASSERT(pig != 0);
 691             if ( pig > 0 )
 692             {
 693                 update<interior, interior, '2', transpose_result>(m_result);
 694                 m_flags |= 1;
 695
 696                 update<boundary, interior, '1', transpose_result>(m_result);
 697                 update<exterior, interior, '2', transpose_result>(m_result);
 698                 m_flags |= 4;
 699             }
 700             else
 701             {
 702                 update<boundary, exterior, '1', transpose_result>(m_result);
 703                 update<interior, exterior, '2', transpose_result>(m_result);
 704                 m_flags |= 2;
 705             }
 706
 707 // TODO: break if all things are set
 708 // also some of them could be checked outside, before the analysis
 709 // In this case we shouldn't relay just on dummy flags
 710 // Flags should be initialized with proper values
 711 // or the result should be checked directly
 712 // THIS IS ALSO TRUE FOR OTHER ANALYSERS! in L/L and L/A
 713
 714             interrupt = m_flags == 7 || m_result.interrupt;
 715         }
 716
 717         template <typename TurnIt>
 718         inline void turns(TurnIt first, TurnIt last)
 719         {
 720             // if those flags are set nothing will change
 721             if ( (m_flags & 6) == 6 )
 722             {
 723                 return;
 724             }
 725
 726             bool found_ii = false;
 727             bool found_uu = false;
 728
 729             for ( TurnIt it = first ; it != last ; ++it )
 730             {
 731                 if ( it->operations[0].operation == overlay::operation_intersection
 732                   && it->operations[1].operation == overlay::operation_intersection )
 733                 {
 734                     found_ii = true;
 735                 }
 736                 else if ( it->operations[0].operation == overlay::operation_union
 737                        && it->operations[1].operation == overlay::operation_union )
 738                 {
 739                     found_uu = true;
 740                 }
 741                 else // ignore
 742                 {
 743                     return; // don't interrupt
 744                 }
 745             }
 746
 747             // only i/i was generated for this ring
 748             if ( found_ii )
 749             {
 750                 update<interior, interior, '2', transpose_result>(m_result);
 751                 m_flags |= 1;
 752
 753                 //update<boundary, boundary, '0', transpose_result>(m_result);
 754
 755                 update<boundary, interior, '1', transpose_result>(m_result);
 756                 update<exterior, interior, '2', transpose_result>(m_result);
 757                 m_flags |= 4;
 758             }
 759
 760             // only u/u was generated for this ring
 761             if ( found_uu )
 762             {
 763                 update<boundary, exterior, '1', transpose_result>(m_result);
 764                 update<interior, exterior, '2', transpose_result>(m_result);
 765                 m_flags |= 2;
 766             }
 767
 768             interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future
 769         }
 770
 771         Geometry const& geometry;
 772         OtherGeometry const& other_geometry;
 773         bool interrupt;
 774
 775     private:
 776         Result & m_result;
 777         PointInArealStrategy const& m_point_in_areal_strategy;
 778         int m_flags;
 779     };
 780
 781     template <std::size_t OpId>
 782     class analyse_uncertain_rings
 783     {
 784     public:
 785         template <typename Analyser, typename TurnIt>
 786         static inline void apply(Analyser & analyser, TurnIt first, TurnIt last)
 787         {
 788             if ( first == last )
 789                 return;
 790
 791             for_preceding_rings(analyser, *first);
 792             //analyser.per_turn(*first);
 793
 794             TurnIt prev = first;
 795             for ( ++first ; first != last ; ++first, ++prev )
 796             {
 797                 // same multi
 798                 if ( prev->operations[OpId].seg_id.multi_index
 799                   == first->operations[OpId].seg_id.multi_index )
 800                 {
 801                     // same ring
 802                     if ( prev->operations[OpId].seg_id.ring_index
 803                       == first->operations[OpId].seg_id.ring_index )
 804                     {
 805                         //analyser.per_turn(*first);
 806                     }
 807                     // same multi, next ring
 808                     else
 809                     {
 810                         //analyser.end_ring(*prev);
 811                         analyser.turns(prev, first);
 812
 813                         //if ( prev->operations[OpId].seg_id.ring_index + 1
 814                         //   < first->operations[OpId].seg_id.ring_index)
 815                         {
 816                             for_no_turns_rings(analyser,
 817                                                *first,
 818                                                prev->operations[OpId].seg_id.ring_index + 1,
 819                                                first->operations[OpId].seg_id.ring_index);
 820                         }
 821
 822                         //analyser.per_turn(*first);
 823                     }
 824                 }
 825                 // next multi
 826                 else
 827                 {
 828                     //analyser.end_ring(*prev);
 829                     analyser.turns(prev, first);
 830                     for_following_rings(analyser, *prev);
 831                     for_preceding_rings(analyser, *first);
 832                     //analyser.per_turn(*first);
 833                 }
 834
 835                 if ( analyser.interrupt )
 836                 {
 837                     return;
 838                 }
 839             }
 840
 841             //analyser.end_ring(*prev);
 842             analyser.turns(prev, first); // first == last
 843             for_following_rings(analyser, *prev);
 844         }
 845
 846     private:
 847         template <typename Analyser, typename Turn>
 848         static inline void for_preceding_rings(Analyser & analyser, Turn const& turn)
 849         {
 850             segment_identifier const& seg_id = turn.operations[OpId].seg_id;
 851
 852             for_no_turns_rings(analyser, turn, -1, seg_id.ring_index);
 853         }
 854
 855         template <typename Analyser, typename Turn>
 856         static inline void for_following_rings(Analyser & analyser, Turn const& turn)
 857         {
 858             segment_identifier const& seg_id = turn.operations[OpId].seg_id;
 859
 860             signed_size_type
 861                 count = boost::numeric_cast<signed_size_type>(
 862                             geometry::num_interior_rings(
 863                                 detail::single_geometry(analyser.geometry, seg_id)));
 864
 865             for_no_turns_rings(analyser, turn, seg_id.ring_index + 1, count);
 866         }
 867
 868         template <typename Analyser, typename Turn>
 869         static inline void for_no_turns_rings(Analyser & analyser,
 870                                               Turn const& turn,
 871                                               signed_size_type first,
 872                                               signed_size_type last)
 873         {
 874             segment_identifier seg_id = turn.operations[OpId].seg_id;
 875
 876             for ( seg_id.ring_index = first ; seg_id.ring_index < last ; ++seg_id.ring_index )
 877             {
 878                 analyser.no_turns(seg_id);
 879             }
 880         }
 881     };
 882 };
 883
 884 }} // namespace detail::relate
 885 #endif // DOXYGEN_NO_DETAIL
 886
 887 }} // namespace boost::geometry
 888
 889 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP