diff options
Diffstat (limited to 'boost/geometry/algorithms/detail/relate')
7 files changed, 486 insertions, 123 deletions
diff --git a/boost/geometry/algorithms/detail/relate/areal_areal.hpp b/boost/geometry/algorithms/detail/relate/areal_areal.hpp index 31d206ac99..2859841de4 100644 --- a/boost/geometry/algorithms/detail/relate/areal_areal.hpp +++ b/boost/geometry/algorithms/detail/relate/areal_areal.hpp @@ -2,19 +2,21 @@ // 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. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // 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 <boost/geometry/core/topological_dimension.hpp> + +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/range.hpp> #include <boost/geometry/algorithms/num_interior_rings.hpp> @@ -197,7 +199,7 @@ struct areal_areal // The result should be FFFFFFFFF relate::set<exterior, exterior, result_dimension<Geometry2>::value>(result);// FFFFFFFFd, d in [1,9] or T - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; // get and analyse turns @@ -207,17 +209,17 @@ struct areal_areal interrupt_policy_areal_areal<Result> interrupt_policy(geometry1, geometry2, result); turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; no_turns_aa_pred<Geometry2, Result, false> pred1(geometry2, result); for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; no_turns_aa_pred<Geometry1, Result, true> pred2(geometry1, result); for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; if ( turns.empty() ) @@ -242,7 +244,7 @@ struct areal_areal turns_analyser<turn_type, 0> analyser; analyse_each_turn(result, analyser, turns.begin(), turns.end()); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; } @@ -257,7 +259,7 @@ struct areal_areal 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 ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; } } @@ -281,7 +283,7 @@ struct areal_areal turns_analyser<turn_type, 1> analyser; analyse_each_turn(result, analyser, turns.begin(), turns.end()); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; } @@ -549,7 +551,7 @@ struct areal_areal { analyser.apply(res, it); - if ( res.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(res.interrupt) ) return; } @@ -573,8 +575,7 @@ struct areal_areal { // check which relations must be analysed - if ( ! may_update<interior, interior, '2', transpose_result>(m_result) - && ! may_update<boundary, interior, '1', transpose_result>(m_result) ) + if ( ! may_update<interior, interior, '2', transpose_result>(m_result) ) { m_flags |= 1; } @@ -595,7 +596,7 @@ struct areal_areal inline void no_turns(segment_identifier const& seg_id) { // if those flags are set nothing will change - if ( (m_flags & 3) == 3 ) + if ( m_flags == 7 ) { return; } @@ -614,15 +615,18 @@ struct areal_areal // 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) + // O(N) - running it in a loop 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<boundary, interior, '1', transpose_result>(m_result); update<interior, interior, '2', transpose_result>(m_result); m_flags |= 1; + + update<boundary, interior, '1', transpose_result>(m_result); + update<exterior, interior, '2', transpose_result>(m_result); + m_flags |= 4; } else { @@ -696,12 +700,6 @@ struct areal_areal update<boundary, exterior, '1', transpose_result>(m_result); update<interior, exterior, '2', transpose_result>(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<exterior, boundary, '1', transpose_result>(m_result); - update<exterior, interior, '2', transpose_result>(m_result); } interrupt = m_flags == 7 || m_result.interrupt; // interrupt if the result won't be changed in the future @@ -795,7 +793,8 @@ struct areal_areal { segment_identifier const& seg_id = turn.operations[OpId].seg_id; - int count = boost::numeric_cast<int>( + signed_index_type + count = boost::numeric_cast<signed_index_type>( geometry::num_interior_rings( detail::single_geometry(analyser.geometry, seg_id))); @@ -803,7 +802,10 @@ struct areal_areal } template <typename Analyser, typename Turn> - static inline void for_no_turns_rings(Analyser & analyser, Turn const& turn, int first, int last) + static inline void for_no_turns_rings(Analyser & analyser, + Turn const& turn, + signed_index_type first, + signed_index_type last) { segment_identifier seg_id = turn.operations[OpId].seg_id; diff --git a/boost/geometry/algorithms/detail/relate/follow_helpers.hpp b/boost/geometry/algorithms/detail/relate/follow_helpers.hpp index 78fa03798d..2c44b009e7 100644 --- a/boost/geometry/algorithms/detail/relate/follow_helpers.hpp +++ b/boost/geometry/algorithms/detail/relate/follow_helpers.hpp @@ -14,6 +14,7 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_FOLLOW_HELPERS_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_FOLLOW_HELPERS_HPP +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/range.hpp> //#include <boost/geometry/algorithms/detail/sub_range.hpp> @@ -98,9 +99,9 @@ struct for_each_disjoint_geometry_if<OpId, Geometry, Tag, true> std::vector<bool> detected_intersections(count, false); for ( TurnIt it = first ; it != last ; ++it ) { - int multi_index = it->operations[OpId].seg_id.multi_index; + signed_index_type multi_index = it->operations[OpId].seg_id.multi_index; BOOST_ASSERT(multi_index >= 0); - std::size_t index = static_cast<std::size_t>(multi_index); + std::size_t const index = static_cast<std::size_t>(multi_index); BOOST_ASSERT(index < count); detected_intersections[index] = true; } @@ -116,8 +117,8 @@ struct for_each_disjoint_geometry_if<OpId, Geometry, Tag, true> if ( *it == false ) { found = true; - bool cont = pred(range::at(geometry, - std::distance(detected_intersections.begin(), it))); + std::size_t const index = std::size_t(std::distance(detected_intersections.begin(), it)); + bool cont = pred(range::at(geometry, index)); if ( !cont ) break; } @@ -375,14 +376,14 @@ static inline bool is_ip_on_boundary(IntersectionPoint const& ip, bool res = false; // IP on the last point of the linestring - if ( (BoundaryQuery == boundary_back || BoundaryQuery == boundary_any) + if ( BOOST_GEOMETRY_CONDITION(BoundaryQuery == boundary_back || BoundaryQuery == boundary_any) && operation_info.position == overlay::position_back ) { // check if this point is a boundary res = boundary_checker.template is_endpoint_boundary<boundary_back>(ip); } // IP on the last point of the linestring - else if ( (BoundaryQuery == boundary_front || BoundaryQuery == boundary_any) + else if ( BOOST_GEOMETRY_CONDITION(BoundaryQuery == boundary_front || BoundaryQuery == boundary_any) && operation_info.position == overlay::position_front ) { // check if this point is a boundary diff --git a/boost/geometry/algorithms/detail/relate/linear_areal.hpp b/boost/geometry/algorithms/detail/relate/linear_areal.hpp index 3159ab329d..7d85a1d9a1 100644 --- a/boost/geometry/algorithms/detail/relate/linear_areal.hpp +++ b/boost/geometry/algorithms/detail/relate/linear_areal.hpp @@ -2,21 +2,23 @@ // 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. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // 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_LINEAR_AREAL_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LINEAR_AREAL_HPP #include <boost/core/ignore_unused.hpp> #include <boost/geometry/core/topological_dimension.hpp> + +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/range.hpp> #include <boost/geometry/algorithms/num_interior_rings.hpp> @@ -193,6 +195,33 @@ struct linear_areal typedef typename geometry::point_type<Geometry1>::type point1_type; typedef typename geometry::point_type<Geometry2>::type point2_type; + + template <typename Geometry> + struct is_multi + : boost::is_base_of + < + multi_tag, + typename tag<Geometry>::type + > + {}; + + template <typename Geom1, typename Geom2> + struct multi_turn_info + : turns::get_turns<Geom1, Geom2>::turn_info + { + multi_turn_info() : priority(0) {} + int priority; // single-geometry sorting priority + }; + + template <typename Geom1, typename Geom2> + struct turn_info_type + : boost::mpl::if_c + < + is_multi<Geometry2>::value, + multi_turn_info<Geom1, Geom2>, + typename turns::get_turns<Geom1, Geom2>::turn_info + > + {}; template <typename Result> static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Result & result) @@ -202,17 +231,17 @@ struct linear_areal // The result should be FFFFFFFFF relate::set<exterior, exterior, result_dimension<Geometry2>::value, TransposeResult>(result);// FFFFFFFFd, d in [1,9] or T - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; // get and analyse turns - typedef typename turns::get_turns<Geometry1, Geometry2>::turn_info turn_type; + typedef typename turn_info_type<Geometry1, Geometry2>::type turn_type; std::vector<turn_type> turns; interrupt_policy_linear_areal<Geometry2, Result> interrupt_policy(geometry2, result); turns::get_turns<Geometry1, Geometry2>::apply(turns, geometry1, geometry2, interrupt_policy); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; boundary_checker<Geometry1> boundary_checker1(geometry1); @@ -224,12 +253,12 @@ struct linear_areal TransposeResult > pred1(geometry2, result, boundary_checker1); for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; no_turns_la_areal_pred<Result, !TransposeResult> pred2(result); for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; if ( turns.empty() ) @@ -238,14 +267,11 @@ struct linear_areal // This is set here because in the case if empty Areal geometry were passed // those shouldn't be set relate::set<exterior, interior, '2', TransposeResult>(result);// FFFFFF2Fd - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; { - // for different multi or same ring id: x, u, i, c - // for same multi and different ring id: c, i, u, x - typedef turns::less<0, turns::less_op_linear_areal<0> > less; - std::sort(turns.begin(), turns.end(), less()); + sort_dispatch(turns.begin(), turns.end(), is_multi<Geometry2>()); turns_analyser<turn_type> analyser; analyse_each_turn(result, analyser, @@ -253,7 +279,7 @@ struct linear_areal geometry1, geometry2, boundary_checker1); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; } @@ -297,7 +323,7 @@ struct linear_areal // if there was some previous ring if ( prev_seg_id_ptr != NULL ) { - int const next_ring_index = prev_seg_id_ptr->ring_index + 1; + signed_index_type const next_ring_index = prev_seg_id_ptr->ring_index + 1; BOOST_ASSERT(next_ring_index >= 0); // if one of the last rings of previous single geometry was ommited @@ -368,7 +394,7 @@ struct linear_areal // if there was some previous ring if ( prev_seg_id_ptr != NULL ) { - int const next_ring_index = prev_seg_id_ptr->ring_index + 1; + signed_index_type const next_ring_index = prev_seg_id_ptr->ring_index + 1; BOOST_ASSERT(next_ring_index >= 0); // if one of the last rings of previous single geometry was ommited @@ -383,6 +409,127 @@ struct linear_areal } } + template <typename It, typename Pred, typename Comp> + static void for_each_equal_range(It first, It last, Pred pred, Comp comp) + { + if ( first == last ) + return; + + It first_equal = first; + It prev = first; + for ( ++first ; ; ++first, ++prev ) + { + if ( first == last || !comp(*prev, *first) ) + { + pred(first_equal, first); + first_equal = first; + } + + if ( first == last ) + break; + } + } + + struct same_ip + { + template <typename Turn> + bool operator()(Turn const& left, Turn const& right) const + { + return left.operations[0].seg_id == right.operations[0].seg_id + && left.operations[0].fraction == right.operations[0].fraction; + } + }; + + struct same_ip_and_multi_index + { + template <typename Turn> + bool operator()(Turn const& left, Turn const& right) const + { + return same_ip()(left, right) + && left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index; + } + }; + + template <typename OpToPriority> + struct set_turns_group_priority + { + template <typename TurnIt> + void operator()(TurnIt first, TurnIt last) const + { + BOOST_ASSERT(first != last); + static OpToPriority op_to_priority; + // find the operation with the least priority + int least_priority = op_to_priority(first->operations[0]); + for ( TurnIt it = first + 1 ; it != last ; ++it ) + { + int priority = op_to_priority(it->operations[0]); + if ( priority < least_priority ) + least_priority = priority; + } + // set the least priority for all turns of the group + for ( TurnIt it = first ; it != last ; ++it ) + { + it->priority = least_priority; + } + } + }; + + template <typename SingleLess> + struct sort_turns_group + { + struct less + { + template <typename Turn> + bool operator()(Turn const& left, Turn const& right) const + { + return left.operations[1].seg_id.multi_index == right.operations[1].seg_id.multi_index ? + SingleLess()(left, right) : + left.priority < right.priority; + } + }; + + template <typename TurnIt> + void operator()(TurnIt first, TurnIt last) const + { + std::sort(first, last, less()); + } + }; + + template <typename TurnIt> + static void sort_dispatch(TurnIt first, TurnIt last, boost::true_type const& /*is_multi*/) + { + // sort turns by Linear seg_id, then by fraction, then by other multi_index + typedef turns::less<0, turns::less_other_multi_index<0> > less; + std::sort(first, last, less()); + + // For the same IP and multi_index - the same other's single geometry + // set priorities as the least operation found for the whole single geometry + // so e.g. single geometries containing 'u' will always be before those only containing 'i' + typedef turns::op_to_int<0,2,3,1,4,0> op_to_int_xuic; + for_each_equal_range(first, last, + set_turns_group_priority<op_to_int_xuic>(), // least operation in xuic order + same_ip_and_multi_index()); // other's multi_index + + // When priorities for single geometries are set now sort turns for the same IP + // if multi_index is the same sort them according to the single-less + // else use priority of the whole single-geometry set earlier + typedef turns::less<0, turns::less_op_linear_areal_single<0> > single_less; + for_each_equal_range(first, last, + sort_turns_group<single_less>(), + same_ip()); + } + + template <typename TurnIt> + static void sort_dispatch(TurnIt first, TurnIt last, boost::false_type const& /*is_multi*/) + { + // sort turns by Linear seg_id, then by fraction, then + // for same ring id: x, u, i, c + // for different ring id: c, i, u, x + typedef turns::less<0, turns::less_op_linear_areal_single<0> > less; + std::sort(first, last, less()); + } + + // interrupt policy which may be passed to get_turns to interrupt the analysis // based on the info in the passed result/mask template <typename Areal, typename Result> @@ -458,6 +605,8 @@ struct linear_areal , m_boundary_counter(0) , m_interior_detected(false) , m_first_interior_other_id_ptr(NULL) + , m_first_from_unknown(false) + , m_first_from_unknown_boundary_detected(false) {} template <typename Result, @@ -485,6 +634,11 @@ struct linear_areal const bool first_in_range = m_seg_watcher.update(seg_id); + // TODO: should apply() for the post-last ip be called if first_in_range ? + // this would unify how last points in ranges are handled + // possibly replacing parts of the code below + // e.g. for is_multi and m_interior_detected + // handle possible exit bool fake_enter_detected = false; if ( m_exit_watcher.get_exit_operation() == overlay::operation_union ) @@ -495,8 +649,24 @@ struct linear_areal { m_exit_watcher.reset_detected_exit(); - // not the last IP update<interior, exterior, '1', TransposeResult>(res); + + // next single geometry + if ( first_in_range && m_previous_turn_ptr ) + { + // NOTE: similar code is in the post-last-ip-apply() + segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id; + + bool const prev_back_b = is_endpoint_on_boundary<boundary_back>( + range::back(sub_range(geometry, prev_seg_id)), + boundary_checker); + + // if there is a boundary on the last point + if ( prev_back_b ) + { + update<boundary, exterior, '0', TransposeResult>(res); + } + } } // fake exit point, reset state else if ( op == overlay::operation_intersection @@ -508,9 +678,12 @@ struct linear_areal } else if ( m_exit_watcher.get_exit_operation() == overlay::operation_blocked ) { - // ignore multiple BLOCKs - if ( op == overlay::operation_blocked ) + // ignore multiple BLOCKs for this same single geometry1 + if ( op == overlay::operation_blocked + && seg_id.multi_index == m_previous_turn_ptr->operations[op_id].seg_id.multi_index ) + { return; + } if ( ( op == overlay::operation_intersection || op == overlay::operation_continue ) @@ -522,11 +695,39 @@ struct linear_areal m_exit_watcher.reset_detected_exit(); } + if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) + && m_first_from_unknown ) + { + // For MultiPolygon many x/u operations may be generated as a first IP + // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString + // then we know that the LineString is outside + // Similar with the u/u turns, if it was the first one it doesn't mean that the + // Linestring came from the exterior + if ( ( m_previous_operation == overlay::operation_blocked + && ( op != overlay::operation_blocked // operation different than block + || seg_id.multi_index != m_previous_turn_ptr->operations[op_id].seg_id.multi_index ) ) // or the next single-geometry + || ( m_previous_operation == overlay::operation_union + && ! turn_on_the_same_ip<op_id>(*m_previous_turn_ptr, *it) ) + ) + { + update<interior, exterior, '1', TransposeResult>(res); + if ( m_first_from_unknown_boundary_detected ) + { + update<boundary, exterior, '0', TransposeResult>(res); + } + + m_first_from_unknown = false; + m_first_from_unknown_boundary_detected = false; + } + } + // NOTE: THE WHOLE m_interior_detected HANDLING IS HERE BECAUSE WE CAN'T EFFICIENTLY SORT TURNS (CORRECTLY) // BECAUSE THE SAME IP MAY BE REPRESENTED BY TWO SEGMENTS WITH DIFFERENT DISTANCES // IT WOULD REQUIRE THE CALCULATION OF MAX DISTANCE // TODO: WE COULD GET RID OF THE TEST IF THE DISTANCES WERE NORMALIZED +// UPDATE: THEY SHOULD BE NORMALIZED NOW + // TODO: THIS IS POTENTIALLY ERROREOUS! // THIS ALGORITHM DEPENDS ON SOME SPECIFIC SEQUENCE OF OPERATIONS // IT WOULD GIVE WRONG RESULTS E.G. @@ -540,6 +741,30 @@ struct linear_areal { update<interior, interior, '1', TransposeResult>(res); m_interior_detected = false; + + // new range detected - reset previous state and check the boundary + if ( first_in_range ) + { + // actually it should be != NULL if m_interior_detected + // so an assert could be checked here + if ( m_previous_turn_ptr ) + { + segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id; + + bool const prev_back_b = is_endpoint_on_boundary<boundary_back>( + range::back(sub_range(geometry, prev_seg_id)), + boundary_checker); + + // if there is a boundary on the last point + if ( prev_back_b ) + { + update<boundary, interior, '0', TransposeResult>(res); + } + } + + // The exit_watcher is reset below + // m_exit_watcher.reset(); + } } // fake interior overlap else if ( op == overlay::operation_continue ) @@ -560,10 +785,20 @@ struct linear_areal } } + // NOTE: If post-last-ip apply() was called this wouldn't be needed + if ( first_in_range ) + { + m_exit_watcher.reset(); + m_boundary_counter = 0; + m_first_from_unknown = false; + m_first_from_unknown_boundary_detected = false; + } + // i/u, c/u if ( op == overlay::operation_intersection || op == overlay::operation_continue ) // operation_boundary/operation_boundary_intersection { + bool const first_point = first_in_range || m_first_from_unknown; bool no_enters_detected = m_exit_watcher.is_outside(); m_exit_watcher.enter(*it); @@ -592,7 +827,7 @@ struct linear_areal { // don't add to the count for all met boundaries // only if this is the "new" boundary - if ( first_in_range || !it->operations[op_id].is_collinear ) + if ( first_point || !it->operations[op_id].is_collinear ) ++m_boundary_counter; update<interior, boundary, '1', TransposeResult>(res); @@ -619,7 +854,7 @@ struct linear_areal && it->operations[op_id].position != overlay::position_front ) { // TODO: calculate_from_inside() is only needed if the current Linestring is not closed - bool const from_inside = first_in_range + bool const from_inside = first_point && calculate_from_inside(geometry, other_geometry, *it); @@ -630,7 +865,7 @@ struct linear_areal update<interior, exterior, '1', TransposeResult>(res); // if it's the first IP then the first point is outside - if ( first_in_range ) + if ( first_point ) { bool const front_b = is_endpoint_on_boundary<boundary_front>( range::front(sub_range(geometry, seg_id)), @@ -647,6 +882,12 @@ struct linear_areal } } } + + if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) ) + { + m_first_from_unknown = false; + m_first_from_unknown_boundary_detected = false; + } } // u/u, x/u else if ( op == overlay::operation_union || op == overlay::operation_blocked ) @@ -709,7 +950,11 @@ struct linear_areal if ( it->operations[op_id].position != overlay::position_front ) { // TODO: calculate_from_inside() is only needed if the current Linestring is not closed - bool const first_from_inside = first_in_range + // NOTE: this is not enough for MultiPolygon and operation_blocked + // For LS/MultiPolygon multiple x/u turns may be generated + // the first checked Polygon may be the one which LS is outside for. + bool const first_point = first_in_range || m_first_from_unknown; + bool const first_from_inside = first_point && calculate_from_inside(geometry, other_geometry, *it); @@ -719,14 +964,26 @@ struct linear_areal // notify the exit_watcher that we started inside m_exit_watcher.enter(*it); + // and reset unknown flags since we know that we started inside + m_first_from_unknown = false; + m_first_from_unknown_boundary_detected = false; } else { - update<interior, exterior, '1', TransposeResult>(res); + if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) + /*&& ( op == overlay::operation_blocked + || op == overlay::operation_union )*/ ) // if we're here it's u or x + { + m_first_from_unknown = true; + } + else + { + update<interior, exterior, '1', TransposeResult>(res); + } } // first IP on the last segment point - this means that the first point is outside or inside - if ( first_in_range && ( !this_b || op_blocked ) ) + if ( first_point && ( !this_b || op_blocked ) ) { bool const front_b = is_endpoint_on_boundary<boundary_front>( range::front(sub_range(geometry, seg_id)), @@ -736,9 +993,23 @@ struct linear_areal if ( front_b ) { if ( first_from_inside ) + { update<boundary, interior, '0', TransposeResult>(res); + } else - update<boundary, exterior, '0', TransposeResult>(res); + { + if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) + /*&& ( op == overlay::operation_blocked + || op == overlay::operation_union )*/ ) // if we're here it's u or x + { + BOOST_ASSERT(m_first_from_unknown); + m_first_from_unknown_boundary_detected = true; + } + else + { + update<boundary, exterior, '0', TransposeResult>(res); + } + } } } } @@ -774,6 +1045,23 @@ struct linear_areal boost::ignore_unused(first, last); //BOOST_ASSERT( first != last ); + // For MultiPolygon many x/u operations may be generated as a first IP + // if for all turns x/u was generated and any of the Polygons doesn't contain the LineString + // then we know that the LineString is outside + if ( BOOST_GEOMETRY_CONDITION( is_multi<OtherGeometry>::value ) + && m_first_from_unknown ) + { + update<interior, exterior, '1', TransposeResult>(res); + if ( m_first_from_unknown_boundary_detected ) + { + update<boundary, exterior, '0', TransposeResult>(res); + } + + // done below + //m_first_from_unknown = false; + //m_first_from_unknown_boundary_detected = false; + } + // here, the possible exit is the real one // we know that we entered and now we exit if ( /*m_exit_watcher.get_exit_operation() == overlay::operation_union // THIS CHECK IS REDUNDANT @@ -822,12 +1110,37 @@ struct linear_areal } } - BOOST_ASSERT_MSG(m_previous_operation != overlay::operation_continue, - "Unexpected operation! Probably the error in get_turns(L,A) or relate(L,A)"); + // This condition may be false if the Linestring is lying on the Polygon's collinear spike + // if Polygon's spikes are not handled in get_turns() or relate() (they currently aren't) + //BOOST_ASSERT_MSG(m_previous_operation != overlay::operation_continue, + // "Unexpected operation! Probably the error in get_turns(L,A) or relate(L,A)"); + // Currently one c/c turn is generated for the exit + // when a Linestring is going out on a collinear spike + // When a Linestring is going in on a collinear spike + // the turn is not generated for the entry + // So assume it's the former + if ( m_previous_operation == overlay::operation_continue ) + { + update<interior, exterior, '1', TransposeResult>(res); + + segment_identifier const& prev_seg_id = m_previous_turn_ptr->operations[op_id].seg_id; + + bool const prev_back_b = is_endpoint_on_boundary<boundary_back>( + range::back(sub_range(geometry, prev_seg_id)), + boundary_checker); + + // if there is a boundary on the last point + if ( prev_back_b ) + { + update<boundary, exterior, '0', TransposeResult>(res); + } + } // Reset exit watcher before the analysis of the next Linestring m_exit_watcher.reset(); m_boundary_counter = 0; + m_first_from_unknown = false; + m_first_from_unknown_boundary_detected = false; } // check if the passed turn's segment of Linear geometry arrived @@ -852,8 +1165,8 @@ struct linear_areal BOOST_ASSERT(s2 > 2); std::size_t const seg_count2 = s2 - 1; - std::size_t const p_seg_ij = turn.operations[op_id].seg_id.segment_index; - std::size_t const q_seg_ij = turn.operations[other_op_id].seg_id.segment_index; + std::size_t const p_seg_ij = static_cast<std::size_t>(turn.operations[op_id].seg_id.segment_index); + std::size_t const q_seg_ij = static_cast<std::size_t>(turn.operations[other_op_id].seg_id.segment_index); BOOST_ASSERT(p_seg_ij + 1 < boost::size(range1)); BOOST_ASSERT(q_seg_ij + 1 < s2); @@ -876,7 +1189,7 @@ struct linear_areal // TODO: the following function should return immediately, however the worst case complexity is O(N) // It would be good to replace it with some O(1) mechanism range2_iterator qk_it = find_next_non_duplicated(boost::begin(range2), - boost::begin(range2) + q_seg_jk, + range::pos(range2, q_seg_jk), boost::end(range2)); // Will this sequence of points be always correct? @@ -945,6 +1258,8 @@ struct linear_areal unsigned m_boundary_counter; bool m_interior_detected; const segment_identifier * m_first_interior_other_id_ptr; + bool m_first_from_unknown; + bool m_first_from_unknown_boundary_detected; }; // call analyser.apply() for each turn in range @@ -971,7 +1286,7 @@ struct linear_areal geometry, other_geometry, boundary_checker); - if ( res.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( res.interrupt ) ) return; } @@ -1017,8 +1332,8 @@ struct linear_areal if ( first == last ) return last; - int const multi_index = first->operations[1].seg_id.multi_index; - int const ring_index = first->operations[1].seg_id.ring_index; + signed_index_type const multi_index = first->operations[1].seg_id.multi_index; + signed_index_type const ring_index = first->operations[1].seg_id.ring_index; fun(*first); ++first; diff --git a/boost/geometry/algorithms/detail/relate/linear_linear.hpp b/boost/geometry/algorithms/detail/relate/linear_linear.hpp index 263c82de56..20a22c3018 100644 --- a/boost/geometry/algorithms/detail/relate/linear_linear.hpp +++ b/boost/geometry/algorithms/detail/relate/linear_linear.hpp @@ -2,8 +2,8 @@ // 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. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 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 @@ -16,6 +16,7 @@ #include <boost/core/ignore_unused.hpp> +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/range.hpp> #include <boost/geometry/algorithms/detail/sub_range.hpp> @@ -117,7 +118,7 @@ struct linear_linear { // The result should be FFFFFFFFF relate::set<exterior, exterior, result_dimension<Geometry1>::value>(result);// FFFFFFFFd, d in [1,9] or T - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; // get and analyse turns @@ -133,19 +134,19 @@ struct linear_linear detail::get_turns::get_turn_info_type<Geometry1, Geometry2, turns::assign_policy<true> > >::apply(turns, geometry1, geometry2, interrupt_policy); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; boundary_checker<Geometry1> boundary_checker1(geometry1); disjoint_linestring_pred<Result, boundary_checker<Geometry1>, false> pred1(result, boundary_checker1); for_each_disjoint_geometry_if<0, Geometry1>::apply(turns.begin(), turns.end(), geometry1, pred1); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; boundary_checker<Geometry2> boundary_checker2(geometry2); disjoint_linestring_pred<Result, boundary_checker<Geometry2>, true> pred2(result, boundary_checker2); for_each_disjoint_geometry_if<1, Geometry2>::apply(turns.begin(), turns.end(), geometry2, pred2); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; if ( turns.empty() ) @@ -171,7 +172,7 @@ struct linear_linear boundary_checker1, boundary_checker2); } - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( result.interrupt ) ) return; if ( may_update<interior, interior, '1', true>(result) @@ -250,6 +251,7 @@ struct linear_linear : m_previous_turn_ptr(NULL) , m_previous_operation(overlay::operation_none) , m_degenerated_turn_ptr(NULL) + , m_collinear_spike_exit(false) {} template <typename Result, @@ -383,7 +385,8 @@ struct linear_linear // if we didn't enter in the past, we were outside if ( was_outside && ! fake_enter_detected - && it->operations[op_id].position != overlay::position_front ) + && it->operations[op_id].position != overlay::position_front + && ! m_collinear_spike_exit ) { update<interior, exterior, '1', transpose_result>(res); @@ -402,6 +405,8 @@ struct linear_linear } } } + + m_collinear_spike_exit = false; } // u/i, u/u, u/x, x/i, x/u, x/x else if ( op == overlay::operation_union || op == overlay::operation_blocked ) @@ -418,6 +423,11 @@ struct linear_linear if ( !was_outside && is_collinear ) { m_exit_watcher.exit(*it, false); + // if the position is not set to back it must be a spike + if ( it->operations[op_id].position != overlay::position_back ) + { + m_collinear_spike_exit = true; + } } bool const op_blocked = op == overlay::operation_blocked; @@ -456,6 +466,7 @@ struct linear_linear // if we are truly outside if ( was_outside && it->operations[op_id].position != overlay::position_front + && ! m_collinear_spike_exit /*&& !is_collinear*/ ) { update<interior, exterior, '1', transpose_result>(res); @@ -526,6 +537,7 @@ struct linear_linear && ( !this_b || op_blocked ) && was_outside && it->operations[op_id].position != overlay::position_front + && ! m_collinear_spike_exit /*&& !is_collinear*/ ) { bool const front_b = is_endpoint_on_boundary<boundary_front>( @@ -607,6 +619,10 @@ struct linear_linear m_previous_turn_ptr = NULL; m_previous_operation = overlay::operation_none; m_degenerated_turn_ptr = NULL; + + // actually if this is set to true here there is some error + // in get_turns_ll or relate_ll, an assert could be checked here + m_collinear_spike_exit = false; } template <typename Result, @@ -724,6 +740,7 @@ struct linear_linear const TurnInfo * m_previous_turn_ptr; overlay::operation_type m_previous_operation; const TurnInfo * m_degenerated_turn_ptr; + bool m_collinear_spike_exit; }; template <typename Result, @@ -750,7 +767,7 @@ struct linear_linear geometry, other_geometry, boundary_checker, other_boundary_checker); - if ( res.interrupt ) + if ( BOOST_GEOMETRY_CONDITION( res.interrupt ) ) return; } diff --git a/boost/geometry/algorithms/detail/relate/point_geometry.hpp b/boost/geometry/algorithms/detail/relate/point_geometry.hpp index 86c236d357..62ab100919 100644 --- a/boost/geometry/algorithms/detail/relate/point_geometry.hpp +++ b/boost/geometry/algorithms/detail/relate/point_geometry.hpp @@ -2,15 +2,15 @@ // 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. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // 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_POINT_GEOMETRY_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_POINT_GEOMETRY_HPP @@ -20,6 +20,8 @@ #include <boost/geometry/algorithms/detail/relate/topology_check.hpp> +#include <boost/geometry/util/condition.hpp> + namespace boost { namespace geometry { @@ -54,7 +56,7 @@ struct point_geometry set<exterior, exterior, result_dimension<Point>::value, Transpose>(result); - if ( result.interrupt ) + if ( BOOST_GEOMETRY_CONDITION(result.interrupt) ) return; // the point is on the boundary diff --git a/boost/geometry/algorithms/detail/relate/result.hpp b/boost/geometry/algorithms/detail/relate/result.hpp index 1bcb5275d2..e26bda67f2 100644 --- a/boost/geometry/algorithms/detail/relate/result.hpp +++ b/boost/geometry/algorithms/detail/relate/result.hpp @@ -2,15 +2,15 @@ // 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. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // 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_RESULT_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_RESULT_HPP @@ -24,6 +24,7 @@ #include <boost/mpl/vector_c.hpp> #include <boost/geometry/core/topological_dimension.hpp> +#include <boost/geometry/util/condition.hpp> // TEMP - move this header to geometry/detail #include <boost/geometry/index/detail/tuples.hpp> @@ -235,17 +236,17 @@ struct interrupt_dispatch<Mask, true> static inline bool apply(Mask const& mask) { char m = mask.template get<F1, F2>(); - return check<V>(m); + return check_element<V>(m); } template <char V> - static inline bool check(char m) + static inline bool check_element(char m) { - if ( V >= '0' && V <= '9' ) + if ( BOOST_GEOMETRY_CONDITION(V >= '0' && V <= '9') ) { return m == 'F' || ( m < V && m >= '0' && m <= '9' ); } - else if ( V == 'T' ) + else if ( BOOST_GEOMETRY_CONDITION(V == 'T') ) { return m == 'F'; } @@ -394,7 +395,7 @@ inline bool may_update(Mask const& mask, Matrix const& matrix) ::template apply<F1, F2, D>(mask, matrix); } -// check() +// check_matrix() template <typename Mask> struct check_dispatch @@ -485,7 +486,7 @@ struct check_dispatch< boost::tuples::cons<Head, Tail> > }; template <typename Mask, typename Matrix> -inline bool check(Mask const& mask, Matrix const& matrix) +inline bool check_matrix(Mask const& mask, Matrix const& matrix) { return check_dispatch<Mask>::apply(mask, matrix); } @@ -546,7 +547,7 @@ public: result_type result() const { return !interrupt - && check(m_mask, static_cast<base_t const&>(*this)); + && check_matrix(m_mask, static_cast<base_t const&>(*this)); } template <field F1, field F2, char D> @@ -964,7 +965,7 @@ struct static_check_dispatch<StaticMask, true> }; template <typename StaticMask> -struct static_check +struct static_check_matrix { template <typename Matrix> static inline bool apply(Matrix const& matrix) @@ -997,7 +998,7 @@ public: result_type result() const { return (!Interrupt || !interrupt) - && static_check<StaticMask>:: + && static_check_matrix<StaticMask>:: apply(static_cast<base_t const&>(*this)); } diff --git a/boost/geometry/algorithms/detail/relate/turns.hpp b/boost/geometry/algorithms/detail/relate/turns.hpp index a2e56a8882..636c9756d8 100644 --- a/boost/geometry/algorithms/detail/relate/turns.hpp +++ b/boost/geometry/algorithms/detail/relate/turns.hpp @@ -1,16 +1,17 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2013, 2014. -// Modifications copyright (c) 2013-2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // 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_TURNS_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_TURNS_HPP @@ -19,8 +20,12 @@ #include <boost/geometry/algorithms/detail/overlay/get_turns.hpp> #include <boost/geometry/algorithms/detail/overlay/get_turn_info.hpp> +#include <boost/geometry/policies/robustness/get_rescale_policy.hpp> +#include <boost/geometry/policies/robustness/no_rescale_policy.hpp> + #include <boost/type_traits/is_base_of.hpp> + namespace boost { namespace geometry { #ifndef DOXYGEN_NO_DETAIL @@ -35,10 +40,16 @@ struct assign_policy // GET_TURNS -template <typename Geometry1, - typename Geometry2, - typename GetTurnPolicy - = detail::get_turns::get_turn_info_type<Geometry1, Geometry2, assign_policy<> > > +template +< + typename Geometry1, + typename Geometry2, + typename GetTurnPolicy = detail::get_turns::get_turn_info_type + < + Geometry1, Geometry2, assign_policy<> + >, + typename RobustPolicy = detail::no_rescale_policy +> struct get_turns { typedef typename geometry::point_type<Geometry1>::type point1_type; @@ -46,11 +57,14 @@ struct get_turns typedef overlay::turn_info < point1_type, - typename segment_ratio_type<point1_type, detail::no_rescale_policy>::type, + typename segment_ratio_type<point1_type, RobustPolicy>::type, typename detail::get_turns::turn_operation_type < Geometry1, Geometry2, - typename segment_ratio_type<point1_type, detail::no_rescale_policy>::type + typename segment_ratio_type + < + point1_type, RobustPolicy + >::type >::type > turn_info; @@ -73,6 +87,12 @@ struct get_turns static const bool reverse1 = detail::overlay::do_reverse<geometry::point_order<Geometry1>::value>::value; static const bool reverse2 = detail::overlay::do_reverse<geometry::point_order<Geometry2>::value>::value; + RobustPolicy robust_policy = geometry::get_rescale_policy + < + RobustPolicy + >(geometry1, geometry2); + + dispatch::get_turns < typename geometry::tag<Geometry1>::type, @@ -83,7 +103,7 @@ struct get_turns reverse2, GetTurnPolicy >::apply(0, geometry1, 1, geometry2, - detail::no_rescale_policy(), turns, interrupt_policy); + robust_policy, turns, interrupt_policy); } }; @@ -125,7 +145,7 @@ struct less_op_linear_linear {}; template <std::size_t OpId> -struct less_op_linear_areal +struct less_op_linear_areal_single { template <typename Turn> inline bool operator()(Turn const& left, Turn const& right) const @@ -137,27 +157,19 @@ struct less_op_linear_areal segment_identifier const& left_other_seg_id = left.operations[other_op_id].seg_id; segment_identifier const& right_other_seg_id = right.operations[other_op_id].seg_id; - if ( left_other_seg_id.multi_index == right_other_seg_id.multi_index ) - { - typedef typename Turn::turn_operation_type operation_type; - operation_type const& left_operation = left.operations[OpId]; - operation_type const& right_operation = right.operations[OpId]; + typedef typename Turn::turn_operation_type operation_type; + operation_type const& left_operation = left.operations[OpId]; + operation_type const& right_operation = right.operations[OpId]; - if ( left_other_seg_id.ring_index == right_other_seg_id.ring_index ) - { - return op_to_int_xuic(left_operation) - < op_to_int_xuic(right_operation); - } - else - { - return op_to_int_xiuc(left_operation) - < op_to_int_xiuc(right_operation); - } + if ( left_other_seg_id.ring_index == right_other_seg_id.ring_index ) + { + return op_to_int_xuic(left_operation) + < op_to_int_xuic(right_operation); } else { - //return op_to_int_xuic(left.operations[OpId]) < op_to_int_xuic(right.operations[OpId]); - return left_other_seg_id.multi_index < right_other_seg_id.multi_index; + return op_to_int_xiuc(left_operation) + < op_to_int_xiuc(right_operation); } } }; @@ -217,6 +229,19 @@ struct less_op_areal_areal } }; +template <std::size_t OpId> +struct less_other_multi_index +{ + static const std::size_t other_op_id = (OpId + 1) % 2; + + template <typename Turn> + inline bool operator()(Turn const& left, Turn const& right) const + { + return left.operations[other_op_id].seg_id.multi_index + < right.operations[other_op_id].seg_id.multi_index; + } +}; + // sort turns by G1 - source_index == 0 by: // seg_id -> distance -> operation template <std::size_t OpId = 0, |