// Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. // This file was modified by Oracle on 2013, 2014. // Modifications copyright (c) 2013-2014 Oracle and/or its affiliates. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP #include #include #include #include #include #include #include #include #include #include namespace boost { namespace geometry { #ifndef DOXYGEN_NO_DETAIL namespace detail { namespace relate { // WARNING! // TODO: In the worst case calling this Pred in a loop for MultiPolygon/MultiPolygon may take O(NM) // Use the rtree in this case! // may be used to set EI and EB for an Areal geometry for which no turns were generated template class no_turns_aa_pred { public: no_turns_aa_pred(OtherAreal const& other_areal, Result & res) : m_result(res) , m_other_areal(other_areal) , m_flags(0) { // check which relations must be analysed if ( ! may_update(m_result) && ! may_update(m_result) && ! may_update(m_result) ) { m_flags |= 1; } if ( ! may_update(m_result) && ! may_update(m_result) ) { m_flags |= 2; } } template bool operator()(Areal const& areal) { // if those flags are set nothing will change if ( m_flags == 3 ) { return false; } typedef typename geometry::point_type::type point_type; point_type pt; bool const ok = boost::geometry::point_on_border(pt, areal); // TODO: for now ignore, later throw an exception? if ( !ok ) { return true; } // check if the areal is inside the other_areal // TODO: This is O(N) // Run in a loop O(NM) - optimize! int const pig = detail::within::point_in_geometry(pt, m_other_areal); //BOOST_ASSERT( pig != 0 ); // inside if ( pig > 0 ) { update(m_result); update(m_result); update(m_result); m_flags |= 1; // TODO: OPTIMIZE! // Only the interior rings of other ONE single geometry must be checked // NOT all geometries // Check if any interior ring is outside ring_identifier ring_id(0, -1, 0); int const irings_count = boost::numeric_cast( geometry::num_interior_rings(areal) ); for ( ; ring_id.ring_index < irings_count ; ++ring_id.ring_index ) { typename detail::sub_range_return_type::type range_ref = detail::sub_range(areal, ring_id); if ( boost::empty(range_ref) ) { // TODO: throw exception? continue; // ignore } // TODO: O(N) // Optimize! int const hpig = detail::within::point_in_geometry(range::front(range_ref), m_other_areal); // hole outside if ( hpig < 0 ) { update(m_result); update(m_result); m_flags |= 2; break; } } } // outside else { update(m_result); update(m_result); m_flags |= 2; // Check if any interior ring is inside ring_identifier ring_id(0, -1, 0); int const irings_count = boost::numeric_cast( geometry::num_interior_rings(areal) ); for ( ; ring_id.ring_index < irings_count ; ++ring_id.ring_index ) { typename detail::sub_range_return_type::type range_ref = detail::sub_range(areal, ring_id); if ( boost::empty(range_ref) ) { // TODO: throw exception? continue; // ignore } // TODO: O(N) // Optimize! int const hpig = detail::within::point_in_geometry(range::front(range_ref), m_other_areal); // hole inside if ( hpig > 0 ) { update(m_result); update(m_result); update(m_result); m_flags |= 1; break; } } } return m_flags != 3 && !m_result.interrupt; } private: Result & m_result; OtherAreal const& m_other_areal; int m_flags; }; // The implementation of an algorithm calculating relate() for A/A template struct areal_areal { // check Linear / Areal BOOST_STATIC_ASSERT(topological_dimension::value == 2 && topological_dimension::value == 2); static const bool interruption_enabled = true; typedef typename geometry::point_type::type point1_type; typedef typename geometry::point_type::type point2_type; template static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result) { // TODO: If Areal geometry may have infinite size, change the following line: // The result should be FFFFFFFFF relate::set::value>(result);// FFFFFFFFd, d in [1,9] or T if ( result.interrupt ) return; // get and analyse turns typedef typename turns::get_turns::turn_info turn_type; std::vector turns; interrupt_policy_areal_areal interrupt_policy(geometry1, geometry2, result); turns::get_turns::apply(turns, geometry1, geometry2, interrupt_policy); if ( result.interrupt ) return; no_turns_aa_pred pred1(geometry2, result); for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1); if ( result.interrupt ) return; no_turns_aa_pred pred2(geometry1, result); for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2); if ( result.interrupt ) return; if ( turns.empty() ) return; if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) || may_update(result) ) { // sort turns typedef turns::less<0, turns::less_op_areal_areal<0> > less; std::sort(turns.begin(), turns.end(), less()); /*if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) )*/ { // analyse sorted turns turns_analyser analyser; analyse_each_turn(result, analyser, turns.begin(), turns.end()); if ( result.interrupt ) return; } if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) || may_update(result) ) { // analyse rings for which turns were not generated // or only i/i or u/u was generated uncertain_rings_analyser<0, Result, Geometry1, Geometry2> rings_analyser(result, geometry1, geometry2); analyse_uncertain_rings<0>::apply(rings_analyser, turns.begin(), turns.end()); if ( result.interrupt ) return; } } if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) || may_update(result) ) { // sort turns typedef turns::less<1, turns::less_op_areal_areal<1> > less; std::sort(turns.begin(), turns.end(), less()); /*if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) )*/ { // analyse sorted turns turns_analyser analyser; analyse_each_turn(result, analyser, turns.begin(), turns.end()); if ( result.interrupt ) return; } if ( may_update(result) || may_update(result) || may_update(result) || may_update(result) || may_update(result) ) { // analyse rings for which turns were not generated // or only i/i or u/u was generated uncertain_rings_analyser<1, Result, Geometry2, Geometry1> rings_analyser(result, geometry2, geometry1); analyse_uncertain_rings<1>::apply(rings_analyser, turns.begin(), turns.end()); //if ( result.interrupt ) // return; } } } // interrupt policy which may be passed to get_turns to interrupt the analysis // based on the info in the passed result/mask template class interrupt_policy_areal_areal { public: static bool const enabled = true; interrupt_policy_areal_areal(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result) : m_result(result) , m_geometry1(geometry1) , m_geometry2(geometry2) {} template inline bool apply(Range const& turns) { typedef typename boost::range_iterator::type iterator; for (iterator it = boost::begin(turns) ; it != boost::end(turns) ; ++it) { per_turn<0>(*it); per_turn<1>(*it); } return m_result.interrupt; } private: template inline void per_turn(Turn const& turn) { static const std::size_t other_op_id = (OpId + 1) % 2; static const bool transpose_result = OpId != 0; overlay::operation_type const op = turn.operations[OpId].operation; if ( op == overlay::operation_union ) { // ignore u/u /*if ( turn.operations[other_op_id].operation != overlay::operation_union ) { update(m_result); update(m_result); }*/ update(m_result); } else if ( op == overlay::operation_intersection ) { // ignore i/i if ( turn.operations[other_op_id].operation != overlay::operation_intersection ) { update(m_result); //update(m_result); } update(m_result); } else if ( op == overlay::operation_continue ) { update(m_result); update(m_result); } else if ( op == overlay::operation_blocked ) { update(m_result); update(m_result); } } Result & m_result; Geometry1 const& m_geometry1; Geometry2 const& m_geometry2; }; // This analyser should be used like Input or SinglePass Iterator // IMPORTANT! It should be called also for the end iterator - last template class turns_analyser { typedef typename TurnInfo::point_type turn_point_type; static const std::size_t op_id = OpId; static const std::size_t other_op_id = (OpId + 1) % 2; static const bool transpose_result = OpId != 0; public: turns_analyser() : m_previous_turn_ptr(0) , m_previous_operation(overlay::operation_none) , m_enter_detected(false) , m_exit_detected(false) {} template void apply(Result & result, TurnIt it) { //BOOST_ASSERT( it != last ); overlay::operation_type const op = it->operations[op_id].operation; if ( op != overlay::operation_union && op != overlay::operation_intersection && op != overlay::operation_blocked && op != overlay::operation_continue ) { return; } segment_identifier const& seg_id = it->operations[op_id].seg_id; //segment_identifier const& other_id = it->operations[other_op_id].seg_id; const bool first_in_range = m_seg_watcher.update(seg_id); if ( m_previous_turn_ptr ) { if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ ) { // real exit point - may be multiple if ( first_in_range || ! turn_on_the_same_ip(*m_previous_turn_ptr, *it) ) { update_exit(result); m_exit_detected = false; } // fake exit point, reset state else if ( op != overlay::operation_union ) { m_exit_detected = false; } } /*else*/ if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ ) { // real entry point if ( first_in_range || ! turn_on_the_same_ip(*m_previous_turn_ptr, *it) ) { update_enter(result); m_enter_detected = false; } // fake entry point, reset state else if ( op != overlay::operation_intersection ) { m_enter_detected = false; } } } if ( op == overlay::operation_union ) { // already set in interrupt policy //update(m_result); // ignore u/u //if ( it->operations[other_op_id].operation != overlay::operation_union ) { m_exit_detected = true; } } else if ( op == overlay::operation_intersection ) { // ignore i/i if ( it->operations[other_op_id].operation != overlay::operation_intersection ) { // already set in interrupt policy //update(result); //update(result); m_enter_detected = true; } } else if ( op == overlay::operation_blocked ) { // already set in interrupt policy } else // if ( op == overlay::operation_continue ) { // already set in interrupt policy } // store ref to previously analysed (valid) turn m_previous_turn_ptr = boost::addressof(*it); // and previously analysed (valid) operation m_previous_operation = op; } // it == last template void apply(Result & result) { //BOOST_ASSERT( first != last ); if ( m_exit_detected /*m_previous_operation == overlay::operation_union*/ ) { update_exit(result); m_exit_detected = false; } if ( m_enter_detected /*m_previous_operation == overlay::operation_intersection*/ ) { update_enter(result); m_enter_detected = false; } } template static inline void update_exit(Result & result) { update(result); update(result); } template static inline void update_enter(Result & result) { update(result); update(result); } private: segment_watcher m_seg_watcher; TurnInfo * m_previous_turn_ptr; overlay::operation_type m_previous_operation; bool m_enter_detected; bool m_exit_detected; }; // call analyser.apply() for each turn in range // IMPORTANT! The analyser is also called for the end iterator - last template static inline void analyse_each_turn(Result & res, Analyser & analyser, TurnIt first, TurnIt last) { if ( first == last ) return; for ( TurnIt it = first ; it != last ; ++it ) { analyser.apply(res, it); if ( res.interrupt ) return; } analyser.apply(res); } template class uncertain_rings_analyser { static const bool transpose_result = OpId != 0; static const int other_id = (OpId + 1) % 2; public: inline uncertain_rings_analyser(Result & result, Geometry const& geom, OtherGeometry const& other_geom) : geometry(geom), other_geometry(other_geom) , interrupt(result.interrupt) // just in case, could be false as well , m_result(result) , m_flags(0) { // check which relations must be analysed if ( ! may_update(m_result) && ! may_update(m_result) ) { m_flags |= 1; } if ( ! may_update(m_result) && ! may_update(m_result) ) { m_flags |= 2; } if ( ! may_update(m_result) && ! may_update(m_result) ) { m_flags |= 4; } } inline void no_turns(segment_identifier const& seg_id) { // if those flags are set nothing will change if ( (m_flags & 3) == 3 ) { return; } typename detail::sub_range_return_type::type range_ref = detail::sub_range(geometry, seg_id); if ( boost::empty(range_ref) ) { // TODO: throw an exception? return; // ignore } // TODO: possible optimization // if the range is an interior ring we may use other IPs generated for this single geometry // to know which other single geometries should be checked // TODO: optimize! e.g. use spatial index // O(N) - running it in a loop would gives O(NM) int const pig = detail::within::point_in_geometry(range::front(range_ref), other_geometry); //BOOST_ASSERT(pig != 0); if ( pig > 0 ) { update(m_result); update(m_result); m_flags |= 1; } else { update(m_result); update(m_result); m_flags |= 2; } // TODO: break if all things are set // also some of them could be checked outside, before the analysis // In this case we shouldn't relay just on dummy flags // Flags should be initialized with proper values // or the result should be checked directly // THIS IS ALSO TRUE FOR OTHER ANALYSERS! in L/L and L/A interrupt = m_flags == 7 || m_result.interrupt; } template inline void turns(TurnIt first, TurnIt last) { // if those flags are set nothing will change if ( (m_flags & 6) == 6 ) { return; } bool found_ii = false; bool found_uu = false; for ( TurnIt it = first ; it != last ; ++it ) { if ( it->operations[0].operation == overlay::operation_intersection && it->operations[1].operation == overlay::operation_intersection ) { // ignore exterior ring if ( it->operations[OpId].seg_id.ring_index >= 0 ) { found_ii = true; } } else if ( it->operations[0].operation == overlay::operation_union && it->operations[1].operation == overlay::operation_union ) { // ignore if u/u is for holes //if ( it->operations[OpId].seg_id.ring_index >= 0 // && it->operations[other_id].seg_id.ring_index >= 0 ) { found_uu = true; } } else // ignore { return; // don't interrupt } } // only i/i was generated for this ring if ( found_ii ) { //update(m_result); //update(m_result); update(m_result); update(m_result); m_flags |= 4; } // only u/u was generated for this ring if ( found_uu ) { update(m_result); update(m_result); m_flags |= 2; // not necessary since this will be checked in the next iteration // but increases the pruning strength // WARNING: this is not reflected in flags update(m_result); update(m_result); } interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future } Geometry const& geometry; OtherGeometry const& other_geometry; bool interrupt; private: Result & m_result; int m_flags; }; template class analyse_uncertain_rings { public: template static inline void apply(Analyser & analyser, TurnIt first, TurnIt last) { if ( first == last ) return; for_preceding_rings(analyser, *first); //analyser.per_turn(*first); TurnIt prev = first; for ( ++first ; first != last ; ++first, ++prev ) { // same multi if ( prev->operations[OpId].seg_id.multi_index == first->operations[OpId].seg_id.multi_index ) { // same ring if ( prev->operations[OpId].seg_id.ring_index == first->operations[OpId].seg_id.ring_index ) { //analyser.per_turn(*first); } // same multi, next ring else { //analyser.end_ring(*prev); analyser.turns(prev, first); //if ( prev->operations[OpId].seg_id.ring_index + 1 // < first->operations[OpId].seg_id.ring_index) { for_no_turns_rings(analyser, *first, prev->operations[OpId].seg_id.ring_index + 1, first->operations[OpId].seg_id.ring_index); } //analyser.per_turn(*first); } } // next multi else { //analyser.end_ring(*prev); analyser.turns(prev, first); for_following_rings(analyser, *prev); for_preceding_rings(analyser, *first); //analyser.per_turn(*first); } if ( analyser.interrupt ) { return; } } //analyser.end_ring(*prev); analyser.turns(prev, first); // first == last for_following_rings(analyser, *prev); } private: template static inline void for_preceding_rings(Analyser & analyser, Turn const& turn) { segment_identifier const& seg_id = turn.operations[OpId].seg_id; for_no_turns_rings(analyser, turn, -1, seg_id.ring_index); } template static inline void for_following_rings(Analyser & analyser, Turn const& turn) { segment_identifier const& seg_id = turn.operations[OpId].seg_id; int count = boost::numeric_cast( geometry::num_interior_rings( detail::single_geometry(analyser.geometry, seg_id))); for_no_turns_rings(analyser, turn, seg_id.ring_index + 1, count); } template static inline void for_no_turns_rings(Analyser & analyser, Turn const& turn, int first, int last) { segment_identifier seg_id = turn.operations[OpId].seg_id; for ( seg_id.ring_index = first ; seg_id.ring_index < last ; ++seg_id.ring_index ) { analyser.no_turns(seg_id); } } }; }; }} // namespace detail::relate #endif // DOXYGEN_NO_DETAIL }} // namespace boost::geometry #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_AREAL_AREAL_HPP