diff options
author | DongHun Kwak <dh0128.kwak@samsung.com> | 2016-10-06 10:38:45 +0900 |
---|---|---|
committer | DongHun Kwak <dh0128.kwak@samsung.com> | 2016-10-06 10:39:52 +0900 |
commit | 5cde13f21d36c7224b0e13d11c4b49379ae5210d (patch) | |
tree | e8269ac85a4b0f7d416e2565fa4f451b5cb41351 /boost/geometry/algorithms/detail | |
parent | d9ec475d945d3035377a0d89ed42e382d8988891 (diff) | |
download | boost-5cde13f21d36c7224b0e13d11c4b49379ae5210d.tar.gz boost-5cde13f21d36c7224b0e13d11c4b49379ae5210d.tar.bz2 boost-5cde13f21d36c7224b0e13d11c4b49379ae5210d.zip |
Imported Upstream version 1.61.0
Change-Id: I96a1f878d1e6164f01e9aadd5147f38fca448d90
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>
Diffstat (limited to 'boost/geometry/algorithms/detail')
34 files changed, 2948 insertions, 1876 deletions
diff --git a/boost/geometry/algorithms/detail/andoyer_inverse.hpp b/boost/geometry/algorithms/detail/andoyer_inverse.hpp index eb8485b8e1..66ad6446e9 100644 --- a/boost/geometry/algorithms/detail/andoyer_inverse.hpp +++ b/boost/geometry/algorithms/detail/andoyer_inverse.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2015 Oracle and/or its affiliates. +// Copyright (c) 2015-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -49,23 +49,17 @@ struct andoyer_inverse T2 const& lat2, Spheroid const& spheroid) { + CT const c0 = CT(0); + CT const c1 = CT(1); + CT const pi = math::pi<CT>(); + result_type result; // coordinates in radians - if ( math::equals(lon1, lon2) - && math::equals(lat1, lat2) ) + if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) ) { - result.set(CT(0), CT(0)); - return result; - } - - CT const pi_half = math::pi<CT>() / CT(2); - - if ( math::equals(math::abs(lat1), pi_half) - && math::equals(math::abs(lat2), pi_half) ) - { - result.set(CT(0), CT(0)); + result.set(c0, c0); return result; } @@ -80,22 +74,16 @@ struct andoyer_inverse // H,G,T = infinity if cos_d = 1 or cos_d = -1 // lat1 == +-90 && lat2 == +-90 // lat1 == lat2 && lon1 == lon2 - CT const cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; - CT const d = acos(cos_d); - CT const sin_d = sin(d); - - // just in case since above lat1 and lat2 is checked - // the check below is equal to cos_d == 1 || cos_d == -1 || d == 0 - if ( math::equals(sin_d, CT(0)) ) - { - result.set(CT(0), CT(0)); - return result; - } - - // if the function returned before this place - // and endpoints were on the poles +-90 deg - // in this case the azimuth could either be 0 or +-pi - + CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; + // on some platforms cos_d may be outside valid range + if (cos_d < -c1) + cos_d = -c1; + else if (cos_d > c1) + cos_d = c1; + + CT const d = acos(cos_d); // [0, pi] + CT const sin_d = sin(d); // [-1, 1] + CT const f = detail::flattening<CT>(spheroid); if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) @@ -103,11 +91,18 @@ struct andoyer_inverse CT const K = math::sqr(sin_lat1-sin_lat2); CT const L = math::sqr(sin_lat1+sin_lat2); CT const three_sin_d = CT(3) * sin_d; - // H or G = infinity if cos_d = 1 or cos_d = -1 - CT const H = (d+three_sin_d)/(CT(1)-cos_d); - CT const G = (d-three_sin_d)/(CT(1)+cos_d); - // for e.g. lat1=-90 && lat2=90 here we have G*L=INF*0 + CT const one_minus_cos_d = c1 - cos_d; + CT const one_plus_cos_d = c1 + cos_d; + // cos_d = 1 or cos_d = -1 means that the points are antipodal + + CT const H = math::equals(one_minus_cos_d, c0) ? + c0 : + (d + three_sin_d) / one_minus_cos_d; + CT const G = math::equals(one_plus_cos_d, c0) ? + c0 : + (d - three_sin_d) / one_plus_cos_d; + CT const dd = -(f/CT(4))*(H*K+G*L); CT const a = get_radius<0>(spheroid); @@ -116,41 +111,88 @@ struct andoyer_inverse } else { - result.distance = CT(0); + result.distance = c0; } if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) ) { - CT A = CT(0); - CT U = CT(0); - if ( ! math::equals(cos_lat2, CT(0)) ) + // sin_d = 0 <=> antipodal points + if (math::equals(sin_d, c0)) { - CT const tan_lat2 = sin_lat2/cos_lat2; - CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; - A = atan2(sin_dlon, M); - CT const sin_2A = sin(CT(2)*A); - U = (f/CT(2))*math::sqr(cos_lat1)*sin_2A; + // T = inf + // dA = inf + // azimuth = -inf + if (lat1 <= lat2) + result.azimuth = c0; + else + result.azimuth = pi; } - - CT V = CT(0); - if ( ! math::equals(cos_lat1, CT(0)) ) + else { - CT const tan_lat1 = sin_lat1/cos_lat1; - CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; - CT const B = atan2(sin_dlon, N); - CT const sin_2B = sin(CT(2)*B); - V = (f/CT(2))*math::sqr(cos_lat2)*sin_2B; + CT const c2 = CT(2); + + CT A = c0; + CT U = c0; + if ( ! math::equals(cos_lat2, c0) ) + { + CT const tan_lat2 = sin_lat2/cos_lat2; + CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; + A = atan2(sin_dlon, M); + CT const sin_2A = sin(c2*A); + U = (f/ c2)*math::sqr(cos_lat1)*sin_2A; + } + + CT V = c0; + if ( ! math::equals(cos_lat1, c0) ) + { + CT const tan_lat1 = sin_lat1/cos_lat1; + CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; + CT const B = atan2(sin_dlon, N); + CT const sin_2B = sin(c2*B); + V = (f/ c2)*math::sqr(cos_lat2)*sin_2B; + } + + CT const T = d / sin_d; + CT const dA = V*T-U; + + result.azimuth = A - dA; + + // even with sin_d == 0 checked above if the second point + // is somewhere in the antipodal area T may still be great + // therefore dA may be great and the resulting azimuth + // may be some more or less arbitrary angle + if (A >= c0) // A indicates Eastern hemisphere + { + if (dA >= c0) // A altered towards 0 + { + if ((result.azimuth) < c0) + result.azimuth = c0; + } + else // dA < 0, A altered towards pi + { + if (result.azimuth > pi) + result.azimuth = pi; + } + } + else // A indicates Western hemisphere + { + if (dA <= c0) // A altered towards 0 + { + if (result.azimuth > c0) + result.azimuth = c0; + } + else // dA > 0, A altered towards -pi + { + CT const minus_pi = -pi; + if ((result.azimuth) < minus_pi) + result.azimuth = minus_pi; + } + } } - - // infinity if sin_d = 0, so cos_d = 1 or cos_d = -1 - CT const T = d / sin_d; - CT const dA = V*T-U; - - result.azimuth = A - dA; } else { - result.azimuth = CT(0); + result.azimuth = c0; } return result; diff --git a/boost/geometry/algorithms/detail/as_range.hpp b/boost/geometry/algorithms/detail/as_range.hpp index d0dfb07e43..fef5156eef 100644 --- a/boost/geometry/algorithms/detail/as_range.hpp +++ b/boost/geometry/algorithms/detail/as_range.hpp @@ -15,8 +15,6 @@ #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_AS_RANGE_HPP -#include <boost/type_traits.hpp> - #include <boost/geometry/core/exterior_ring.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> diff --git a/boost/geometry/algorithms/detail/assign_values.hpp b/boost/geometry/algorithms/detail/assign_values.hpp index 5e4a1795b5..ef5691e97c 100644 --- a/boost/geometry/algorithms/detail/assign_values.hpp +++ b/boost/geometry/algorithms/detail/assign_values.hpp @@ -23,7 +23,6 @@ #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/bounds.hpp> #include <boost/numeric/conversion/cast.hpp> -#include <boost/type_traits.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/algorithms/append.hpp> diff --git a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp index 606726f338..8447532a6c 100644 --- a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp @@ -918,6 +918,10 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator collection.discard_rings(); collection.block_turns(); collection.enrich(); + + // phase 1: turns (after enrichment/clustering) + visit_pieces_policy.apply(const_collection, 1); + collection.traverse(); // Reverse all offsetted rings / traversed rings if: @@ -943,8 +947,8 @@ inline void buffer_inserter(GeometryInput const& geometry_input, OutputIterator collection.template assign<GeometryOutput>(out); // Visit collection again - // phase 1: rings (after discarding and traversing) - visit_pieces_policy.apply(const_collection, 1); + // phase 2: rings (after traversing) + visit_pieces_policy.apply(const_collection, 2); } template diff --git a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp index 255ac797dc..c0fa7b0e33 100644 --- a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp @@ -46,15 +46,26 @@ class backtrack_for_buffer public : typedef detail::overlay::backtrack_state state_type; - template <typename Operation, typename Rings, typename Turns, typename Geometry, typename RobustPolicy> + template + < + typename Operation, + typename Rings, + typename Turns, + typename Geometry, + typename RobustPolicy, + typename Visitor + > static inline void apply(std::size_t size_at_start, Rings& rings, typename boost::range_value<Rings>::type& ring, - Turns& turns, Operation& operation, - std::string const& /*reason*/, + Turns& turns, + typename boost::range_value<Turns>::type const& turn, + Operation& operation, + detail::overlay::traverse_error_type traverse_error, Geometry const& , Geometry const& , RobustPolicy const& , - state_type& state + state_type& state, + Visitor& visitor ) { #if defined(BOOST_GEOMETRY_COUNT_BACKTRACK_WARNINGS) @@ -62,7 +73,7 @@ extern int g_backtrack_warning_count; g_backtrack_warning_count++; #endif //std::cout << "!"; -//std::cout << "WARNING " << reason << std::endl; +//std::cout << "WARNING " << traverse_error_string(traverse_error) << std::endl; state.m_good = false; @@ -78,6 +89,41 @@ g_backtrack_warning_count++; } }; +struct buffer_overlay_visitor +{ +public : + void print(char const* header) + { + } + + template <typename Turns> + void print(char const* header, Turns const& turns, int turn_index) + { + } + + template <typename Turns> + void print(char const* header, Turns const& turns, int turn_index, int op_index) + { + } + + template <typename Turns> + void visit_turns(int , Turns const& ) {} + + template <typename Clusters, typename Turns> + void visit_clusters(Clusters const& , Turns const& ) {} + + template <typename Turns, typename Turn, typename Operation> + void visit_traverse(Turns const& turns, Turn const& turn, Operation const& op, const char* header) + { + } + + template <typename Turns, typename Turn, typename Operation> + void visit_traverse_reject(Turns const& , Turn const& , Operation const& , + detail::overlay::traverse_error_type ) + {} +}; + + // Should follow traversal-turn-concept (enrichment, visit structure) // and adds index in piece vector to find it back template <typename Point, typename SegmentRatio> diff --git a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp index b580cf5b9b..3fc3c2347a 100644 --- a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp @@ -2,6 +2,10 @@ // Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2016. +// Modifications copyright (c) 2016 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) @@ -255,8 +259,11 @@ struct buffered_piece_collection { geometry::envelope(m_ring, m_box); - // create monotonic sections in y-dimension - typedef boost::mpl::vector_c<std::size_t, 1> dimensions; + // create monotonic sections in x-dimension + // The dimension is critical because the direction is later used + // in the optimization for within checks using winding strategy + // and this strategy is scanning in x direction. + typedef boost::mpl::vector_c<std::size_t, 0> dimensions; geometry::sectionalize<false, dimensions>(m_ring, detail::no_rescale_policy(), m_sections); } @@ -286,6 +293,15 @@ struct buffered_piece_collection typedef geometry::sections<robust_box_type, 2> sections_type; sections_type monotonic_sections; + // Define the clusters, mapping cluster_id -> turns + typedef std::map + < + signed_size_type, + std::set<signed_size_type> + > cluster_type; + + cluster_type m_clusters; + RobustPolicy const& m_robust_policy; @@ -479,7 +495,8 @@ struct buffered_piece_collection // Within can have in rare cases a rounding issue. We don't discard this // point, so it can be used to continue started rings in traversal. But // will never start a new ring from this type of points. - it->selectable_start = false; + it->operations[0].enriched.startable = false; + it->operations[1].enriched.startable = false; } #endif } @@ -931,8 +948,12 @@ struct buffered_piece_collection piece pc; pc.type = type; pc.index = static_cast<signed_size_type>(boost::size(m_pieces)); + + current_segment_id.piece_index = pc.index; + pc.first_seg_id = current_segment_id; + // Assign left/right (for first/last piece per ring they will be re-assigned later) pc.left_index = pc.index - 1; pc.right_index = pc.index + 1; @@ -1196,7 +1217,7 @@ struct buffered_piece_collection >::type side_strategy_type; enrich_intersection_points<false, false, overlay_union>(m_turns, - detail::overlay::operation_union, + m_clusters, detail::overlay::operation_union, offsetted_rings, offsetted_rings, m_robust_policy, side_strategy_type()); } @@ -1213,7 +1234,7 @@ struct buffered_piece_collection offsetted_rings[it->operations[0].seg_id.multi_index].has_discarded_intersections = true; offsetted_rings[it->operations[1].seg_id.multi_index].has_discarded_intersections = true; } - else if (! it->both(detail::overlay::operation_union)) + else { offsetted_rings[it->operations[0].seg_id.multi_index].has_accepted_intersections = true; offsetted_rings[it->operations[1].seg_id.multi_index].has_accepted_intersections = true; @@ -1331,13 +1352,15 @@ struct buffered_piece_collection false, false, buffered_ring_collection<buffered_ring<Ring> >, buffered_ring_collection<buffered_ring<Ring > >, + detail::overlay::operation_union, backtrack_for_buffer > traverser; traversed_rings.clear(); + buffer_overlay_visitor visitor; traverser::apply(offsetted_rings, offsetted_rings, - detail::overlay::operation_union, - m_robust_policy, m_turns, traversed_rings); + m_robust_policy, m_turns, traversed_rings, + m_clusters, visitor); } inline void reverse() diff --git a/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp b/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp index 29a618b923..db6136e1ae 100644 --- a/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp +++ b/boost/geometry/algorithms/detail/buffer/buffered_ring.hpp @@ -149,6 +149,14 @@ struct ring_type } + +template <> +struct single_tag_of<detail::buffer::buffered_ring_collection_tag> +{ + typedef ring_tag type; +}; + + namespace dispatch { @@ -213,6 +221,21 @@ struct within }; +template <typename Geometry> +struct is_empty<Geometry, detail::buffer::buffered_ring_collection_tag> + : detail::is_empty::multi_is_empty<detail::is_empty::range_is_empty> +{}; + + +template <typename Geometry> +struct envelope<Geometry, detail::buffer::buffered_ring_collection_tag> + : detail::envelope::envelope_multi_range + < + detail::envelope::envelope_range + > +{}; + + } // namespace dispatch namespace detail { namespace overlay diff --git a/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp b/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp index 05fc6df8ff..da1d084d32 100644 --- a/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp +++ b/boost/geometry/algorithms/detail/buffer/turn_in_original_visitor.hpp @@ -2,6 +2,10 @@ // Copyright (c) 2014 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2016. +// Modifications copyright (c) 2016 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) @@ -103,32 +107,32 @@ template typename Iterator > inline bool point_in_section(Strategy& strategy, State& state, - Point const& point, CoordinateType const& point_y, + Point const& point, CoordinateType const& point_x, Iterator begin, Iterator end, int direction) { if (direction == 0) { - // Not a monotonic section, or no change in Y-direction + // Not a monotonic section, or no change in X-direction return point_in_range(strategy, state, point, begin, end); } - // We're in a monotonic section in y-direction + // We're in a monotonic section in x-direction Iterator it = begin; for (Iterator previous = it++; it != end; ++previous, ++it) { // Depending on sections.direction we can quit for this section - CoordinateType const previous_y = geometry::get<1>(*previous); + CoordinateType const previous_x = geometry::get<0>(*previous); - if (direction == 1 && point_y < previous_y) + if (direction == 1 && point_x < previous_x) { - // Section goes upwards, y increases, point is is below section + // Section goes upwards, x increases, point is is below section return true; } - else if (direction == -1 && point_y > previous_y) + else if (direction == -1 && point_x > previous_x) { - // Section goes downwards, y decreases, point is above section + // Section goes downwards, x decreases, point is above section return true; } @@ -145,6 +149,9 @@ inline bool point_in_section(Strategy& strategy, State& state, template <typename Point, typename Original> inline int point_in_original(Point const& point, Original const& original) { + // The winding strategy is scanning in x direction + // therefore it's critical to pass direction calculated + // for x dimension below. typedef strategy::within::winding<Point> strategy_type; typename strategy_type::state_type state; @@ -166,7 +173,7 @@ inline int point_in_original(Point const& point, Original const& original) typedef typename boost::range_value<sections_type const>::type section_type; typedef typename geometry::coordinate_type<Point>::type coordinate_type; - coordinate_type const point_y = geometry::get<1>(point); + coordinate_type const point_x = geometry::get<0>(point); // Walk through all monotonic sections of this original for (iterator_type it = boost::begin(original.m_sections); @@ -177,11 +184,11 @@ inline int point_in_original(Point const& point, Original const& original) if (! section.duplicate && section.begin_index < section.end_index - && point_y >= geometry::get<min_corner, 1>(section.bounding_box) - && point_y <= geometry::get<max_corner, 1>(section.bounding_box)) + && point_x >= geometry::get<min_corner, 0>(section.bounding_box) + && point_x <= geometry::get<max_corner, 0>(section.bounding_box)) { - // y-coordinate of point overlaps with section - if (! point_in_section(strategy, state, point, point_y, + // x-coordinate of point overlaps with section + if (! point_in_section(strategy, state, point, point_x, boost::begin(original.m_ring) + section.begin_index, boost::begin(original.m_ring) + section.end_index + 1, section.directions[0])) diff --git a/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp b/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp index 0aca21ce88..d7146befb9 100644 --- a/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp +++ b/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp @@ -2,6 +2,10 @@ // Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2016. +// Modifications copyright (c) 2016 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) @@ -218,16 +222,16 @@ public : BOOST_GEOMETRY_ASSERT(! piece.sections.empty()); - coordinate_type const point_y = geometry::get<1>(turn.robust_point); + coordinate_type const point_x = geometry::get<0>(turn.robust_point); for (std::size_t s = 0; s < piece.sections.size(); s++) { section_type const& section = piece.sections[s]; - // If point within vertical range of monotonic section: + // If point within horizontal range of monotonic section: if (! section.duplicate && section.begin_index < section.end_index - && point_y >= geometry::get<min_corner, 1>(section.bounding_box) - 1 - && point_y <= geometry::get<max_corner, 1>(section.bounding_box) + 1) + && point_x >= geometry::get<min_corner, 0>(section.bounding_box) - 1 + && point_x <= geometry::get<max_corner, 0>(section.bounding_box) + 1) { for (signed_size_type i = section.begin_index + 1; i <= section.end_index; i++) { @@ -239,21 +243,21 @@ public : // First check if it is in range - if it is not, the // expensive side_of_intersection does not need to be // applied - coordinate_type y1 = geometry::get<1>(previous); - coordinate_type y2 = geometry::get<1>(current); + coordinate_type x1 = geometry::get<0>(previous); + coordinate_type x2 = geometry::get<0>(current); - if (y1 > y2) + if (x1 > x2) { - std::swap(y1, y2); + std::swap(x1, x2); } - if (point_y >= y1 - 1 && point_y <= y2 + 1) + if (point_x >= x1 - 1 && point_x <= x2 + 1) { segment_type const r(previous, current); int const side = strategy::side::side_of_intersection::apply(p, q, r, turn.robust_point); - // Sections are monotonic in y-dimension + // Sections are monotonic in x-dimension if (side == 1) { // Left on segment diff --git a/boost/geometry/algorithms/detail/direction_code.hpp b/boost/geometry/algorithms/detail/direction_code.hpp new file mode 100644 index 0000000000..26d53ab4e5 --- /dev/null +++ b/boost/geometry/algorithms/detail/direction_code.hpp @@ -0,0 +1,79 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015 Oracle and/or its affiliates. + +// 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) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECITON_CODE_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECITON_CODE_HPP + + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/util/math.hpp> + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <std::size_t Index, typename Point1, typename Point2> +inline int sign_of_difference(Point1 const& point1, Point2 const& point2) +{ + return + math::equals(geometry::get<Index>(point1), geometry::get<Index>(point2)) + ? + 0 + : + (geometry::get<Index>(point1) > geometry::get<Index>(point2) ? 1 : -1); +} + + +// Gives sense of direction for point p, collinear w.r.t. segment (a,b) +// Returns -1 if p goes backward w.r.t (a,b), so goes from b in direction of a +// Returns 1 if p goes forward, so extends (a,b) +// Returns 0 if p is equal with b, or if (a,b) is degenerate +// Note that it does not do any collinearity test, that should be done before +template <typename Point1, typename Point2> +inline int direction_code(Point1 const& segment_a, Point1 const& segment_b, + const Point2& p) +{ + // Suppose segment = (4 3,4 4) and p =(4 2) + // Then sign_a1 = 1 and sign_p1 = 1 -> goes backward -> return -1 + + int const sign_a0 = sign_of_difference<0>(segment_b, segment_a); + int const sign_a1 = sign_of_difference<1>(segment_b, segment_a); + + if (sign_a0 == 0 && sign_a1 == 0) + { + return 0; + } + + int const sign_p0 = sign_of_difference<0>(segment_b, p); + int const sign_p1 = sign_of_difference<1>(segment_b, p); + + if (sign_p0 == 0 && sign_p1 == 0) + { + return 0; + } + + return sign_a0 == sign_p0 && sign_a1 == sign_p1 ? -1 : 1; +} + + +} // namespace detail +#endif //DOXYGEN_NO_DETAIL + + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_DIRECITON_CODE_HPP diff --git a/boost/geometry/algorithms/detail/disjoint/box_box.hpp b/boost/geometry/algorithms/detail/disjoint/box_box.hpp index 6074af982b..3b81755e20 100644 --- a/boost/geometry/algorithms/detail/disjoint/box_box.hpp +++ b/boost/geometry/algorithms/detail/disjoint/box_box.hpp @@ -5,8 +5,8 @@ // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. // Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2013-2015. -// Modifications copyright (c) 2013-2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2013-2016. +// Modifications copyright (c) 2013-2016, 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 @@ -28,6 +28,9 @@ #include <boost/geometry/algorithms/dispatch/disjoint.hpp> +#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + namespace boost { namespace geometry { @@ -39,7 +42,13 @@ namespace detail { namespace disjoint template < typename Box1, typename Box2, - std::size_t Dimension, std::size_t DimensionCount + std::size_t Dimension = 0, + std::size_t DimensionCount = dimension<Box1>::value, + typename CSTag = typename tag_cast + < + typename cs_tag<Box1>::type, + spherical_tag + >::type > struct box_box { @@ -62,8 +71,8 @@ struct box_box }; -template <typename Box1, typename Box2, std::size_t DimensionCount> -struct box_box<Box1, Box2, DimensionCount, DimensionCount> +template <typename Box1, typename Box2, std::size_t DimensionCount, typename CSTag> +struct box_box<Box1, Box2, DimensionCount, DimensionCount, CSTag> { static inline bool apply(Box1 const& , Box2 const& ) { @@ -72,6 +81,51 @@ struct box_box<Box1, Box2, DimensionCount, DimensionCount> }; +template <typename Box1, typename Box2, std::size_t DimensionCount> +struct box_box<Box1, Box2, 0, DimensionCount, spherical_tag> +{ + static inline bool apply(Box1 const& box1, Box2 const& box2) + { + typedef typename geometry::select_most_precise + < + typename coordinate_type<Box1>::type, + typename coordinate_type<Box2>::type + >::type calc_t; + typedef typename coordinate_system<Box1>::type::units units_t; + typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; + + calc_t const b1_min = get<min_corner, 0>(box1); + calc_t const b1_max = get<max_corner, 0>(box1); + calc_t const b2_min = get<min_corner, 0>(box2); + calc_t const b2_max = get<max_corner, 0>(box2); + + // min <= max <=> diff >= 0 + calc_t const diff1 = b1_max - b1_min; + calc_t const diff2 = b2_max - b2_min; + + // check the intersection if neither box cover the whole globe + if (diff1 < constants::period() && diff2 < constants::period()) + { + // calculate positive longitude translation with b1_min as origin + calc_t const diff_min = math::longitude_distance_unsigned<units_t>(b1_min, b2_min); + calc_t const b2_min_transl = b1_min + diff_min; // always right of b1_min + + if (b2_min_transl > b1_max // b2_min right of b1_max + && b2_min_transl - constants::period() + diff2 < b1_min) // b2_max left of b1_min + { + return true; + } + } + + return box_box + < + Box1, Box2, + 1, DimensionCount + >::apply(box1, box2); + } +}; + + /*! \brief Internal utility function to detect if boxes are disjoint \note Is used from other algorithms, declared separately @@ -80,11 +134,7 @@ struct box_box<Box1, Box2, DimensionCount, DimensionCount> template <typename Box1, typename Box2> inline bool disjoint_box_box(Box1 const& box1, Box2 const& box2) { - return box_box - < - Box1, Box2, - 0, dimension<Box1>::type::value - >::apply(box1, box2); + return box_box<Box1, Box2>::apply(box1, box2); } diff --git a/boost/geometry/algorithms/detail/disjoint/point_box.hpp b/boost/geometry/algorithms/detail/disjoint/point_box.hpp index 73b7b70990..2f1085ada9 100644 --- a/boost/geometry/algorithms/detail/disjoint/point_box.hpp +++ b/boost/geometry/algorithms/detail/disjoint/point_box.hpp @@ -5,8 +5,8 @@ // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. // Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland -// This file was modified by Oracle on 2013-2015. -// Modifications copyright (c) 2013-2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2013-2016. +// Modifications copyright (c) 2013-2016, 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 @@ -28,6 +28,7 @@ #include <boost/geometry/core/tags.hpp> #include <boost/geometry/algorithms/dispatch/disjoint.hpp> +#include <boost/geometry/strategies/cartesian/point_in_box.hpp> namespace boost { namespace geometry { @@ -37,49 +38,19 @@ namespace detail { namespace disjoint { -template -< - typename Point, typename Box, - std::size_t Dimension, std::size_t DimensionCount -> -struct point_box -{ - static inline bool apply(Point const& point, Box const& box) - { - if (get<Dimension>(point) < get<min_corner, Dimension>(box) - || get<Dimension>(point) > get<max_corner, Dimension>(box)) - { - return true; - } - return point_box - < - Point, Box, - Dimension + 1, DimensionCount - >::apply(point, box); - } -}; - - -template <typename Point, typename Box, std::size_t DimensionCount> -struct point_box<Point, Box, DimensionCount, DimensionCount> -{ - static inline bool apply(Point const& , Box const& ) - { - return false; - } -}; - /*! \brief Internal utility function to detect if point/box are disjoint */ template <typename Point, typename Box> inline bool disjoint_point_box(Point const& point, Box const& box) { - return detail::disjoint::point_box - < - Point, Box, - 0, dimension<Point>::type::value - >::apply(point, box); + // ! covered_by(point, box) + return ! strategy::within::relate_point_box_loop + < + strategy::within::covered_by_range, + Point, Box, + 0, dimension<Point>::type::value + >::apply(point, box); } @@ -94,8 +65,18 @@ namespace dispatch template <typename Point, typename Box, std::size_t DimensionCount> struct disjoint<Point, Box, DimensionCount, point_tag, box_tag, false> - : detail::disjoint::point_box<Point, Box, 0, DimensionCount> -{}; +{ + static inline bool apply(Point const& point, Box const& box) + { + // ! covered_by(point, box) + return ! strategy::within::relate_point_box_loop + < + strategy::within::covered_by_range, + Point, Box, + 0, DimensionCount + >::apply(point, box); + } +}; } // namespace dispatch diff --git a/boost/geometry/algorithms/detail/disjoint/segment_box.hpp b/boost/geometry/algorithms/detail/disjoint/segment_box.hpp index 5368432ed4..cc0c7949e3 100644 --- a/boost/geometry/algorithms/detail/disjoint/segment_box.hpp +++ b/boost/geometry/algorithms/detail/disjoint/segment_box.hpp @@ -128,7 +128,9 @@ struct disjoint_segment_box_impl && t_min.first > ti_max) || (geometry::math::equals(t_max.second, 0) - && t_max.first < ti_min) ) + && t_max.first < ti_min) + || + (math::sign(ti_min) * math::sign(ti_max) > 0) ) { return true; } @@ -197,6 +199,11 @@ struct disjoint_segment_box_impl { if ( geometry::math::equals(t_min.first, 0) ) { t_min.first = -1; } if ( geometry::math::equals(t_max.first, 0) ) { t_max.first = 1; } + + if (math::sign(t_min.first) * math::sign(t_max.first) > 0) + { + return true; + } } if ( t_min.first > diff || t_max.first < 0 ) diff --git a/boost/geometry/algorithms/detail/envelope/segment.hpp b/boost/geometry/algorithms/detail/envelope/segment.hpp index 570f0e1a43..6186e72a3a 100644 --- a/boost/geometry/algorithms/detail/envelope/segment.hpp +++ b/boost/geometry/algorithms/detail/envelope/segment.hpp @@ -4,10 +4,11 @@ // Copyright (c) 2008-2015 Bruno Lalande, Paris, France. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2015. -// Modifications copyright (c) 2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2015, 2016. +// Modifications copyright (c) 2015-2016, Oracle and/or its affiliates. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at @@ -84,6 +85,7 @@ class compute_mbr_of_segment private: // computes the azimuths of the segment with endpoints (lon1, lat1) // and (lon2, lat2) + // radians template <typename CalculationType> static inline void azimuths(CalculationType const& lon1, CalculationType const& lat1, @@ -92,7 +94,7 @@ private: CalculationType& a1, CalculationType& a2) { - BOOST_GEOMETRY_ASSERT(math::smaller(lon1, lon2)); + BOOST_GEOMETRY_ASSERT(lon1 <= lon2); CalculationType dlon = lon2 - lon1; CalculationType sin_dlon = sin(dlon); @@ -108,8 +110,9 @@ private: a2 = atan2(-sin_dlon * cos_lat1, cos_lat2 * sin_lat1 - sin_lat2 * cos_lat1 * cos_dlon); a2 += math::pi<CalculationType>(); - } + } + // degrees or radians template <typename CalculationType> static inline void swap(CalculationType& lon1, CalculationType& lat1, @@ -120,6 +123,7 @@ private: std::swap(lat1, lat2); } + // radians template <typename CalculationType> static inline bool contains_pi_half(CalculationType const& a1, CalculationType const& a2) @@ -134,13 +138,20 @@ private: : (a1 > pi_half && pi_half > a2); } - template <typename CoordinateType> + // radians or degrees + template <typename Units, typename CoordinateType> static inline bool crosses_antimeridian(CoordinateType const& lon1, CoordinateType const& lon2) { - return math::larger(math::abs(lon1 - lon2), math::pi<CoordinateType>()); + typedef math::detail::constants_on_spheroid + < + CoordinateType, Units + > constants; + + return math::abs(lon1 - lon2) > constants::half_period(); // > pi } + // radians template <typename CalculationType> static inline CalculationType max_latitude(CalculationType const& azimuth, CalculationType const& latitude) @@ -149,39 +160,46 @@ private: return acos( math::abs(cos(latitude) * sin(azimuth)) ); } - template <typename CalculationType> + // degrees or radians + template <typename Units, typename CalculationType> static inline void compute_box_corners(CalculationType& lon1, CalculationType& lat1, CalculationType& lon2, - CalculationType& lat2, - CalculationType const& a1, - CalculationType const& a2) + CalculationType& lat2) { // coordinates are assumed to be in radians - BOOST_GEOMETRY_ASSERT(! math::larger(lon1, lon2)); + BOOST_GEOMETRY_ASSERT(lon1 <= lon2); - if (math::equals(a1, a2)) + CalculationType lon1_rad = math::as_radian<Units>(lon1); + CalculationType lat1_rad = math::as_radian<Units>(lat1); + CalculationType lon2_rad = math::as_radian<Units>(lon2); + CalculationType lat2_rad = math::as_radian<Units>(lat2); + + CalculationType a1 = 0, a2 = 0; + azimuths(lon1_rad, lat1_rad, lon2_rad, lat2_rad, a1, a2); + + if (lat1 > lat2) { - // the segment must lie on the equator; nothing to do - BOOST_GEOMETRY_ASSERT(math::equals(lat1, CalculationType(0))); - BOOST_GEOMETRY_ASSERT(math::equals(lat2, CalculationType(0))); - return; + std::swap(lat1, lat2); + std::swap(lat1_rad, lat2_rad); } - if (math::larger(lat1, lat2)) + if (math::equals(a1, a2)) { - std::swap(lat1, lat2); + // the segment must lie on the equator or is very short + return; } if (contains_pi_half(a1, a2)) { - CalculationType mid_lat = lat1 + lat2; + CalculationType const mid_lat = lat1 + lat2; if (mid_lat < 0) { // update using min latitude - CalculationType lat_min = -max_latitude(a1, lat1); + CalculationType const lat_min_rad = -max_latitude(a1, lat1_rad); + CalculationType const lat_min = math::from_radian<Units>(lat_min_rad); - if (math::larger(lat1, lat_min)) + if (lat1 > lat_min) { lat1 = lat_min; } @@ -189,9 +207,10 @@ private: else if (mid_lat > 0) { // update using max latitude - CalculationType lat_max = max_latitude(a1, lat1); + CalculationType const lat_max_rad = max_latitude(a1, lat1_rad); + CalculationType const lat_max = math::from_radian<Units>(lat_max_rad); - if (math::smaller(lat2, lat_max)) + if (lat2 < lat_max) { lat2 = lat_max; } @@ -199,24 +218,27 @@ private: } } - template <typename CalculationType> + template <typename Units, typename CalculationType> static inline void apply(CalculationType& lon1, CalculationType& lat1, CalculationType& lon2, CalculationType& lat2) { - CalculationType const half_pi = math::half_pi<CalculationType>(); + typedef math::detail::constants_on_spheroid + < + CalculationType, Units + > constants; - bool is_pole1 = math::equals(math::abs(lat1), half_pi); - bool is_pole2 = math::equals(math::abs(lat2), half_pi); + bool is_pole1 = math::equals(math::abs(lat1), constants::max_latitude()); + bool is_pole2 = math::equals(math::abs(lat2), constants::max_latitude()); if (is_pole1 && is_pole2) { // both points are poles; nothing more to do: // longitudes are already normalized to 0 - BOOST_GEOMETRY_ASSERT(lon1 == CalculationType(0) - && - lon2 == CalculationType(0)); + // but just in case + lon1 = 0; + lon2 = 0; } else if (is_pole1 && !is_pole2) { @@ -233,10 +255,10 @@ private: lon2 = lon1; } - if (math::equals(lon1, lon2)) + if (lon1 == lon2) { // segment lies on a meridian - if (math::larger(lat1, lat2)) + if (lat1 > lat2) { std::swap(lat1, lat2); } @@ -245,25 +267,22 @@ private: BOOST_GEOMETRY_ASSERT(!is_pole1 && !is_pole2); - if (math::larger(lon1, lon2)) + if (lon1 > lon2) { swap(lon1, lat1, lon2, lat2); } - if (crosses_antimeridian(lon1, lon2)) + if (crosses_antimeridian<Units>(lon1, lon2)) { - lon1 += math::two_pi<CalculationType>(); + lon1 += constants::period(); swap(lon1, lat1, lon2, lat2); } - CalculationType a1 = 0, a2 = 0; - azimuths(lon1, lat1, lon2, lat2, a1, a2); - - compute_box_corners(lon1, lat1, lon2, lat2, a1, a2); + compute_box_corners<Units>(lon1, lat1, lon2, lat2); } public: - template <typename CalculationType, typename Box> + template <typename Units, typename CalculationType, typename Box> static inline void apply(CalculationType lon1, CalculationType lat1, CalculationType lon2, @@ -274,12 +293,12 @@ public: typedef typename helper_geometry < - Box, box_coordinate_type, radian + Box, box_coordinate_type, Units >::type helper_box_type; helper_box_type radian_mbr; - apply(lon1, lat1, lon2, lat2); + apply<Units>(lon1, lat1, lon2, lat2); geometry::set < @@ -316,11 +335,14 @@ struct envelope_segment_on_sphere Point p1_normalized = detail::return_normalized<Point>(p1); Point p2_normalized = detail::return_normalized<Point>(p2); - compute_mbr_of_segment::apply(geometry::get_as_radian<0>(p1_normalized), - geometry::get_as_radian<1>(p1_normalized), - geometry::get_as_radian<0>(p2_normalized), - geometry::get_as_radian<1>(p2_normalized), - mbr); + typedef typename coordinate_system<Point>::type::units units_type; + + compute_mbr_of_segment::template apply<units_type>( + geometry::get<0>(p1_normalized), + geometry::get<1>(p1_normalized), + geometry::get<0>(p2_normalized), + geometry::get<1>(p2_normalized), + mbr); // now compute the envelope range for coordinates of // dimension 2 and higher diff --git a/boost/geometry/algorithms/detail/intersection/interface.hpp b/boost/geometry/algorithms/detail/intersection/interface.hpp index 2af618d974..fd989866dd 100644 --- a/boost/geometry/algorithms/detail/intersection/interface.hpp +++ b/boost/geometry/algorithms/detail/intersection/interface.hpp @@ -54,7 +54,7 @@ struct intersection < Geometry1, Geometry2, OneOut, overlay_intersection - >::apply(geometry1, geometry2, robust_policy, std::back_inserter(geometry_out), strategy); + >::apply(geometry1, geometry2, robust_policy, range::back_inserter(geometry_out), strategy); return true; } diff --git a/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp b/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp index 285edfdd6c..03c06c28d1 100644 --- a/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp +++ b/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp @@ -72,7 +72,7 @@ inline void append_no_dups_or_spikes(Range& range, Point const& point, << geometry::get<0>(point) << ", " << geometry::get<1>(point) << ")" << std::endl; #endif - // The code below thies condition checks all spikes/dups + // The code below this condition checks all spikes/dups // for geometries >= 3 points. // So we have to check the first potential duplicate differently if (boost::size(range) == 1 diff --git a/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp b/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp index 90901dee70..a8171e1482 100644 --- a/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp +++ b/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp @@ -48,7 +48,6 @@ inline void clear_visit_info(Turns& turns) { op_it->visited.clear(); } - it->discarded = false; } } @@ -62,6 +61,34 @@ struct backtrack_state }; +enum traverse_error_type +{ + traverse_error_none, + traverse_error_no_next_ip_at_start, + traverse_error_no_next_ip, + traverse_error_dead_end_at_start, + traverse_error_dead_end, + traverse_error_visit_again, + traverse_error_endless_loop +}; + +inline std::string traverse_error_string(traverse_error_type error) +{ + switch (error) + { + case traverse_error_none : return ""; + case traverse_error_no_next_ip_at_start : return "No next IP at start"; + case traverse_error_no_next_ip : return "No next IP"; + case traverse_error_dead_end_at_start : return "Dead end at start"; + case traverse_error_dead_end : return "Dead end"; + case traverse_error_visit_again : return "Visit again"; + case traverse_error_endless_loop : return "Endless loop"; + default : return ""; + } + return ""; +} + + template < typename Geometry1, @@ -73,23 +100,28 @@ class backtrack_check_self_intersections { bool m_checked; inline state() - : m_checked() + : m_checked(true) {} }; public : typedef state state_type; - template <typename Operation, typename Rings, typename Ring, typename Turns, typename RobustPolicy> + template <typename Operation, typename Rings, typename Ring, typename Turns, typename RobustPolicy, typename Visitor> static inline void apply(std::size_t size_at_start, Rings& rings, Ring& ring, - Turns& turns, Operation& operation, - std::string const& , + Turns& turns, + typename boost::range_value<Turns>::type const& turn, + Operation& operation, + traverse_error_type traverse_error, Geometry1 const& geometry1, Geometry2 const& geometry2, RobustPolicy const& robust_policy, - state_type& state + state_type& state, + Visitor& visitor ) { + visitor.visit_traverse_reject(turns, turn, operation, traverse_error); + state.m_good = false; // Check self-intersections and throw exception if appropriate diff --git a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp index dbf3770357..790819917c 100644 --- a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp +++ b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp @@ -22,7 +22,9 @@ #include <boost/geometry/algorithms/convert.hpp> #include <boost/geometry/geometries/concepts/check.hpp> #include <boost/geometry/util/range.hpp> -#include <boost/geometry/views/detail/normalized_view.hpp> +#include <boost/geometry/iterators/ever_circling_iterator.hpp> +#include <boost/geometry/views/closeable_view.hpp> +#include <boost/geometry/views/reversible_view.hpp> namespace boost { namespace geometry @@ -38,25 +40,33 @@ template <typename Range, bool Reverse, typename SegmentIdentifier, typename Poi struct copy_segment_point_range { static inline bool apply(Range const& range, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point) { - detail::normalized_view<Range const> view(range); + typedef typename closeable_view + < + Range const, + closure<Range>::value + >::type cview_type; - signed_size_type const n = boost::size(view); - signed_size_type index = seg_id.segment_index; - if (second) + typedef typename reversible_view + < + cview_type const, + Reverse ? iterate_reverse : iterate_forward + >::type rview_type; + + cview_type cview(range); + rview_type view(cview); + + typedef typename boost::range_iterator<rview_type>::type iterator; + geometry::ever_circling_iterator<iterator> it(boost::begin(view), boost::end(view), + boost::begin(view) + seg_id.segment_index, true); + + for (signed_size_type i = 0; i < offset; ++i, ++it) { - index++; - if (index >= n) - { - index = 0; - } } - BOOST_GEOMETRY_ASSERT(index >= 0 && index < n); - - geometry::convert(*(boost::begin(view) + index), point); + geometry::convert(*it, point); return true; } }; @@ -66,7 +76,7 @@ template <typename Polygon, bool Reverse, typename SegmentIdentifier, typename P struct copy_segment_point_polygon { static inline bool apply(Polygon const& polygon, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point) { // Call ring-version with the right ring @@ -81,7 +91,7 @@ struct copy_segment_point_polygon seg_id.ring_index < 0 ? geometry::exterior_ring(polygon) : range::at(geometry::interior_rings(polygon), seg_id.ring_index), - seg_id, second, + seg_id, offset, point ); } @@ -92,11 +102,11 @@ template <typename Box, bool Reverse, typename SegmentIdentifier, typename Point struct copy_segment_point_box { static inline bool apply(Box const& box, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point) { signed_size_type index = seg_id.segment_index; - if (second) + for (int i = 0; i < offset; i++) { index++; } @@ -119,7 +129,7 @@ template struct copy_segment_point_multi { static inline bool apply(MultiGeometry const& multi, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point) { @@ -130,7 +140,7 @@ struct copy_segment_point_multi ); // Call the single-version - return Policy::apply(range::at(multi, seg_id.multi_index), seg_id, second, point); + return Policy::apply(range::at(multi, seg_id.multi_index), seg_id, offset, point); } }; @@ -271,7 +281,7 @@ struct copy_segment_point */ template<bool Reverse, typename Geometry, typename SegmentIdentifier, typename PointOut> inline bool copy_segment_point(Geometry const& geometry, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point_out) { concept::check<Geometry const>(); @@ -283,7 +293,7 @@ inline bool copy_segment_point(Geometry const& geometry, Reverse, SegmentIdentifier, PointOut - >::apply(geometry, seg_id, second, point_out); + >::apply(geometry, seg_id, offset, point_out); } @@ -300,7 +310,7 @@ template typename PointOut > inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geometry2, - SegmentIdentifier const& seg_id, bool second, + SegmentIdentifier const& seg_id, int offset, PointOut& point_out) { concept::check<Geometry1 const>(); @@ -317,7 +327,7 @@ inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geom Reverse1, SegmentIdentifier, PointOut - >::apply(geometry1, seg_id, second, point_out); + >::apply(geometry1, seg_id, offset, point_out); } else if (seg_id.source_index == 1) { @@ -328,7 +338,7 @@ inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geom Reverse2, SegmentIdentifier, PointOut - >::apply(geometry2, seg_id, second, point_out); + >::apply(geometry2, seg_id, offset, point_out); } // Exception? return false; @@ -354,10 +364,33 @@ inline bool copy_segment_points(Geometry1 const& geometry1, Geometry2 const& geo concept::check<Geometry1 const>(); concept::check<Geometry2 const>(); - return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, false, point1) - && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, true, point2); + return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 0, point1) + && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 1, point2); } +/*! + \brief Helper function, copies three points: two from the specified segment + (from, to) and the next one + \ingroup overlay + */ +template +< + bool Reverse1, bool Reverse2, + typename Geometry1, typename Geometry2, + typename SegmentIdentifier, + typename PointOut +> +inline bool copy_segment_points(Geometry1 const& geometry1, Geometry2 const& geometry2, + SegmentIdentifier const& seg_id, + PointOut& point1, PointOut& point2, PointOut& point3) +{ + concept::check<Geometry1 const>(); + concept::check<Geometry2 const>(); + + return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 0, point1) + && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 1, point2) + && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 2, point3); +} diff --git a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp index bc84286241..c79c9cf500 100644 --- a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp +++ b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp @@ -24,17 +24,20 @@ #include <boost/range.hpp> +#include <boost/geometry/iterators/ever_circling_iterator.hpp> #include <boost/geometry/algorithms/detail/ring_identifier.hpp> #include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp> #include <boost/geometry/algorithms/detail/overlay/handle_colocations.hpp> -#include <boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp> +#include <boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp> #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp> +#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> #include <boost/geometry/policies/robustness/robust_type.hpp> #include <boost/geometry/strategies/side.hpp> #ifdef BOOST_GEOMETRY_DEBUG_ENRICH # include <boost/geometry/algorithms/detail/overlay/check_enrich.hpp> #endif + namespace boost { namespace geometry { @@ -42,183 +45,6 @@ namespace boost { namespace geometry namespace detail { namespace overlay { -// Wraps "turn_operation" from turn_info.hpp, -// giving it extra information -template <typename TurnOperation> -struct indexed_turn_operation -{ - typedef TurnOperation type; - - std::size_t turn_index; - std::size_t operation_index; - bool discarded; - // use pointers to avoid copies, const& is not possible because of usage in vector - segment_identifier const* other_seg_id; // segment id of other segment of intersection of two segments - TurnOperation const* subject; - - inline indexed_turn_operation(std::size_t ti, std::size_t oi, - TurnOperation const& sub, - segment_identifier const& oid) - : turn_index(ti) - , operation_index(oi) - , discarded(false) - , other_seg_id(&oid) - , subject(boost::addressof(sub)) - {} - -}; - -template <typename IndexedTurnOperation> -struct remove_discarded -{ - inline bool operator()(IndexedTurnOperation const& operation) const - { - return operation.discarded; - } -}; - - -template -< - typename TurnPoints, - typename Indexed, - typename Geometry1, typename Geometry2, - typename RobustPolicy, - bool Reverse1, bool Reverse2, - typename Strategy -> -struct sort_on_segment_and_ratio -{ - inline sort_on_segment_and_ratio(TurnPoints const& turn_points - , Geometry1 const& geometry1 - , Geometry2 const& geometry2 - , RobustPolicy const& robust_policy - , Strategy const& strategy - , bool* clustered) - : m_turn_points(turn_points) - , m_geometry1(geometry1) - , m_geometry2(geometry2) - , m_robust_policy(robust_policy) - , m_strategy(strategy) - , m_clustered(clustered) - { - } - -private : - - TurnPoints const& m_turn_points; - Geometry1 const& m_geometry1; - Geometry2 const& m_geometry2; - RobustPolicy const& m_robust_policy; - Strategy const& m_strategy; - mutable bool* m_clustered; - - typedef typename geometry::point_type<Geometry1>::type point_type; - - inline bool default_order(Indexed const& left, Indexed const& right) const - { - // We've nothing to sort on. Take the indexes - return left.turn_index < right.turn_index; - } - - inline bool consider_relative_order(Indexed const& left, - Indexed const& right) const - { - point_type pi, pj, ri, rj, si, sj; - - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - left.subject->seg_id, - pi, pj); - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - *left.other_seg_id, - ri, rj); - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - *right.other_seg_id, - si, sj); - - typedef typename strategy::side::services::default_strategy - < - typename cs_tag<point_type>::type - >::type strategy; - - int const side_rj_p = strategy::apply(pi, pj, rj); - int const side_sj_p = strategy::apply(pi, pj, sj); - - // Put the one turning left (1; right == -1) as last - if (side_rj_p != side_sj_p) - { - return side_rj_p < side_sj_p; - } - - int const side_sj_r = strategy::apply(ri, rj, sj); - int const side_rj_s = strategy::apply(si, sj, rj); - - // If they both turn left: the most left as last - // If they both turn right: this is not relevant, but take also here most left - if (side_rj_s != side_sj_r) - { - return side_rj_s < side_sj_r; - } - - return default_order(left, right); - } - -public : - - // Note that left/right do NOT correspond to m_geometry1/m_geometry2 - // but to the "indexed_turn_operation" - inline bool operator()(Indexed const& left, Indexed const& right) const - { - if (! (left.subject->seg_id == right.subject->seg_id)) - { - return left.subject->seg_id < right.subject->seg_id; - } - - // Both left and right are located on the SAME segment. - - if (! (left.subject->fraction == right.subject->fraction)) - { - return left.subject->fraction < right.subject->fraction; - } - - - // First check "real" intersection (crosses) - // -> distance zero due to precision, solve it by sorting - if (m_turn_points[left.turn_index].method == method_crosses - && m_turn_points[right.turn_index].method == method_crosses) - { - return consider_relative_order(left, right); - } - - // If that is not the case, cluster it later on. - // Indicate that this is necessary. - *m_clustered = true; - - return default_order(left, right); - } -}; - - -template<typename Turns, typename Operations> -inline void update_discarded(Turns& turn_points, Operations& operations) -{ - // Vice-versa, set discarded to true for discarded operations; - // AND set discarded points to true - for (typename boost::range_iterator<Operations>::type it = boost::begin(operations); - it != boost::end(operations); - ++it) - { - if (turn_points[it->turn_index].discarded) - { - it->discarded = true; - } - else if (it->discarded) - { - turn_points[it->turn_index].discarded = true; - } - } -} - // Sorts IP-s of this ring on segment-identifier, and if on same segment, // on distance. @@ -227,148 +53,88 @@ inline void update_discarded(Turns& turn_points, Operations& operations) // (might be on another segment) template < - typename IndexType, bool Reverse1, bool Reverse2, - typename Container, - typename TurnPoints, + typename Operations, + typename Turns, typename Geometry1, typename Geometry2, typename RobustPolicy, typename Strategy > -inline void enrich_sort(Container& operations, - TurnPoints& turn_points, +inline void enrich_sort(Operations& operations, + Turns const& turns, operation_type for_operation, - Geometry1 const& geometry1, Geometry2 const& geometry2, + Geometry1 const& geometry1, + Geometry2 const& geometry2, RobustPolicy const& robust_policy, Strategy const& strategy) { - typedef typename IndexType::type operations_type; - - bool clustered = false; std::sort(boost::begin(operations), - boost::end(operations), - sort_on_segment_and_ratio - < - TurnPoints, - IndexType, - Geometry1, Geometry2, - RobustPolicy, - Reverse1, Reverse2, - Strategy - >(turn_points, geometry1, geometry2, robust_policy, strategy, &clustered)); - - // DONT'T discard xx / (for union) ix / ii / (for intersection) ux / uu here - // It would give way to "lonely" ui turn points, traveling all - // the way round. See #105 - - if (clustered) - { - typedef typename boost::range_iterator<Container>::type nc_iterator; - nc_iterator it = boost::begin(operations); - nc_iterator begin_cluster = boost::end(operations); - for (nc_iterator prev = it++; - it != boost::end(operations); - prev = it++) - { - operations_type& prev_op = turn_points[prev->turn_index] - .operations[prev->operation_index]; - operations_type& op = turn_points[it->turn_index] - .operations[it->operation_index]; - - if (prev_op.seg_id == op.seg_id - && (turn_points[prev->turn_index].method != method_crosses - || turn_points[it->turn_index].method != method_crosses) - && prev_op.fraction == op.fraction) - { - if (begin_cluster == boost::end(operations)) - { - begin_cluster = prev; - } - } - else if (begin_cluster != boost::end(operations)) - { - handle_cluster<IndexType, Reverse1, Reverse2>(begin_cluster, it, turn_points, - for_operation, geometry1, geometry2, robust_policy, strategy); - begin_cluster = boost::end(operations); - } - } - if (begin_cluster != boost::end(operations)) - { - handle_cluster<IndexType, Reverse1, Reverse2>(begin_cluster, it, turn_points, - for_operation, geometry1, geometry2, robust_policy, strategy); - } - } - - update_discarded(turn_points, operations); + boost::end(operations), + less_by_segment_ratio + < + Turns, + typename boost::range_value<Operations>::type, + Geometry1, Geometry2, + RobustPolicy, + Reverse1, Reverse2 + >(turns, for_operation, geometry1, geometry2, robust_policy)); } -template -< - typename IndexType, - typename Container, - typename TurnPoints -> -inline void enrich_discard(Container& operations, TurnPoints& turn_points) -{ - update_discarded(turn_points, operations); - - // Then delete discarded operations from vector - remove_discarded<IndexType> predicate; - operations.erase( - std::remove_if(boost::begin(operations), - boost::end(operations), - predicate), - boost::end(operations)); -} - -template -< - typename IndexType, - typename Container, - typename TurnPoints, - typename Geometry1, - typename Geometry2, - typename Strategy -> -inline void enrich_assign(Container& operations, - TurnPoints& turn_points, - operation_type , - Geometry1 const& , Geometry2 const& , - Strategy const& ) +template <typename Operations, typename Turns> +inline void enrich_assign(Operations& operations, Turns& turns) { - typedef typename IndexType::type operations_type; - typedef typename boost::range_iterator<Container const>::type iterator_type; + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::turn_operation_type op_type; + typedef typename boost::range_iterator<Operations>::type iterator_type; if (operations.size() > 0) { // Assign travel-to-vertex/ip index for each turning point. - // Because IP's are circular, PREV starts at the very last one, - // being assigned from the first one. - // "next ip on same segment" should not be considered circular. - bool first = true; - iterator_type it = boost::begin(operations); - for (iterator_type prev = it + (boost::size(operations) - 1); - it != boost::end(operations); - prev = it++) + // Iterator "next" is circular + + geometry::ever_circling_range_iterator<Operations const> next(operations); + ++next; + + for (iterator_type it = boost::begin(operations); + it != boost::end(operations); ++it) { - operations_type& prev_op - = turn_points[prev->turn_index].operations[prev->operation_index]; - operations_type& op - = turn_points[it->turn_index].operations[it->operation_index]; - - prev_op.enriched.travels_to_ip_index - = static_cast<signed_size_type>(it->turn_index); - prev_op.enriched.travels_to_vertex_index - = it->subject->seg_id.segment_index; - - if (! first - && prev_op.seg_id.segment_index == op.seg_id.segment_index) + turn_type& turn = turns[it->turn_index]; + op_type& op = turn.operations[it->operation_index]; + + // Normal behaviour: next should point at next turn: + if (it->turn_index == next->turn_index) + { + ++next; + } + + // Cluster behaviour: next should point after cluster, unless + // their seg_ids are not the same + while (turn.cluster_id != -1 + && it->turn_index != next->turn_index + && turn.cluster_id == turns[next->turn_index].cluster_id + && op.seg_id == turns[next->turn_index].operations[next->operation_index].seg_id) + { + ++next; + } + + turn_type const& next_turn = turns[next->turn_index]; + op_type const& next_op = next_turn.operations[next->operation_index]; + + op.enriched.travels_to_ip_index + = static_cast<signed_size_type>(next->turn_index); + op.enriched.travels_to_vertex_index + = next->subject->seg_id.segment_index; + + if (op.seg_id.segment_index == next_op.seg_id.segment_index + && op.fraction < next_op.fraction) { - prev_op.enriched.next_ip_index = static_cast<signed_size_type>(it->turn_index); + // Next turn is located further on same segment + // assign next_ip_index + // (this is one not circular therefore fraction is considered) + op.enriched.next_ip_index = static_cast<signed_size_type>(next->turn_index); } - first = false; } } @@ -379,20 +145,20 @@ inline void enrich_assign(Container& operations, it != boost::end(operations); ++it) { - operations_type& op = turn_points[it->turn_index] + op_type& op = turns[it->turn_index] .operations[it->operation_index]; std::cout << it->turn_index - << " meth: " << method_char(turn_points[it->turn_index].method) - << " seg: " << op.seg_id - << " dst: " << op.fraction // needs define - << " op: " << operation_char(turn_points[it->turn_index].operations[0].operation) - << operation_char(turn_points[it->turn_index].operations[1].operation) - << " dsc: " << (turn_points[it->turn_index].discarded ? "T" : "F") - << " ->vtx " << op.enriched.travels_to_vertex_index - << " ->ip " << op.enriched.travels_to_ip_index - << " ->nxt ip " << op.enriched.next_ip_index - //<< " vis: " << visited_char(op.visited) + << " cl=" << turns[it->turn_index].cluster_id + << " meth=" << method_char(turns[it->turn_index].method) + << " seg=" << op.seg_id + << " dst=" << op.fraction // needs define + << " op=" << operation_char(turns[it->turn_index].operations[0].operation) + << operation_char(turns[it->turn_index].operations[1].operation) + << " (" << operation_char(op.operation) << ")" + << " nxt=" << op.enriched.next_ip_index + << " / " << op.enriched.travels_to_ip_index + << " [vx " << op.enriched.travels_to_vertex_index << "]" << std::endl; ; } @@ -403,43 +169,57 @@ inline void enrich_assign(Container& operations, } -template <typename IndexedType, typename TurnPoints, typename MappedVector> -inline void create_map(TurnPoints const& turn_points, MappedVector& mapped_vector) +template <typename Turns, typename MappedVector> +inline void create_map(Turns const& turns, + detail::overlay::operation_type for_operation, + MappedVector& mapped_vector) { - typedef typename boost::range_value<TurnPoints>::type turn_point_type; - typedef typename turn_point_type::container_type container_type; + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::container_type container_type; + typedef typename MappedVector::mapped_type mapped_type; + typedef typename boost::range_value<mapped_type>::type indexed_type; std::size_t index = 0; - for (typename boost::range_iterator<TurnPoints const>::type - it = boost::begin(turn_points); - it != boost::end(turn_points); + for (typename boost::range_iterator<Turns const>::type + it = boost::begin(turns); + it != boost::end(turns); ++it, ++index) { - // Add operations on this ring, but skip discarded ones - if (! it->discarded) + // Add all (non discarded) operations on this ring + // Blocked operations or uu on clusters (for intersection) + // should be included, to block potential paths in clusters + turn_type const& turn = *it; + if (turn.discarded) { - std::size_t op_index = 0; - for (typename boost::range_iterator<container_type const>::type - op_it = boost::begin(it->operations); - op_it != boost::end(it->operations); - ++op_it, ++op_index) - { - // We should NOT skip blocked operations here - // because they can be relevant for "the other side" - // NOT if (op_it->operation != operation_blocked) - - ring_identifier ring_id - ( - op_it->seg_id.source_index, - op_it->seg_id.multi_index, - op_it->seg_id.ring_index - ); - mapped_vector[ring_id].push_back - ( - IndexedType(index, op_index, *op_it, - it->operations[1 - op_index].seg_id) - ); - } + continue; + } + + if (for_operation == operation_intersection + && turn.cluster_id == -1 + && turn.both(operation_union)) + { + // Only include uu turns if part of cluster (to block potential paths), + // otherwise they can block possibly viable paths + continue; + } + + std::size_t op_index = 0; + for (typename boost::range_iterator<container_type const>::type + op_it = boost::begin(turn.operations); + op_it != boost::end(turn.operations); + ++op_it, ++op_index) + { + ring_identifier const ring_id + ( + op_it->seg_id.source_index, + op_it->seg_id.multi_index, + op_it->seg_id.ring_index + ); + mapped_vector[ring_id].push_back + ( + indexed_type(index, op_index, *op_it, + it->operations[1 - op_index].seg_id) + ); } } } @@ -453,12 +233,12 @@ inline void create_map(TurnPoints const& turn_points, MappedVector& mapped_vecto /*! \brief All intersection points are enriched with successor information \ingroup overlay -\tparam TurnPoints type of intersection container +\tparam Turns type of intersection container (e.g. vector of "intersection/turn point"'s) \tparam Geometry1 \tparam_geometry \tparam Geometry2 \tparam_geometry \tparam Strategy side strategy type -\param turn_points container containing intersectionpoints +\param turns container containing intersectionpoints \param for_operation operation_type (union or intersection) \param geometry1 \param_geometry \param geometry2 \param_geometry @@ -469,22 +249,24 @@ template < bool Reverse1, bool Reverse2, overlay_type OverlayType, - typename TurnPoints, + typename Turns, + typename Clusters, typename Geometry1, typename Geometry2, typename RobustPolicy, typename Strategy > -inline void enrich_intersection_points(TurnPoints& turn_points, +inline void enrich_intersection_points(Turns& turns, + Clusters& clusters, detail::overlay::operation_type for_operation, Geometry1 const& geometry1, Geometry2 const& geometry2, RobustPolicy const& robust_policy, Strategy const& strategy) { - typedef typename boost::range_value<TurnPoints>::type turn_point_type; - typedef typename turn_point_type::turn_operation_type turn_operation_type; + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::turn_operation_type op_type; typedef detail::overlay::indexed_turn_operation < - turn_operation_type + op_type > indexed_turn_operation; typedef std::map @@ -493,51 +275,27 @@ inline void enrich_intersection_points(TurnPoints& turn_points, std::vector<indexed_turn_operation> > mapped_vector_type; - // Iterate through turns and discard uu - // and check if there are possible colocations - bool check_colocations = false; - for (typename boost::range_iterator<TurnPoints>::type - it = boost::begin(turn_points); - it != boost::end(turn_points); + bool const has_colocations + = detail::overlay::handle_colocations<Reverse1, Reverse2>(turns, + clusters, geometry1, geometry2); + + // Discard none turns, if any + for (typename boost::range_iterator<Turns>::type + it = boost::begin(turns); + it != boost::end(turns); ++it) { - if (it->both(detail::overlay::operation_union)) - { - // Discard (necessary for a.o. #76). With uu, at all points there - // is the risk that rings are being traversed twice or more. - // Without uu, all rings having only uu will be untouched - // and gathered by assemble - it->discarded = true; - check_colocations = true; - } - else if (it->combination(detail::overlay::operation_union, - detail::overlay::operation_blocked)) - { - check_colocations = true; - } - else if (OverlayType == overlay_difference - && it->both(detail::overlay::operation_intersection)) - { - // For difference operation (u/u -> i/i) - check_colocations = true; - } - else if (it->both(detail::overlay::operation_none)) + if (it->both(detail::overlay::operation_none)) { it->discarded = true; } } - if (check_colocations) - { - detail::overlay::handle_colocations<OverlayType>(turn_points); - } - // Create a map of vectors of indexed operation-types to be able // to sort intersection points PER RING mapped_vector_type mapped_vector; - detail::overlay::create_map<indexed_turn_operation>(turn_points, mapped_vector); - + detail::overlay::create_map(turns, for_operation, mapped_vector); // No const-iterator; contents of mapped copy is temporary, // and changed by enrich @@ -550,8 +308,10 @@ inline void enrich_intersection_points(TurnPoints& turn_points, std::cout << "ENRICH-sort Ring " << mit->first << std::endl; #endif - detail::overlay::enrich_sort<indexed_turn_operation, Reverse1, Reverse2>(mit->second, turn_points, for_operation, - geometry1, geometry2, robust_policy, strategy); + detail::overlay::enrich_sort<Reverse1, Reverse2>( + mit->second, turns, for_operation, + geometry1, geometry2, + robust_policy, strategy); } for (typename mapped_vector_type::iterator mit @@ -560,27 +320,20 @@ inline void enrich_intersection_points(TurnPoints& turn_points, ++mit) { #ifdef BOOST_GEOMETRY_DEBUG_ENRICH - std::cout << "ENRICH-discard Ring " + std::cout << "ENRICH-assign Ring " << mit->first << std::endl; #endif - detail::overlay::enrich_discard<indexed_turn_operation>(mit->second, turn_points); + detail::overlay::enrich_assign(mit->second, turns); } - for (typename mapped_vector_type::iterator mit - = mapped_vector.begin(); - mit != mapped_vector.end(); - ++mit) + if (has_colocations) { -#ifdef BOOST_GEOMETRY_DEBUG_ENRICH - std::cout << "ENRICH-assign Ring " - << mit->first << std::endl; -#endif - detail::overlay::enrich_assign<indexed_turn_operation>(mit->second, turn_points, for_operation, - geometry1, geometry2, strategy); + detail::overlay::assign_startable_in_clusters<Reverse1, Reverse2>( + clusters, turns, for_operation, geometry1, geometry2); } #ifdef BOOST_GEOMETRY_DEBUG_ENRICH - //detail::overlay::check_graph(turn_points, for_operation); + //detail::overlay::check_graph(turns, for_operation); #endif } diff --git a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp index 5964f3ba44..f5c5cc6a2e 100644 --- a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp @@ -32,6 +32,9 @@ struct enrichment_info : travels_to_vertex_index(-1) , travels_to_ip_index(-1) , next_ip_index(-1) + , startable(true) + , count_left(0) + , count_right(0) {} // vertex to which is free travel after this IP, @@ -44,6 +47,12 @@ struct enrichment_info // index of next IP on this segment, -1 if there is no one signed_size_type next_ip_index; + + bool startable; // Can be used to start in traverse + + // Counts if polygons left/right of this operation + std::size_t count_left; + std::size_t count_right; }; diff --git a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp index 6f332ddff2..ad75127fab 100644 --- a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp +++ b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp @@ -16,6 +16,7 @@ #include <boost/range.hpp> #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp> +#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> #include <boost/geometry/algorithms/detail/ring_identifier.hpp> #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp> @@ -34,6 +35,29 @@ namespace boost { namespace geometry namespace detail { namespace overlay { +template <typename SegmentRatio> +struct segment_fraction +{ + segment_identifier seg_id; + SegmentRatio fraction; + + segment_fraction(segment_identifier const& id, SegmentRatio const& fr) + : seg_id(id) + , fraction(fr) + {} + + segment_fraction() + {} + + bool operator<(segment_fraction<SegmentRatio> const& other) const + { + return seg_id == other.seg_id + ? fraction < other.fraction + : seg_id < other.seg_id; + } + +}; + struct turn_operation_index { turn_operation_index(signed_size_type ti = -1, @@ -43,22 +67,22 @@ struct turn_operation_index {} signed_size_type turn_index; - signed_size_type op_index; // basically only 0,1 + signed_size_type op_index; // only 0,1 }; -template <typename TurnPoints> +template <typename Turns> struct less_by_fraction_and_type { - inline less_by_fraction_and_type(TurnPoints const& turn_points) - : m_turns(turn_points) + inline less_by_fraction_and_type(Turns const& turns) + : m_turns(turns) { } inline bool operator()(turn_operation_index const& left, turn_operation_index const& right) const { - typedef typename boost::range_value<TurnPoints>::type turn_type; + typedef typename boost::range_value<Turns>::type turn_type; typedef typename turn_type::turn_operation_type turn_operation_type; turn_type const& left_turn = m_turns[left.turn_index]; @@ -74,6 +98,29 @@ struct less_by_fraction_and_type return left_op.fraction < right_op.fraction; } + // Order xx first - used to discard any following colocated turn + bool const left_both_xx = left_turn.both(operation_blocked); + bool const right_both_xx = right_turn.both(operation_blocked); + if (left_both_xx && ! right_both_xx) + { + return true; + } + if (! left_both_xx && right_both_xx) + { + return false; + } + + bool const left_both_uu = left_turn.both(operation_union); + bool const right_both_uu = right_turn.both(operation_union); + if (left_both_uu && ! right_both_uu) + { + return true; + } + if (! left_both_uu && right_both_uu) + { + return false; + } + turn_operation_type const& left_other_op = left_turn.operations[1 - left.op_index]; @@ -86,88 +133,219 @@ struct less_by_fraction_and_type } private: - TurnPoints const& m_turns; + Turns const& m_turns; }; -template <overlay_type OverlayType, typename TurnPoints, typename OperationVector> -inline void handle_colocation_cluster(TurnPoints& turn_points, - segment_identifier const& current_ring_seg_id, - OperationVector const& vec) +template <typename Operation, typename ClusterPerSegment> +inline signed_size_type get_cluster_id(Operation const& op, ClusterPerSegment const& cluster_per_segment) +{ + typedef typename ClusterPerSegment::key_type segment_fraction_type; + + segment_fraction_type seg_frac(op.seg_id, op.fraction); + typename ClusterPerSegment::const_iterator it + = cluster_per_segment.find(seg_frac); + + if (it == cluster_per_segment.end()) + { + return -1; + } + return it->second; +} + +template <typename Operation, typename ClusterPerSegment> +inline void add_cluster_id(Operation const& op, + ClusterPerSegment& cluster_per_segment, signed_size_type id) { - typedef typename boost::range_value<TurnPoints>::type turn_type; + typedef typename ClusterPerSegment::key_type segment_fraction_type; + + segment_fraction_type seg_frac(op.seg_id, op.fraction); + + cluster_per_segment[seg_frac] = id; +} + +template <typename Turn, typename ClusterPerSegment> +inline signed_size_type add_turn_to_cluster(Turn const& turn, + ClusterPerSegment& cluster_per_segment, signed_size_type& cluster_id) +{ + signed_size_type cid0 = get_cluster_id(turn.operations[0], cluster_per_segment); + signed_size_type cid1 = get_cluster_id(turn.operations[1], cluster_per_segment); + + if (cid0 == -1 && cid1 == -1) + { + ++cluster_id; + add_cluster_id(turn.operations[0], cluster_per_segment, cluster_id); + add_cluster_id(turn.operations[1], cluster_per_segment, cluster_id); + return cluster_id; + } + else if (cid0 == -1 && cid1 != -1) + { + add_cluster_id(turn.operations[0], cluster_per_segment, cid1); + return cid1; + } + else if (cid0 != -1 && cid1 == -1) + { + add_cluster_id(turn.operations[1], cluster_per_segment, cid0); + return cid0; + } + else if (cid0 == cid1) + { + // Both already added to same cluster, no action + return cid0; + } + + // Both operations.seg_id/fraction were already part of any cluster, and + // these clusters are not the same. Merge of two clusters is necessary + std::cout << " TODO: merge " << cid0 << " and " << cid1 << std::endl; + return cid0; +} + +template +< + bool Reverse1, bool Reverse2, + typename Turns, + typename ClusterPerSegment, + typename Operations, + typename Geometry1, + typename Geometry2 +> +inline void handle_colocation_cluster(Turns& turns, + signed_size_type& cluster_id, + ClusterPerSegment& cluster_per_segment, + Operations const& operations, + Geometry1 const& geometry1, Geometry2 const& geometry2) +{ + typedef typename boost::range_value<Turns>::type turn_type; typedef typename turn_type::turn_operation_type turn_operation_type; - std::vector<turn_operation_index>::const_iterator vit = vec.begin(); + std::vector<turn_operation_index>::const_iterator vit = operations.begin(); - turn_type cluster_turn = turn_points[vit->turn_index]; - turn_operation_type cluster_op - = cluster_turn.operations[vit->op_index]; - segment_identifier cluster_other_id - = cluster_turn.operations[1 - vit->op_index].seg_id; - bool const discard_colocated - = cluster_turn.both(operation_union) - || cluster_turn.combination(operation_blocked, operation_union); + turn_operation_index ref_toi = *vit; + signed_size_type ref_id = -1; - for (++vit; vit != vec.end(); ++vit) + for (++vit; vit != operations.end(); ++vit) { + turn_type& ref_turn = turns[ref_toi.turn_index]; + turn_operation_type const& ref_op + = ref_turn.operations[ref_toi.op_index]; + turn_operation_index const& toi = *vit; - turn_type& turn = turn_points[toi.turn_index]; + turn_type& turn = turns[toi.turn_index]; turn_operation_type const& op = turn.operations[toi.op_index]; - segment_identifier const& other_id - = turn.operations[1 - toi.op_index].seg_id; - if (cluster_op.fraction == op.fraction) + BOOST_ASSERT(ref_op.seg_id == op.seg_id); + + if (ref_op.fraction == op.fraction) { - // Two turns of current ring with same source are colocated, - // one is from exterior ring, one from interior ring - bool const colocated_ext_int - = cluster_other_id.multi_index == other_id.multi_index - && cluster_other_id.ring_index == -1 - && other_id.ring_index >= 0; - - // Turn of current interior ring with other interior ring - bool const touch_int_int - = current_ring_seg_id.ring_index >= 0 - && other_id.ring_index >= 0; - - if (discard_colocated && colocated_ext_int) + turn_operation_type const& other_op = turn.operations[1 - toi.op_index]; + + if (ref_id == -1) { - // If the two turns on this same segment are a - // colocation with two different segments on the - // other geometry, of the same polygon but with - // the outer (u/u or u/x) and the inner ring (non u/u), - // that turn with inner ring should be discarded - turn.discarded = true; - turn.colocated = true; + ref_id = add_turn_to_cluster(ref_turn, cluster_per_segment, cluster_id); } - else if (cluster_turn.colocated - && touch_int_int - && turn.both(operation_intersection)) + BOOST_ASSERT(ref_id != -1); + + // ref_turn (both operations) are already added to cluster, + // so also "op" is already added to cluster, + // We only need to add other_op + signed_size_type id = get_cluster_id(other_op, cluster_per_segment); + if (id != -1 && id != ref_id) { - // Two holes touch each other at a point where the - // exterior ring also touches - turn.discarded = true; - turn.colocated = true; } - else if (OverlayType == overlay_difference - && turn.both(operation_intersection) - && colocated_ext_int) + else if (id == -1) + { + // Add to same cluster + add_cluster_id(other_op, cluster_per_segment, ref_id); + id = ref_id; + } + + // In case of colocated xx turns, all other turns may NOT be + // followed at all. xx cannot be discarded (otherwise colocated + // turns are followed). + if (ref_turn.both(operation_blocked)) { - // For difference (polygon inside out) we need to - // discard i/i instead, in case of colocations turn.discarded = true; - turn.colocated = true; + // We can either set or not set colocated because it is not effective on blocked turns } } else { // Not on same fraction on this segment - // assign for next potential cluster - cluster_turn = turn; - cluster_op = op; - cluster_other_id = other_id; + // assign for next + ref_toi = toi; + ref_id = -1; } + } +} + +template +< + typename Turns, + typename Clusters, + typename ClusterPerSegment +> +inline void assign_cluster_to_turns(Turns& turns, + Clusters& clusters, + ClusterPerSegment const& cluster_per_segment) +{ + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::turn_operation_type turn_operation_type; + typedef typename ClusterPerSegment::key_type segment_fraction_type; + signed_size_type turn_index = 0; + for (typename boost::range_iterator<Turns>::type it = turns.begin(); + it != turns.end(); ++it, ++turn_index) + { + turn_type& turn = *it; + + if (turn.discarded) + { + // They were processed (to create proper map) but will not be added + // This might leave a cluster with only 1 turn, which will be fixed + // afterwards + continue; + } + + for (int i = 0; i < 2; i++) + { + turn_operation_type const& op = turn.operations[i]; + segment_fraction_type seg_frac(op.seg_id, op.fraction); + typename ClusterPerSegment::const_iterator it = cluster_per_segment.find(seg_frac); + if (it != cluster_per_segment.end()) + { + if (turn.cluster_id != -1 + && turn.cluster_id != it->second) + { + std::cout << " CONFLICT " << std::endl; + } + turn.cluster_id = it->second; + clusters[turn.cluster_id].insert(turn_index); + } + } + } +} + +template +< + typename Turns, + typename Clusters +> +inline void remove_clusters(Turns& turns, Clusters& clusters) +{ + typename Clusters::iterator it = clusters.begin(); + while (it != clusters.end()) + { + // Hold iterator and increase. We can erase cit, this keeps the + // iterator valid (cf The standard associative-container erase idiom) + typename Clusters::iterator current_it = it; + ++it; + + std::set<signed_size_type> const& turn_indices = current_it->second; + if (turn_indices.size() == 1) + { + signed_size_type turn_index = *turn_indices.begin(); + turns[turn_index].cluster_id = -1; + clusters.erase(current_it); + } } } @@ -179,8 +357,16 @@ inline void handle_colocation_cluster(TurnPoints& turn_points, // This function can be extended to replace handle_tangencies: at each // colocation incoming and outgoing vectors should be inspected -template <overlay_type OverlayType, typename TurnPoints> -inline void handle_colocations(TurnPoints& turn_points) +template +< + bool Reverse1, bool Reverse2, + typename Turns, + typename Clusters, + typename Geometry1, + typename Geometry2 +> +inline bool handle_colocations(Turns& turns, Clusters& clusters, + Geometry1 const& geometry1, Geometry2 const& geometry2) { typedef std::map < @@ -195,9 +381,9 @@ inline void handle_colocations(TurnPoints& turn_points) map_type map; int index = 0; - for (typename boost::range_iterator<TurnPoints>::type - it = boost::begin(turn_points); - it != boost::end(turn_points); + for (typename boost::range_iterator<Turns>::type + it = boost::begin(turns); + it != boost::end(turns); ++it, ++index) { map[it->operations[0].seg_id].push_back(turn_operation_index(index, 0)); @@ -220,27 +406,43 @@ inline void handle_colocations(TurnPoints& turn_points) if (! colocations) { - return; + return false; } // Sort all vectors, per same segment - less_by_fraction_and_type<TurnPoints> less(turn_points); + less_by_fraction_and_type<Turns> less(turns); for (typename map_type::iterator it = map.begin(); it != map.end(); ++it) { std::sort(it->second.begin(), it->second.end(), less); } + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::segment_ratio_type segment_ratio_type; + + typedef std::map + < + segment_fraction<segment_ratio_type>, + signed_size_type + > cluster_per_segment_type; + + cluster_per_segment_type cluster_per_segment; + signed_size_type cluster_id = 0; + for (typename map_type::const_iterator it = map.begin(); it != map.end(); ++it) { - if (it->second.size() > 1) + if (it->second.size() > 1u) { - handle_colocation_cluster<OverlayType>(turn_points, - it->first, it->second); + handle_colocation_cluster<Reverse1, Reverse2>(turns, cluster_id, + cluster_per_segment, it->second, + geometry1, geometry2); } } + assign_cluster_to_turns(turns, clusters, cluster_per_segment); + remove_clusters(turns, clusters); + #if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS) std::cout << "*** Colocations " << map.size() << std::endl; for (typename map_type::const_iterator it = map.begin(); @@ -251,21 +453,123 @@ inline void handle_colocations(TurnPoints& turn_points) = it->second.begin(); vit != it->second.end(); ++vit) { turn_operation_index const& toi = *vit; - std::cout << geometry::wkt(turn_points[toi.turn_index].point) + std::cout << geometry::wkt(turns[toi.turn_index].point) << std::boolalpha - << " discarded=" << turn_points[toi.turn_index].discarded - << " colocated=" << turn_points[toi.turn_index].colocated - << " " << operation_char(turn_points[toi.turn_index].operations[0].operation) - << " " << turn_points[toi.turn_index].operations[0].seg_id - << " " << turn_points[toi.turn_index].operations[0].fraction - << " // " << operation_char(turn_points[toi.turn_index].operations[1].operation) - << " " << turn_points[toi.turn_index].operations[1].seg_id - << " " << turn_points[toi.turn_index].operations[1].fraction + << " discarded=" << turns[toi.turn_index].discarded + << " colocated=" << turns[toi.turn_index].colocated + << " " << operation_char(turns[toi.turn_index].operations[0].operation) + << " " << turns[toi.turn_index].operations[0].seg_id + << " " << turns[toi.turn_index].operations[0].fraction + << " // " << operation_char(turns[toi.turn_index].operations[1].operation) + << " " << turns[toi.turn_index].operations[1].seg_id + << " " << turns[toi.turn_index].operations[1].fraction << std::endl; } } #endif // DEBUG + return true; +} + + +struct is_turn_index +{ + is_turn_index(signed_size_type index) + : m_index(index) + {} + + template <typename Indexed> + inline bool operator()(Indexed const& indexed) const + { + // Indexed is a indexed_turn_operation<Operation> + return indexed.turn_index == m_index; + } + + std::size_t m_index; +}; + + +template +< + bool Reverse1, bool Reverse2, + typename Turns, + typename Clusters, + typename Geometry1, + typename Geometry2 +> +inline void assign_startable_in_clusters(Clusters& clusters, Turns& turns, + operation_type for_operation, + Geometry1 const& geometry1, Geometry2 const& geometry2) +{ + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::point_type point_type; + typedef typename turn_type::turn_operation_type turn_operation_type; + + // Define sorter, sorting counter-clockwise such that polygons are on the + // right side + typedef sort_by_side::side_sorter + < + Reverse1, Reverse2, point_type, std::less<int> + > sbs_type; + + for (typename Clusters::iterator mit = clusters.begin(); + mit != clusters.end(); ++mit) + { + std::set<signed_size_type> const& ids = mit->second; + if (ids.empty()) + { + continue; + } + + sbs_type sbs; + point_type turn_point; // should be all the same for all turns in cluster + + bool first = true; + for (typename std::set<signed_size_type>::const_iterator sit = ids.begin(); + sit != ids.end(); ++sit) + { + signed_size_type turn_index = *sit; + turn_type const& turn = turns[turn_index]; + if (first) + { + turn_point = turn.point; + } + for (int i = 0; i < 2; i++) + { + turn_operation_type const& op = turn.operations[i]; + sbs.add(op, turn_index, i, geometry1, geometry2, first); + first = false; + } + } + sbs.apply(turn_point); + + sbs.find_open(); + + // Unset the startable flag for all 'closed' spaces + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + const typename sbs_type::rp& ranked = sbs.m_ranked_points[i]; + turn_type& turn = turns[ranked.turn_index]; + turn_operation_type& op = turn.operations[ranked.op_index]; + + if (ranked.index != sort_by_side::index_to) + { + continue; + } + + op.enriched.count_left = ranked.left_count; + op.enriched.count_right = ranked.right_count; + + if ((for_operation == operation_union + && ranked.left_count != 0) + || (for_operation == operation_intersection + && ranked.right_count != 2)) + { + op.enriched.startable = false; + } + } + + } } diff --git a/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp b/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp deleted file mode 100644 index 277a223d47..0000000000 --- a/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp +++ /dev/null @@ -1,786 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) - -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. - -// 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) - -#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TANGENCIES_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TANGENCIES_HPP - -#include <algorithm> - -#include <boost/geometry/algorithms/detail/ring_identifier.hpp> -#include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp> -#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> -#include <boost/geometry/algorithms/detail/recalculate.hpp> - -#include <boost/geometry/policies/robustness/robust_point_type.hpp> -#include <boost/geometry/policies/robustness/segment_ratio_type.hpp> -#include <boost/geometry/policies/robustness/robust_type.hpp> - -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) -#include <boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp> -#endif - -#include <boost/geometry/geometries/point.hpp> -#include <boost/geometry/geometries/segment.hpp> - - -// TODO: the approach below should be completely replaced by the new -// get_left_turns, to keep the outgoing vector which has open space one of its -// sides. - - -namespace boost { namespace geometry -{ - -#ifndef DOXYGEN_NO_DETAIL -namespace detail { namespace overlay -{ - - -template -< - typename TurnPoints, - typename Indexed, - typename Geometry1, typename Geometry2, - typename RobustPolicy, - bool Reverse1, bool Reverse2, - typename Strategy -> -struct sort_in_cluster -{ - inline sort_in_cluster(TurnPoints const& turn_points - , Geometry1 const& geometry1 - , Geometry2 const& geometry2 - , RobustPolicy const& robust_policy - , Strategy const& strategy) - : m_turn_points(turn_points) - , m_geometry1(geometry1) - , m_geometry2(geometry2) - , m_rescale_policy(robust_policy) - , m_strategy(strategy) - {} - -private : - - TurnPoints const& m_turn_points; - Geometry1 const& m_geometry1; - Geometry2 const& m_geometry2; - RobustPolicy const& m_rescale_policy; - Strategy const& m_strategy; - - typedef typename Indexed::type turn_operation_type; - typedef typename geometry::point_type<Geometry1>::type point_type; - - typedef typename geometry::robust_point_type - < - point_type, - RobustPolicy - >::type robust_point_type; - - inline bool default_order(Indexed const& left, Indexed const& right) const - { - // We've nothing to sort on. Take the indexes - return left.turn_index < right.turn_index; - } - - // Still necessary in some situations, - // for example #case_102_multi, #case_107_multi, #case_recursive_boxes_3 - inline void get_situation_map(Indexed const& left, Indexed const& right, - robust_point_type& pi_rob, robust_point_type& pj_rob, - robust_point_type& ri_rob, robust_point_type& rj_rob, - robust_point_type& si_rob, robust_point_type& sj_rob) const - { - point_type pi, pj, ri, rj, si, sj; - - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - left.subject->seg_id, - pi, pj); - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - *left.other_seg_id, - ri, rj); - geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, - *right.other_seg_id, - si, sj); - - geometry::recalculate(pi_rob, pi, m_rescale_policy); - geometry::recalculate(pj_rob, pj, m_rescale_policy); - geometry::recalculate(ri_rob, ri, m_rescale_policy); - geometry::recalculate(rj_rob, rj, m_rescale_policy); - geometry::recalculate(si_rob, si, m_rescale_policy); - geometry::recalculate(sj_rob, sj, m_rescale_policy); - } - -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - // This method was still called but did no effect whatsoever on the results, - // with or without robustness-rescaling. - // Probable cause: we rescale in this file ourselves, ignoring passed policy - // TODO: check this more. - // Besides this, it currently does not compile for yet unknown reasons - // (does not find specialization for segment_ratio_type) - // It is currently only called from the Unit Test "multi_intersection.cpp" - - // Determine how p/r and p/s are located. - inline void overlap_info( - robust_point_type const& pi, robust_point_type const& pj, - robust_point_type const& ri, robust_point_type const& rj, - robust_point_type const& si, robust_point_type const& sj, - bool& pr_overlap, bool& ps_overlap, bool& rs_overlap) const - { - // Determine how p/r and p/s are located. - // One of them is coming from opposite direction. - - typedef segment_intersection_points - < - point_type, - typename segment_ratio_type - < - point_type, RobustPolicy - >::type - > intersection_return_type; - - typedef strategy::intersection::relate_cartesian_segments - < - policies::relate::segments_intersection_points - < - intersection_return_type - > - > policy; - - typedef model::referring_segment<robust_point_type const> segment_type; - segment_type p(pi, pj); - segment_type r(ri, rj); - segment_type s(si, sj); - - // Get the intersection point (or two points) - intersection_return_type pr = policy::apply(p, r, m_rescale_policy, pi, pj, ri, rj); - intersection_return_type ps = policy::apply(p, s, m_rescale_policy, pi, pj, si, sj); - intersection_return_type rs = policy::apply(r, s, m_rescale_policy, ri, rj, si, sj); - - // Check on overlap - pr_overlap = pr.count == 2; - ps_overlap = ps.count == 2; - rs_overlap = rs.count == 2; - } -#endif - - -#ifdef BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES - inline void debug_consider(int order, Indexed const& left, - Indexed const& right, std::string const& header, - bool skip = true, - std::string const& extra = "", bool ret = false - ) const - { - if (skip) return; - - std::cout << "Case: " << header << " for " << left.turn_index << " / " << right.turn_index << std::endl; - - robust_point_type pi, pj, ri, rj, si, sj; - get_situation_map(left, right, pi, pj, ri, rj, si, sj); - -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - bool prc = false, psc = false, rsc = false; - overlap_info(pi, pj, ri, rj, si, sj, prc, psc, rsc); -#endif - - int const side_ri_p = m_strategy.apply(pi, pj, ri); - int const side_rj_p = m_strategy.apply(pi, pj, rj); - int const side_si_p = m_strategy.apply(pi, pj, si); - int const side_sj_p = m_strategy.apply(pi, pj, sj); - int const side_si_r = m_strategy.apply(ri, rj, si); - int const side_sj_r = m_strategy.apply(ri, rj, sj); - -#ifdef BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES_MORE - std::cout << " Segment p:" << geometry::wkt(pi) << " .. " << geometry::wkt(pj) << std::endl; - std::cout << " Segment r:" << geometry::wkt(ri) << " .. " << geometry::wkt(rj) << std::endl; - std::cout << " Segment s:" << geometry::wkt(si) << " .. " << geometry::wkt(sj) << std::endl; - - std::cout << " r//p: " << side_ri_p << " / " << side_rj_p << std::endl; - std::cout << " s//p: " << side_si_p << " / " << side_sj_p << std::endl; - std::cout << " s//r: " << side_si_r << " / " << side_sj_r << std::endl; -#endif - - std::cout << header - //<< " order: " << order - << " ops: " << operation_char(left.subject->operation) - << "/" << operation_char(right.subject->operation) - << " ri//p: " << side_ri_p - << " si//p: " << side_si_p - << " si//r: " << side_si_r -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - << " cnts: " << int(prc) << "," << int(psc) << "," << int(rsc) -#endif - //<< " idx: " << left.turn_index << "/" << right.turn_index - ; - - if (! extra.empty()) - { - std::cout << " " << extra << " " << (ret ? "true" : "false"); - } - std::cout << std::endl; - } -#else - inline void debug_consider(int, Indexed const& , - Indexed const& , std::string const& , - bool = true, - std::string const& = "", bool = false - ) const - {} -#endif - - - // ux/ux - inline bool consider_ux_ux(Indexed const& left, - Indexed const& right - , std::string const& // header - ) const - { - bool ret = left.turn_index < right.turn_index; - - // In combination of u/x, x/u: take first union, then blocked. - // Solves #88, #61, #56, #80 - if (left.subject->operation == operation_union - && right.subject->operation == operation_blocked) - { - ret = true; - } - else if (left.subject->operation == operation_blocked - && right.subject->operation == operation_union) - { - ret = false; - } - else - { -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - std::cout << "ux/ux unhandled" << std::endl; -#endif - } - - //debug_consider(0, left, right, header, false, "-> return ", ret); - - return ret; - } - - inline bool consider_iu_ux(Indexed const& left, - Indexed const& right, - int order // 1: iu first, -1: ux first - , std::string const& // header - ) const - { - bool ret = false; - - if (left.subject->operation == operation_union - && right.subject->operation == operation_union) - { - ret = order == 1; - } - else if (left.subject->operation == operation_union - && right.subject->operation == operation_blocked) - { - ret = true; - } - else if (right.subject->operation == operation_union - && left.subject->operation == operation_blocked) - { - ret = false; - } - else if (left.subject->operation == operation_union) - { - ret = true; - } - else if (right.subject->operation == operation_union) - { - ret = false; - } - else - { -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - // this still happens in the traverse.cpp test - std::cout << " iu/ux unhandled" << std::endl; -#endif - ret = order == 1; - } - - //debug_consider(0, left, right, header, false, "-> return", ret); - return ret; - } - - inline bool consider_iu_ix(Indexed const& left, - Indexed const& right, - int order // 1: iu first, -1: ix first - , std::string const& // header - ) const - { - //debug_consider(order, left, right, header, false, "iu/ix"); - - return left.subject->operation == operation_intersection - && right.subject->operation == operation_intersection ? order == 1 - : left.subject->operation == operation_intersection ? false - : right.subject->operation == operation_intersection ? true - : order == 1; - } - - inline bool consider_ix_ix(Indexed const& left, Indexed const& right - , std::string const& // header - ) const - { - // Take first intersection, then blocked. - if (left.subject->operation == operation_intersection - && right.subject->operation == operation_blocked) - { - return true; - } - else if (left.subject->operation == operation_blocked - && right.subject->operation == operation_intersection) - { - return false; - } - - // Default case, should not occur - -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - std::cout << "ix/ix unhandled" << std::endl; -#endif - //debug_consider(0, left, right, header, false, "-> return", ret); - - return default_order(left, right); - } - - - inline bool consider_iu_iu(Indexed const& left, Indexed const& right, - std::string const& header, bool redo = false) const - { - //debug_consider(0, left, right, header); - - // In general, order it like "union, intersection". - if (left.subject->operation == operation_intersection - && right.subject->operation == operation_union) - { - //debug_consider(0, left, right, header, false, "i,u", false); - return false; - } - else if (left.subject->operation == operation_union - && right.subject->operation == operation_intersection) - { - //debug_consider(0, left, right, header, false, "u,i", true); - return true; - } - - robust_point_type pi, pj, ri, rj, si, sj; - get_situation_map(left, right, pi, pj, ri, rj, si, sj); - - int const side_ri_p = m_strategy.apply(pi, pj, ri); - int const side_si_p = m_strategy.apply(pi, pj, si); - int const side_si_r = m_strategy.apply(ri, rj, si); - - // Both located at same side (#58, pie_21_7_21_0_3) - if (side_ri_p * side_si_p == 1 && side_si_r != 0) - { - if (left.subject->operation == operation_union - && right.subject->operation == operation_union) - { - int const side_ri_s = m_strategy.apply(si, sj, ri); - if (side_si_r == side_ri_s) - { - return default_order(left, right); - } - - // Take the most left one - bool const ret = side_si_r == 1; - //debug_consider(0, left, right, header, false, "same side", ret); - return ret; - } - } - - - // Coming from opposite sides (#59, #99) - if (side_ri_p * side_si_p == -1) - { - bool ret = false; - - { - ret = side_ri_p == 1; // #100 - debug_consider(0, left, right, header, false, "opp.", ret); - return ret; - } -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - std::cout << " iu/iu coming from opposite unhandled" << std::endl; -#endif - } - -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - // We need EXTRA information here: are p/r/s overlapping? - bool pr_ov = false, ps_ov = false, rs_ov = false; - overlap_info(pi, pj, ri, rj, si, sj, pr_ov, ps_ov, rs_ov); -#else - // std::cout << "Boost.Geometry: skipping overlap_info" << std::endl; -#endif - - // One coming from right (#83,#90) - // One coming from left (#90, #94, #95) - if (side_si_r != 0 && (side_ri_p != 0 || side_si_p != 0)) - { - int const side_ri_s = m_strategy.apply(si, sj, ri); - if (side_si_r == side_ri_s) - { - return default_order(left, right); - } - - bool ret = false; - -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - if (pr_ov || ps_ov) - { - int r = side_ri_p != 0 ? side_ri_p : side_si_p; - ret = r * side_si_r == 1; - } - else -#endif - { - ret = side_si_r == 1; - } - - debug_consider(0, left, right, header, false, "left or right", ret); - return ret; - } - - // All aligned (#92, #96) - if (side_ri_p == 0 && side_si_p == 0 && side_si_r == 0) - { - // One of them is coming from opposite direction. - - // Take the one NOT overlapping - bool ret = false; - bool found = false; -#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO - if (pr_ov && ! ps_ov) - { - ret = true; - found = true; - } - else if (!pr_ov && ps_ov) - { - ret = false; - found = true; - } -#endif - - debug_consider(0, left, right, header, false, "aligned", ret); - if (found) - { - return ret; - } - } - -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - std::cout << " iu/iu unhandled" << std::endl; - debug_consider(0, left, right, header, false, "unhandled", left.turn_index < right.turn_index); -#endif - if (! redo) - { - // In some cases behaviour is not symmetrical. TODO: fix this properly - // OR: alternatively we might consider calling all these functions one-way anyway - return ! consider_iu_iu(right, left, header, true); - } - - return default_order(left, right); - } - - inline bool consider_ii(Indexed const& left, Indexed const& right, - std::string const& header) const - { - debug_consider(0, left, right, header); - - robust_point_type pi, pj, ri, rj, si, sj; - get_situation_map(left, right, pi, pj, ri, rj, si, sj); - - int const side_ri_p = m_strategy.apply(pi, pj, ri); - int const side_si_p = m_strategy.apply(pi, pj, si); - - // Two other points are (mostly) lying both right of the considered segment - // Take the most left one - int const side_si_r = m_strategy.apply(ri, rj, si); - if (side_ri_p == -1 - && side_si_p == -1 - && side_si_r != 0) - { - bool const ret = side_si_r != 1; - return ret; - } - return default_order(left, right); - } - - -public : - inline bool operator()(Indexed const& left, Indexed const& right) const - { - if ((m_turn_points[left.turn_index].discarded || left.discarded) - && (m_turn_points[right.turn_index].discarded || right.discarded)) - { - return default_order(left, right); - } - else if (m_turn_points[left.turn_index].discarded || left.discarded) - { - // Be careful to sort discarded first, then all others - return true; - } - else if (m_turn_points[right.turn_index].discarded || right.discarded) - { - // See above so return false here such that right (discarded) - // is sorted before left (not discarded) - return false; - } - else if (m_turn_points[left.turn_index].combination(operation_blocked, operation_union) - && m_turn_points[right.turn_index].combination(operation_blocked, operation_union)) - { - // ux/ux - return consider_ux_ux(left, right, "ux/ux"); - } - else if (m_turn_points[left.turn_index].both(operation_union) - && m_turn_points[right.turn_index].both(operation_union)) - { - // uu/uu, Order is arbitrary - // Note: uu/uu is discarded now before so this point will - // not be reached. - return default_order(left, right); - } - else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union) - && m_turn_points[right.turn_index].combination(operation_intersection, operation_union)) - { - return consider_iu_iu(left, right, "iu/iu"); - } - else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_blocked) - && m_turn_points[right.turn_index].combination(operation_intersection, operation_blocked)) - { - return consider_ix_ix(left, right, "ix/ix"); - } - else if (m_turn_points[left.turn_index].both(operation_intersection) - && m_turn_points[right.turn_index].both(operation_intersection)) - { - return consider_ii(left, right, "ii/ii"); - } - else if (m_turn_points[left.turn_index].combination(operation_union, operation_blocked) - && m_turn_points[right.turn_index].combination(operation_intersection, operation_union)) - { - return consider_iu_ux(left, right, -1, "ux/iu"); - } - else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union) - && m_turn_points[right.turn_index].combination(operation_union, operation_blocked)) - { - return consider_iu_ux(left, right, 1, "iu/ux"); - } - else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_blocked) - && m_turn_points[right.turn_index].combination(operation_intersection, operation_union)) - { - return consider_iu_ix(left, right, 1, "ix/iu"); - } - else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union) - && m_turn_points[right.turn_index].combination(operation_intersection, operation_blocked)) - { - return consider_iu_ix(left, right, -1, "iu/ix"); - } - else if (m_turn_points[left.turn_index].method != method_equal - && m_turn_points[right.turn_index].method == method_equal - ) - { - // If one of them was EQUAL or CONTINUES, it should always come first - return false; - } - else if (m_turn_points[left.turn_index].method == method_equal - && m_turn_points[right.turn_index].method != method_equal - ) - { - return true; - } - - // Now we have no clue how to sort. - -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - std::cout << " Consider: " << operation_char(m_turn_points[left.turn_index].operations[0].operation) - << operation_char(m_turn_points[left.turn_index].operations[1].operation) - << "/" << operation_char(m_turn_points[right.turn_index].operations[0].operation) - << operation_char(m_turn_points[right.turn_index].operations[1].operation) - << " " << " Take " << left.turn_index << " < " << right.turn_index - << std::endl; -#endif - - return default_order(left, right); - } -}; - - - -template -< - typename IndexType, - typename Iterator, - typename TurnPoints, - typename Geometry1, - typename Geometry2, - typename Strategy -> -inline void inspect_cluster(Iterator begin_cluster, Iterator end_cluster, - TurnPoints& turn_points, - operation_type , - Geometry1 const& , Geometry2 const& , - Strategy const& ) -{ - int count = 0; - - // Make an analysis about all occuring cases here. - std::map<std::pair<operation_type, operation_type>, int> inspection; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - operation_type first = turn_points[it->turn_index].operations[0].operation; - operation_type second = turn_points[it->turn_index].operations[1].operation; - if (first > second) - { - std::swap(first, second); - } - inspection[std::make_pair(first, second)]++; - count++; - } - - - bool keep_cc = false; - - // Decide about which is going to be discarded here. - if (inspection[std::make_pair(operation_union, operation_union)] == 1 - && inspection[std::make_pair(operation_continue, operation_continue)] == 1) - { - // In case of uu/cc, discard the uu, that indicates a tangency and - // inclusion would disturb the (e.g.) cc-cc-cc ordering - // NOTE: uu is now discarded anyhow. - keep_cc = true; - } - else if (count == 2 - && inspection[std::make_pair(operation_intersection, operation_intersection)] == 1 - && inspection[std::make_pair(operation_union, operation_intersection)] == 1) - { - // In case of ii/iu, discard the iu. The ii should always be visited, - // Because (in case of not discarding iu) correctly ordering of ii/iu appears impossible - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - if (turn_points[it->turn_index].combination(operation_intersection, operation_union)) - { - it->discarded = true; - } - } - } - - // Discard any continue turn, unless it is the only thing left - // (necessary to avoid cc-only rings, all being discarded - // e.g. traversal case #75) - int nd_count= 0, cc_count = 0; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - if (! it->discarded) - { - nd_count++; - if (turn_points[it->turn_index].both(operation_continue)) - { - cc_count++; - } - } - } - - if (nd_count == cc_count) - { - keep_cc = true; - } - - if (! keep_cc) - { - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - if (turn_points[it->turn_index].both(operation_continue)) - { - it->discarded = true; - } - } - } -} - - -template -< - typename IndexType, - bool Reverse1, bool Reverse2, - typename Iterator, - typename TurnPoints, - typename Geometry1, - typename Geometry2, - typename RobustPolicy, - typename Strategy -> -inline void handle_cluster(Iterator begin_cluster, Iterator end_cluster, - TurnPoints& turn_points, - operation_type for_operation, - Geometry1 const& geometry1, Geometry2 const& geometry2, - RobustPolicy& robust_policy, - Strategy const& strategy) -{ - // First inspect and (possibly) discard rows - inspect_cluster<IndexType>(begin_cluster, end_cluster, turn_points, - for_operation, geometry1, geometry2, strategy); - - - // Then sort this range (discarded rows will be ordered first and will be removed in enrich_assign) - std::sort(begin_cluster, end_cluster, - sort_in_cluster - < - TurnPoints, - IndexType, - Geometry1, Geometry2, - RobustPolicy, - Reverse1, Reverse2, - Strategy - >(turn_points, geometry1, geometry2, robust_policy, strategy)); - -#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES) - typedef typename IndexType::type operations_type; - operations_type const& op = turn_points[begin_cluster->turn_index].operations[begin_cluster->operation_index]; - std::cout << std::endl << "Clustered points on equal distance " << op.fraction << std::endl; - - std::cout << "->Indexes "; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - std::cout << " " << it->turn_index; - } - std::cout << std::endl << "->Methods: "; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - std::cout << " " << method_char(turn_points[it->turn_index].method); - } - std::cout << std::endl << "->Operations: "; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - std::cout << " " << operation_char(turn_points[it->turn_index].operations[0].operation) - << operation_char(turn_points[it->turn_index].operations[1].operation); - } - std::cout << std::endl << "->Discarded: "; - for (Iterator it = begin_cluster; it != end_cluster; ++it) - { - std::cout << " " << (it->discarded ? "true" : "false"); - } - std::cout << std::endl; - //<< "\tOn segments: " << prev_op.seg_id << " / " << prev_op.other_id - //<< " and " << op.seg_id << " / " << op.other_id - //<< geometry::distance(turn_points[prev->turn_index].point, turn_points[it->turn_index].point) -#endif - -} - - -}} // namespace detail::overlay -#endif //DOXYGEN_NO_DETAIL - - -}} // namespace boost::geometry - - -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TANGENCIES_HPP diff --git a/boost/geometry/algorithms/detail/overlay/handle_touch.hpp b/boost/geometry/algorithms/detail/overlay/handle_touch.hpp new file mode 100644 index 0000000000..1febefc83a --- /dev/null +++ b/boost/geometry/algorithms/detail/overlay/handle_touch.hpp @@ -0,0 +1,336 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// 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) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TOUCH_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TOUCH_HPP + +#include <cstddef> + +#include <boost/range.hpp> + +#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> +#include <boost/geometry/geometries/concepts/check.hpp> +#include <boost/geometry/algorithms/detail/ring_identifier.hpp> +#include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp> + + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace overlay +{ + + +template <typename Turns, typename Visitor> +class handle_touch_uu +{ +private : + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename boost::range_iterator<Turns>::type turn_iterator; + typedef typename boost::range_iterator<Turns const>::type turn_const_iterator; + + typedef typename boost::range_iterator + < + typename turn_type::container_type const + >::type operation_const_iterator; + +public : + + handle_touch_uu(Visitor& visitor) + : m_visitor(visitor) + {} + + inline void apply(detail::overlay::operation_type operation, Turns& turns) + { + if (! has_uu(turns)) + { + // Performance - if there is no u/u at all, nothing to be done + return; + } + + // Iterate through all u/u points + int turn_index = 0; + for (turn_iterator it = boost::begin(turns); + it != boost::end(turns); + ++it, ++turn_index) + { + turn_type& turn = *it; + if (! turn.both(operation_union)) + { + continue; + } + + m_visitor.print("handle_touch uu:", turns, turn_index); + + bool const traverse = turn_should_be_traversed(turns, turn, turn_index); + bool const start = traverse + && turn_should_be_startable(turns, turn, turn_index); + m_visitor.print("handle_touch, ready ", turns, turn_index); +// << std::boolalpha +// << traverse << " " << start + + if (traverse) + { + // Indicate the sources should switch here to create + // separate rings (outer ring / inner ring) + turn.switch_source = true; + } + // TODO: this is often not correct, fix this + turn.operations[0].enriched.startable = start; + turn.operations[1].enriched.startable = start; + } + } + +private : + + // Generic utility to be moved somewhere else + static inline + ring_identifier ring_id_from_seg_id(const segment_identifier& seg_id) + { + return ring_identifier(seg_id.source_index, + seg_id.multi_index, + seg_id.ring_index); + } + + static inline + ring_identifier ring_id_from_op(const turn_type& turn, + int operation_index) + { + return ring_id_from_seg_id(turn.operations[operation_index].seg_id); + } + + static inline bool in_range(const Turns& turns, signed_size_type index) + { + signed_size_type const turns_size = + static_cast<signed_size_type>(boost::size(turns)); + return index >= 0 && index < turns_size; + } + + static inline bool has_uu(const Turns& turns) + { + for (turn_const_iterator it = boost::begin(turns); + it != boost::end(turns); + ++it) + { + const turn_type& turn = *it; + if (turn.both(operation_union)) + { + return true; + } + } + return false; + } + + static inline + bool turn_should_be_startable(const Turns& turns, + const turn_type& uu_turn, + signed_size_type uu_turn_index) + { + return turn_startable(turns, uu_turn, 0, uu_turn_index) + || turn_startable(turns, uu_turn, 1, uu_turn_index); + + } + + static inline + bool turn_startable(const Turns& turns, + const turn_type& uu_turn, + std::size_t op_index, + signed_size_type original_turn_index, + std::size_t iteration = 0) + { + if (iteration >= boost::size(turns)) + { + // Defensive check to avoid infinite recursion + return false; + } + + signed_size_type const index + = uu_turn.operations[op_index].enriched.travels_to_ip_index; + if (index == original_turn_index) + { + // Completely traveled, having u/u only, via this op_index + return true; + } + + if (! in_range(turns, index)) + { + return false; + } + + const turn_type& new_turn = turns[index]; + + if (new_turn.operations[0].enriched.startable) + { + // Already selectable - no need to select u/u turn too + return false; + } + + // If this u/u turn is traversed normally (without skipping), sources are switched + return turn_startable(turns, new_turn, 1 - op_index, + original_turn_index, iteration + 1); + } + + inline bool turn_should_be_traversed(const Turns& turns, + const turn_type& uu_turn, + signed_size_type uu_turn_index) + { + return turn_should_be_traversed(turns, uu_turn, uu_turn_index, 0) + || turn_should_be_traversed(turns, uu_turn, uu_turn_index, 1); + } + + inline bool turn_should_be_traversed(const Turns& turns, + const turn_type& uu_turn, + signed_size_type uu_turn_index, + int uu_operation_index) + { + // Suppose this is a u/u turn between P and Q + // Examine all other turns on P and check if Q can be reached + // Use one of the operations and check if you can reach the other + signed_size_type const to_turn_index + = uu_turn.operations[uu_operation_index].enriched.travels_to_ip_index; + if (! in_range(turns, to_turn_index)) + { + return false; + } + + m_visitor.print("Examine:", turns, to_turn_index); + ring_identifier const other_ring_id + = ring_id_from_op(uu_turn, 1 - uu_operation_index); + + bool complete = false; + return can_reach(complete, turns, turns[to_turn_index], uu_operation_index, + other_ring_id, uu_turn_index, to_turn_index); + } + + inline bool can_reach(bool& complete, const Turns& turns, + const turn_type& turn, + signed_size_type uu_operation_index, + const ring_identifier& target_ring_id, + signed_size_type uu_turn_index, + signed_size_type to_turn_index, + std::size_t iteration = 0) + { + if (complete) + { + return false; + } + + if (turn.cluster_id >= 0) + { + // Clustered turns are yet not supported + return false; + } + + if (iteration >= boost::size(turns)) + { + m_visitor.print("Too much iterations"); + // Defensive check to avoid infinite recursion + return false; + } + + if (uu_operation_index != -1 && turn.both(operation_union)) + { + // If we end up in a u/u turn, check the way how, for this operation + m_visitor.print("Via u/u"); + return can_reach_via(complete, turns, uu_operation_index, + turn.operations[uu_operation_index], + target_ring_id, + uu_turn_index, to_turn_index, iteration); + } + else + { + // Check if specified ring can be reached via one of both operations + return can_reach_via(complete, turns, 0, turn.operations[0], target_ring_id, + uu_turn_index, to_turn_index, iteration) + || can_reach_via(complete, turns, 1, turn.operations[1], target_ring_id, + uu_turn_index, to_turn_index, iteration); + } + } + + template <typename Operation> + inline bool can_reach_via(bool& complete, const Turns& turns, + signed_size_type operation_index, + const Operation& operation, + const ring_identifier& target_ring_id, + signed_size_type uu_turn_index, + signed_size_type to_turn_index, + std::size_t iteration = 0) + { + if (operation.operation != operation_union + && operation.operation != operation_continue) + { + return false; + } + + signed_size_type const index = operation.enriched.travels_to_ip_index; + if (index == to_turn_index) + { + m_visitor.print("Dead end at", turns, index); + // Completely traveled, the target is not found + return false; + } + if (index == uu_turn_index) + { + // End up where trial was started + m_visitor.print("Travel complete at", turns, index); + complete = true; + return false; + } + + if (! in_range(turns, index)) + { + return false; + } + + m_visitor.print("Now to", turns, index, operation_index); + const turn_type& new_turn = turns[index]; + + if (new_turn.both(operation_union)) + { + ring_identifier const ring_id = ring_id_from_op(new_turn, operation_index); + if (ring_id == target_ring_id) + { + m_visitor.print("Found (at u/u)!"); + return true; + } + } + else + { + ring_identifier const ring_id1 = ring_id_from_op(new_turn, 0); + ring_identifier const ring_id2 = ring_id_from_op(new_turn, 1); + if (ring_id1 == target_ring_id || ring_id2 == target_ring_id) + { + m_visitor.print("Found!"); + return true; + } + } + + // Recursively check this turn + return can_reach(complete, turns, new_turn, operation_index, target_ring_id, + uu_turn_index, to_turn_index, iteration + 1); + } + +private : + Visitor m_visitor; +}; + +template <typename Turns, typename Visitor> +inline void handle_touch(detail::overlay::operation_type operation, + Turns& turns, Visitor& visitor) +{ + handle_touch_uu<Turns, Visitor> handler(visitor); + handler.apply(operation, turns); +} + +}} // namespace detail::overlay +#endif // DOXYGEN_NO_DETAIL + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TOUCH_HPP diff --git a/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp new file mode 100644 index 0000000000..21868a2939 --- /dev/null +++ b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp @@ -0,0 +1,203 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. + +// 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) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_ON_SEGMENT_RATIO_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_ON_SEGMENT_RATIO_HPP + +#include <cstddef> +#include <algorithm> +#include <map> +#include <set> +#include <vector> + +#include <boost/range.hpp> + +#include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp> +#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> +#include <boost/geometry/strategies/side.hpp> + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace overlay +{ + +// Wraps "turn_operation" from turn_info.hpp, +// giving it extra information, necessary for sorting +template <typename TurnOperation> +struct indexed_turn_operation +{ + typedef TurnOperation type; + + std::size_t turn_index; + std::size_t operation_index; + bool skip; + // use pointers to avoid copies, const& is not possible because of usage in vector + segment_identifier const* other_seg_id; // segment id of other segment of intersection of two segments + TurnOperation const* subject; + + inline indexed_turn_operation(std::size_t ti, std::size_t oi, + TurnOperation const& sub, + segment_identifier const& oid) + : turn_index(ti) + , operation_index(oi) + , skip(false) + , other_seg_id(&oid) + , subject(boost::addressof(sub)) + {} + +}; + +template +< + typename Turns, + typename Indexed, + typename Geometry1, typename Geometry2, + typename RobustPolicy, + bool Reverse1, bool Reverse2 +> +struct less_by_segment_ratio +{ + inline less_by_segment_ratio(Turns const& turns + , operation_type for_operation + , Geometry1 const& geometry1 + , Geometry2 const& geometry2 + , RobustPolicy const& robust_policy) + : m_turns(turns) + , m_for_operation(for_operation) + , m_geometry1(geometry1) + , m_geometry2(geometry2) + , m_robust_policy(robust_policy) + { + } + +private : + + Turns const& m_turns; + operation_type m_for_operation; + Geometry1 const& m_geometry1; + Geometry2 const& m_geometry2; + RobustPolicy const& m_robust_policy; + + typedef typename geometry::point_type<Geometry1>::type point_type; + + inline bool default_order(Indexed const& left, Indexed const& right) const + { + // We've nothing to sort on. Take the indexes + return left.turn_index < right.turn_index; + } + + inline bool consider_relative_order(Indexed const& left, + Indexed const& right) const + { + point_type pi, pj, ri, rj, si, sj; + + geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, + left.subject->seg_id, + pi, pj); + geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, + *left.other_seg_id, + ri, rj); + geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2, + *right.other_seg_id, + si, sj); + + typedef typename strategy::side::services::default_strategy + < + typename cs_tag<point_type>::type + >::type strategy; + + int const side_rj_p = strategy::apply(pi, pj, rj); + int const side_sj_p = strategy::apply(pi, pj, sj); + + // Put the one turning left (1; right == -1) as last + if (side_rj_p != side_sj_p) + { + return side_rj_p < side_sj_p; + } + + int const side_sj_r = strategy::apply(ri, rj, sj); + int const side_rj_s = strategy::apply(si, sj, rj); + + // If they both turn left: the most left as last + // If they both turn right: this is not relevant, but take also here most left + if (side_rj_s != side_sj_r) + { + return side_rj_s < side_sj_r; + } + + return default_order(left, right); + } + + +public : + + // Note that left/right do NOT correspond to m_geometry1/m_geometry2 + // but to the "indexed_turn_operation" + inline bool operator()(Indexed const& left, Indexed const& right) const + { + if (! (left.subject->seg_id == right.subject->seg_id)) + { + return left.subject->seg_id < right.subject->seg_id; + } + + // Both left and right are located on the SAME segment. + + if (! (left.subject->fraction == right.subject->fraction)) + { + return left.subject->fraction < right.subject->fraction; + } + + + typedef typename boost::range_value<Turns>::type turn_type; + turn_type const& left_turn = m_turns[left.turn_index]; + turn_type const& right_turn = m_turns[right.turn_index]; + + // First check "real" intersection (crosses) + // -> distance zero due to precision, solve it by sorting + if (left_turn.method == method_crosses + && right_turn.method == method_crosses) + { + return consider_relative_order(left, right); + } + + bool const left_both_xx = left_turn.both(operation_blocked); + bool const right_both_xx = right_turn.both(operation_blocked); + if (left_both_xx && ! right_both_xx) + { + return true; + } + if (! left_both_xx && right_both_xx) + { + return false; + } + + bool const left_both_uu = left_turn.both(operation_union); + bool const right_both_uu = right_turn.both(operation_union); + if (left_both_uu && ! right_both_uu) + { + return true; + } + if (! left_both_uu && right_both_uu) + { + return false; + } + + return default_order(left, right); + } +}; + + +}} // namespace detail::overlay +#endif //DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_ON_SEGMENT_RATIO_HPP diff --git a/boost/geometry/algorithms/detail/overlay/overlay.hpp b/boost/geometry/algorithms/detail/overlay/overlay.hpp index 6eb0b8864c..c3ecaa0b01 100644 --- a/boost/geometry/algorithms/detail/overlay/overlay.hpp +++ b/boost/geometry/algorithms/detail/overlay/overlay.hpp @@ -26,6 +26,7 @@ #include <boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp> #include <boost/geometry/algorithms/detail/overlay/enrichment_info.hpp> #include <boost/geometry/algorithms/detail/overlay/get_turns.hpp> +#include <boost/geometry/algorithms/detail/overlay/handle_touch.hpp> #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp> #include <boost/geometry/algorithms/detail/overlay/traverse.hpp> #include <boost/geometry/algorithms/detail/overlay/traversal_info.hpp> @@ -59,28 +60,50 @@ namespace detail { namespace overlay { -template <typename TurnPoints, typename TurnInfoMap> -inline void get_ring_turn_info(TurnInfoMap& turn_info_map, - TurnPoints const& turn_points) +//! Default visitor for overlay, doing nothing +struct overlay_null_visitor { - typedef typename boost::range_value<TurnPoints>::type turn_point_type; - typedef typename turn_point_type::container_type container_type; + void print(char const* ) {} - for (typename boost::range_iterator<TurnPoints const>::type - it = boost::begin(turn_points); - it != boost::end(turn_points); + template <typename Turns> + void print(char const* , Turns const& , int) {} + + template <typename Turns> + void print(char const* , Turns const& , int , int ) {} + + template <typename Turns> + void visit_turns(int , Turns const& ) {} + + template <typename Clusters, typename Turns> + void visit_clusters(Clusters const& , Turns const& ) {} + + template <typename Turns, typename Turn, typename Operation> + void visit_traverse(Turns const& , Turn const& , Operation const& , char const*) + {} + + template <typename Turns, typename Turn, typename Operation> + void visit_traverse_reject(Turns const& , Turn const& , Operation const& , traverse_error_type ) + {} +}; + +template <typename Turns, typename TurnInfoMap> +inline void get_ring_turn_info(TurnInfoMap& turn_info_map, Turns const& turns) +{ + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::container_type container_type; + + for (typename boost::range_iterator<Turns const>::type + it = boost::begin(turns); + it != boost::end(turns); ++it) { - typename boost::range_value<TurnPoints>::type const& turn_info = *it; - bool both_uu = turn_info.both(operation_union); - bool skip = (turn_info.discarded || both_uu) - && ! turn_info.any_blocked() - && ! turn_info.both(operation_intersection) - ; + typename boost::range_value<Turns>::type const& turn_info = *it; - if (! both_uu && turn_info.colocated) + if (turn_info.discarded + && ! turn_info.any_blocked() + && ! turn_info.colocated) { - skip = true; + continue; } for (typename boost::range_iterator<container_type const>::type @@ -88,21 +111,13 @@ inline void get_ring_turn_info(TurnInfoMap& turn_info_map, op_it != boost::end(turn_info.operations); ++op_it) { - ring_identifier ring_id + ring_identifier const ring_id ( op_it->seg_id.source_index, op_it->seg_id.multi_index, op_it->seg_id.ring_index ); - - if (! skip) - { - turn_info_map[ring_id].has_normal_turn = true; - } - else if (both_uu) - { - turn_info_map[ring_id].has_uu_turn = true; - } + turn_info_map[ring_id].has_normal_turn = true; } } } @@ -164,12 +179,13 @@ template > struct overlay { - template <typename RobustPolicy, typename OutputIterator, typename Strategy> + template <typename RobustPolicy, typename OutputIterator, typename Strategy, typename Visitor> static inline OutputIterator apply( Geometry1 const& geometry1, Geometry2 const& geometry2, RobustPolicy const& robust_policy, OutputIterator out, - Strategy const& ) + Strategy const& , + Visitor& visitor) { bool const is_empty1 = geometry::is_empty(geometry1); bool const is_empty2 = geometry::is_empty(geometry2); @@ -193,14 +209,23 @@ struct overlay point_type, typename geometry::segment_ratio_type<point_type, RobustPolicy>::type > turn_info; - typedef std::deque<turn_info> container_type; + typedef std::deque<turn_info> turn_container_type; typedef std::deque < typename geometry::ring_type<GeometryOut>::type > ring_container_type; - container_type turn_points; + // Define the clusters, mapping cluster_id -> turns + typedef std::map + < + signed_size_type, + std::set<signed_size_type> + > cluster_type; + + cluster_type clusters; + + turn_container_type turns; #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE std::cout << "get turns" << std::endl; @@ -210,20 +235,42 @@ std::cout << "get turns" << std::endl; < Reverse1, Reverse2, detail::overlay::assign_null_policy - >(geometry1, geometry2, robust_policy, turn_points, policy); + >(geometry1, geometry2, robust_policy, turns, policy); + + visitor.visit_turns(1, turns); + + static const operation_type op_type + = OverlayType == overlay_union + ? geometry::detail::overlay::operation_union + : geometry::detail::overlay::operation_intersection; #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE std::cout << "enrich" << std::endl; #endif typename Strategy::side_strategy_type side_strategy; - geometry::enrich_intersection_points<Reverse1, Reverse2, OverlayType>(turn_points, - OverlayType == overlay_union - ? geometry::detail::overlay::operation_union - : geometry::detail::overlay::operation_intersection, + geometry::enrich_intersection_points<Reverse1, Reverse2, OverlayType>(turns, + clusters, op_type, geometry1, geometry2, robust_policy, side_strategy); + visitor.visit_turns(2, turns); + + visitor.visit_clusters(clusters, turns); + + +#if 0 + // TODO: does not work always correctly, move to traverse and fix + if (op_type == geometry::detail::overlay::operation_union) + { + #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE + std::cout << "handle_touch" << std::endl; + #endif + + handle_touch(op_type, turns, visitor); + } +#endif + #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE std::cout << "traverse" << std::endl; #endif @@ -231,18 +278,17 @@ std::cout << "traverse" << std::endl; // Note that these rings are always in clockwise order, even in CCW polygons, // and are marked as "to be reversed" below ring_container_type rings; - traverse<Reverse1, Reverse2, Geometry1, Geometry2>::apply + traverse<Reverse1, Reverse2, Geometry1, Geometry2, op_type>::apply ( geometry1, geometry2, - OverlayType == overlay_union - ? geometry::detail::overlay::operation_union - : geometry::detail::overlay::operation_intersection, robust_policy, - turn_points, rings + turns, rings, + clusters, + visitor ); std::map<ring_identifier, ring_turn_info> turn_info_per_ring; - get_ring_turn_info(turn_info_per_ring, turn_points); + get_ring_turn_info(turn_info_per_ring, turns); typedef ring_properties < @@ -272,6 +318,17 @@ std::cout << "traverse" << std::endl; return add_rings<GeometryOut>(selected_ring_properties, geometry1, geometry2, rings, out); } + + template <typename RobustPolicy, typename OutputIterator, typename Strategy> + static inline OutputIterator apply( + Geometry1 const& geometry1, Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + OutputIterator out, + Strategy const& strategy) + { + overlay_null_visitor visitor; + return apply(geometry1, geometry2, robust_policy, out, strategy, visitor); + } }; diff --git a/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp b/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp index e77a163dd5..5447c8813e 100644 --- a/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp +++ b/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp @@ -38,6 +38,7 @@ struct segment_identifier , multi_index(-1) , ring_index(-1) , segment_index(-1) + , piece_index(-1) {} inline segment_identifier(signed_size_type src, @@ -48,6 +49,7 @@ struct segment_identifier , multi_index(mul) , ring_index(rin) , segment_index(seg) + , piece_index(-1) {} inline bool operator<(segment_identifier const& other) const @@ -85,6 +87,9 @@ struct segment_identifier signed_size_type multi_index; signed_size_type ring_index; signed_size_type segment_index; + + // For buffer - todo: move this to buffer-only + signed_size_type piece_index; }; diff --git a/boost/geometry/algorithms/detail/overlay/select_rings.hpp b/boost/geometry/algorithms/detail/overlay/select_rings.hpp index d18e012b2d..1b3cd866d7 100644 --- a/boost/geometry/algorithms/detail/overlay/select_rings.hpp +++ b/boost/geometry/algorithms/detail/overlay/select_rings.hpp @@ -35,13 +35,11 @@ namespace detail { namespace overlay struct ring_turn_info { - bool has_uu_turn; bool has_normal_turn; bool within_other; ring_turn_info() - : has_uu_turn(false) - , has_normal_turn(false) + : has_normal_turn(false) , within_other(false) {} }; @@ -174,7 +172,7 @@ struct decide<overlay_union> { static bool include(ring_identifier const& , ring_turn_info const& info) { - return info.has_uu_turn || ! info.within_other; + return ! info.within_other; } static bool reversed(ring_identifier const& , ring_turn_info const& ) @@ -190,19 +188,16 @@ struct decide<overlay_difference> { // Difference: A - B - // If this is A (source_index=0) and there is only a u/u turn, - // then the ring is inside B - // If this is B (source_index=1) and there is only a u/u turn, - // then the ring is NOT inside A + // If this is A (source_index=0), then the ring is inside B + // If this is B (source_index=1), then the ring is NOT inside A // If this is A and the ring is within the other geometry, // then we should NOT include it. // If this is B then we SHOULD include it. - bool const is_first = id.source_index == 0; - bool const within_other = info.within_other - || (is_first && info.has_uu_turn); - return is_first ? ! within_other : within_other; + return id.source_index == 0 + ? ! info.within_other + : info.within_other; } static bool reversed(ring_identifier const& id, ring_turn_info const& info) @@ -219,7 +214,7 @@ struct decide<overlay_intersection> { static bool include(ring_identifier const& , ring_turn_info const& info) { - return ! info.has_uu_turn && info.within_other; + return info.within_other; } static bool reversed(ring_identifier const& , ring_turn_info const& ) @@ -266,19 +261,16 @@ inline void update_ring_selection(Geometry1 const& geometry1, continue; } - if (! info.has_uu_turn) + // Check if the ring is within the other geometry, by taking + // a point lying on the ring + switch(id.source_index) { - // Check if the ring is within the other geometry, by taking - // a point lying on the ring - switch(id.source_index) - { - case 0 : - info.within_other = geometry::within(it->second.point, geometry2); - break; - case 1 : - info.within_other = geometry::within(it->second.point, geometry1); - break; - } + case 0 : + info.within_other = geometry::within(it->second.point, geometry2); + break; + case 1 : + info.within_other = geometry::within(it->second.point, geometry1); + break; } if (decide<OverlayType>::include(id, info)) diff --git a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp new file mode 100644 index 0000000000..894cddab87 --- /dev/null +++ b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp @@ -0,0 +1,434 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// 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) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP + +#include <algorithm> +#include <vector> + +#include <boost/geometry/algorithms/detail/direction_code.hpp> +#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> +#include <boost/geometry/strategies/side.hpp> + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace overlay { namespace sort_by_side +{ + +enum index_type { index_unknown = -1, index_from = 0, index_to = 1 }; + +// Point-wrapper, adding some properties +template <typename Point> +struct ranked_point +{ + ranked_point() + : main_rank(0) + , turn_index(-1) + , op_index(-1) + , index(index_unknown) + , left_count(0) + , right_count(0) + , operation(operation_none) + {} + + ranked_point(const Point& p, signed_size_type ti, signed_size_type oi, + index_type i, operation_type op, segment_identifier sid) + : point(p) + , main_rank(0) + , turn_index(ti) + , op_index(oi) + , index(i) + , left_count(0) + , right_count(0) + , operation(op) + , seg_id(sid) + {} + + Point point; + std::size_t main_rank; + signed_size_type turn_index; + signed_size_type op_index; + index_type index; + std::size_t left_count; + std::size_t right_count; + operation_type operation; + segment_identifier seg_id; +}; + +struct less_by_turn_index +{ + template <typename T> + inline bool operator()(const T& first, const T& second) const + { + return first.turn_index == second.turn_index + ? first.index < second.index + : first.turn_index < second.turn_index + ; + } +}; + +struct less_by_index +{ + template <typename T> + inline bool operator()(const T& first, const T& second) const + { + // First order by from/to + if (first.index != second.index) + { + return first.index < second.index; + } + // All the same, order by turn index (we might consider length too) + return first.turn_index < second.turn_index; + } +}; + +struct less_false +{ + template <typename T> + inline bool operator()(const T&, const T& ) const + { + return false; + } +}; + +template <typename Point, typename LessOnSame, typename Compare> +struct less_by_side +{ + typedef typename strategy::side::services::default_strategy + < + typename cs_tag<Point>::type + >::type side; + + less_by_side(const Point& p1, const Point& p2) + : m_p1(p1) + , m_p2(p2) + {} + + template <typename T> + inline bool operator()(const T& first, const T& second) const + { + LessOnSame on_same; + Compare compare; + + int const side_first = side::apply(m_p1, m_p2, first.point); + int const side_second = side::apply(m_p1, m_p2, second.point); + + if (side_first == 0 && side_second == 0) + { + // Both collinear. They might point into different directions: <------*------> + // If so, order the one going backwards as the very first. + + int const first_code = direction_code(m_p1, m_p2, first.point); + int const second_code = direction_code(m_p1, m_p2, second.point); + + // Order by code, backwards first, then forward. + return first_code != second_code + ? first_code < second_code + : on_same(first, second) + ; + } + else if (side_first == 0 + && direction_code(m_p1, m_p2, first.point) == -1) + { + // First collinear and going backwards. + // Order as the very first, so return always true + return true; + } + else if (side_second == 0 + && direction_code(m_p1, m_p2, second.point) == -1) + { + // Second is collinear and going backwards + // Order as very last, so return always false + return false; + } + + // They are not both collinear + + if (side_first != side_second) + { + return compare(side_first, side_second); + } + + // They are both left, both right, and/or both collinear (with each other and/or with p1,p2) + // Check mutual side + int const side_second_wrt_first = side::apply(m_p2, first.point, second.point); + + if (side_second_wrt_first == 0) + { + return on_same(first, second); + } + + int const side_first_wrt_second = -side_second_wrt_first; + + // Both are on same side, and not collinear + // Union: return true if second is right w.r.t. first, so -1, + // so other is 1. union has greater as compare functor + // Intersection: v.v. + return compare(side_first_wrt_second, side_second_wrt_first); + } + +private : + Point m_p1, m_p2; +}; + +template <bool Reverse1, bool Reverse2, typename Point, typename Compare> +struct side_sorter +{ + typedef ranked_point<Point> rp; + + inline void set_origin(Point const& origin) + { + m_origin = origin; + } + + template <typename Operation, typename Geometry1, typename Geometry2> + void add(Operation const& op, signed_size_type turn_index, signed_size_type op_index, + Geometry1 const& geometry1, + Geometry2 const& geometry2, + bool is_origin) + { + Point point1, point2, point3; + geometry::copy_segment_points<Reverse1, Reverse2>(geometry1, geometry2, + op.seg_id, point1, point2, point3); + Point const& point_to = op.fraction.is_one() ? point3 : point2; + + m_ranked_points.push_back(rp(point1, turn_index, op_index, index_from, op.operation, op.seg_id)); + m_ranked_points.push_back(rp(point_to, turn_index, op_index, index_to, op.operation, op.seg_id)); + + if (is_origin) + { + m_origin = point1; + } + } + + void apply(Point const& turn_point) + { + // We need three compare functors: + // 1) to order clockwise (union) or counter clockwise (intersection) + // 2) to order by side, resulting in unique ranks for all points + // 3) to order by side, resulting in non-unique ranks + // to give colinear points + + // Sort by side and assign rank + less_by_side<Point, less_by_index, Compare> less_unique(m_origin, turn_point); + less_by_side<Point, less_false, Compare> less_non_unique(m_origin, turn_point); + + std::sort(m_ranked_points.begin(), m_ranked_points.end(), less_unique); + + std::size_t colinear_rank = 0; + for (std::size_t i = 0; i < m_ranked_points.size(); i++) + { + if (i > 0 + && less_non_unique(m_ranked_points[i - 1], m_ranked_points[i])) + { + // It is not collinear + colinear_rank++; + } + + m_ranked_points[i].main_rank = colinear_rank; + } + } + + template <signed_size_type segment_identifier::*Member, typename Map> + void find_open_generic(Map& handled) + { + for (std::size_t i = 0; i < m_ranked_points.size(); i++) + { + const rp& ranked = m_ranked_points[i]; + if (ranked.index != index_from) + { + continue; + } + + signed_size_type const& index = ranked.seg_id.*Member; + if (! handled[index]) + { + find_polygons_for_source<Member>(index, i); + handled[index] = true; + } + } + } + + void find_open() + { + // TODO: we might pass Buffer as overlay_type, instead on the fly below + bool one_source = true; + for (std::size_t i = 0; i < m_ranked_points.size(); i++) + { + const rp& ranked = m_ranked_points[i]; + signed_size_type const& src = ranked.seg_id.source_index; + if (src != 0) + { + one_source = false; + break; + } + } + + if (one_source) + { + // by multi index + std::map<signed_size_type, bool> handled; + find_open_generic + < + &segment_identifier::piece_index + >(handled); + } + else + { + // by source (there should only source 0,1) TODO assert this + bool handled[2] = {false, false}; + find_open_generic + < + &segment_identifier::source_index + >(handled); + } + } + + void reverse() + { + if (m_ranked_points.empty()) + { + return; + } + + int const last = 1 + m_ranked_points.back().main_rank; + + // Move iterator after main_rank==0 + bool has_first = false; + typename container_type::iterator it = m_ranked_points.begin() + 1; + for (; it != m_ranked_points.end() && it->main_rank == 0; ++it) + { + has_first = true; + } + + if (has_first) + { + // Reverse first part (having main_rank == 0), if any, + // but skip the very first row + std::reverse(m_ranked_points.begin() + 1, it); + for (typename container_type::iterator fit = m_ranked_points.begin(); + fit != it; ++fit) + { + BOOST_ASSERT(fit->main_rank == 0); + } + } + + // Reverse the rest (main rank > 0) + std::reverse(it, m_ranked_points.end()); + for (; it != m_ranked_points.end(); ++it) + { + BOOST_ASSERT(it->main_rank > 0); + it->main_rank = last - it->main_rank; + } + } + +//protected : + + typedef std::vector<rp> container_type; + container_type m_ranked_points; + Point m_origin; + +private : + + + std::size_t move(std::size_t index) const + { + std::size_t const result = index + 1; + return result >= m_ranked_points.size() ? 0 : result; + } + + //! member is pointer to member (source_index or multi_index) + template <signed_size_type segment_identifier::*Member> + std::size_t move(signed_size_type member_index, std::size_t index) const + { + std::size_t result = move(index); + while (m_ranked_points[result].seg_id.*Member != member_index) + { + result = move(result); + } + return result; + } + + void assign_ranks(std::size_t min_rank, std::size_t max_rank, int side_index) + { + for (std::size_t i = 0; i < m_ranked_points.size(); i++) + { + rp& ranked = m_ranked_points[i]; + // Suppose there are 8 ranks, if min=4,max=6: assign 4,5,6 + // if min=5,max=2: assign from 5,6,7,1,2 + bool const in_range + = max_rank >= min_rank + ? ranked.main_rank >= min_rank && ranked.main_rank <= max_rank + : ranked.main_rank >= min_rank || ranked.main_rank <= max_rank + ; + + if (in_range) + { + if (side_index == 1) + { + ranked.left_count++; + } + else if (side_index == 2) + { + ranked.right_count++; + } + } + } + } + + template <signed_size_type segment_identifier::*Member> + void find_polygons_for_source(signed_size_type the_index, + std::size_t start_index) + { + int state = 1; // 'closed', because start_index is "from", arrives at the turn + std::size_t last_from_rank = m_ranked_points[start_index].main_rank; + std::size_t previous_rank = m_ranked_points[start_index].main_rank; + + for (std::size_t index = move<Member>(the_index, start_index); + ; + index = move<Member>(the_index, index)) + { + rp& ranked = m_ranked_points[index]; + + if (ranked.main_rank != previous_rank && state == 0) + { + assign_ranks(last_from_rank, previous_rank - 1, 1); + assign_ranks(last_from_rank + 1, previous_rank, 2); + } + + if (index == start_index) + { + return; + } + + if (ranked.index == index_from) + { + last_from_rank = ranked.main_rank; + state++; + } + else if (ranked.index == index_to) + { + state--; + } + + previous_rank = ranked.main_rank; + } + } +}; + + +}}} // namespace detail::overlay::sort_by_side +#endif //DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP diff --git a/boost/geometry/algorithms/detail/overlay/traverse.hpp b/boost/geometry/algorithms/detail/overlay/traverse.hpp index 45e15d13d0..a8f4232124 100644 --- a/boost/geometry/algorithms/detail/overlay/traverse.hpp +++ b/boost/geometry/algorithms/detail/overlay/traverse.hpp @@ -15,6 +15,7 @@ #include <boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp> #include <boost/geometry/algorithms/detail/overlay/copy_segments.hpp> +#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> #include <boost/geometry/algorithms/num_points.hpp> #include <boost/geometry/core/access.hpp> @@ -65,186 +66,740 @@ inline void debug_traverse(Turn const& , Operation, const char*) #endif -template <typename Info, typename Turn> -inline void set_visited_for_continue(Info& info, Turn const& turn) +//! Metafunction to define side_order (clockwise, ccw) by operation_type +template <operation_type OpType> +struct side_compare {}; + +template <> +struct side_compare<operation_union> { - // On "continue", set "visited" for ALL directions - if (turn.operation == detail::overlay::operation_continue) + typedef std::greater<int> type; +}; + +template <> +struct side_compare<operation_intersection> +{ + typedef std::less<int> type; +}; + + +template +< + bool Reverse1, + bool Reverse2, + operation_type OperationType, + typename Geometry1, + typename Geometry2, + typename Turns, + typename Clusters, + typename RobustPolicy, + typename Visitor, + typename Backtrack +> +struct traversal +{ + typedef typename side_compare<OperationType>::type side_compare_type; + typedef typename boost::range_value<Turns>::type turn_type; + typedef typename turn_type::turn_operation_type turn_operation_type; + + typedef typename geometry::point_type<Geometry1>::type point_type; + typedef sort_by_side::side_sorter + < + Reverse1, Reverse2, + point_type, side_compare_type + > sbs_type; + + inline traversal(Geometry1 const& geometry1, Geometry2 const& geometry2, + Turns& turns, Clusters const& clusters, + RobustPolicy const& robust_policy, Visitor& visitor) + : m_geometry1(geometry1) + , m_geometry2(geometry2) + , m_turns(turns) + , m_clusters(clusters) + , m_robust_policy(robust_policy) + , m_visitor(visitor) + , m_has_uu(false) + , m_has_only_uu(true) + , m_switch_at_uu(true) + {} + + + inline bool select_source(signed_size_type turn_index, + segment_identifier const& seg_id1, + segment_identifier const& seg_id2) { - for (typename boost::range_iterator - < - typename Info::container_type - >::type it = boost::begin(info.operations); - it != boost::end(info.operations); - ++it) + if (OperationType == operation_intersection) { - if (it->visited.none()) + // For intersections always switch sources + return seg_id1.source_index != seg_id2.source_index; + } + else if (OperationType == operation_union) + { + // For uu, only switch sources if indicated + turn_type const& turn = m_turns[turn_index]; + + // TODO: pass this information + bool const is_buffer + = turn.operations[0].seg_id.source_index + == turn.operations[1].seg_id.source_index; + + if (is_buffer) { - it->visited.set_visited(); + // Buffer does not use source_index (always 0) + return turn.switch_source + ? seg_id1.multi_index != seg_id2.multi_index + : seg_id1.multi_index == seg_id2.multi_index; } + + // Temporarily use m_switch_at_uu, which does not solve all cases, + // but the majority of the more simple cases, making the interior + // rings valid + return m_switch_at_uu // turn.switch_source + ? seg_id1.source_index != seg_id2.source_index + : seg_id1.source_index == seg_id2.source_index; } + return false; } -} + inline + signed_size_type get_next_turn_index(turn_operation_type const& op) const + { + return op.enriched.next_ip_index == -1 + ? op.enriched.travels_to_ip_index + : op.enriched.next_ip_index; + } -template -< - bool Reverse1, bool Reverse2, - typename GeometryOut, - typename G1, - typename G2, - typename Turns, - typename IntersectionInfo, - typename RobustPolicy -> -inline bool assign_next_ip(G1 const& g1, G2 const& g2, - Turns& turns, - typename boost::range_iterator<Turns>::type& ip, - GeometryOut& current_output, - IntersectionInfo& info, - segment_identifier& seg_id, - RobustPolicy const& robust_policy) -{ - info.visited.set_visited(); - set_visited_for_continue(*ip, info); + inline bool traverse_possible(signed_size_type turn_index) const + { + if (turn_index == -1) + { + return false; + } + + turn_type const& turn = m_turns[turn_index]; + + // It is not a dead end if there is an operation to continue, or of + // there is a cluster (assuming for now we can get out of the cluster) + return turn.cluster_id >= 0 + || turn.has(OperationType) + || turn.has(operation_continue); + } - // If there is no next IP on this segment - if (info.enriched.next_ip_index < 0) + inline bool select_operation(turn_type& turn, + signed_size_type start_turn_index, + segment_identifier const& seg_id, + int& selected_op_index) { - if (info.enriched.travels_to_vertex_index < 0 - || info.enriched.travels_to_ip_index < 0) + if (turn.discarded) { return false; } - BOOST_GEOMETRY_ASSERT(info.enriched.travels_to_vertex_index >= 0); - BOOST_GEOMETRY_ASSERT(info.enriched.travels_to_ip_index >= 0); + bool result = false; + + typename turn_operation_type::comparable_distance_type + max_remaining_distance = 0; - if (info.seg_id.source_index == 0) + selected_op_index = -1; + for (int i = 0; i < 2; i++) { - geometry::copy_segments<Reverse1>(g1, info.seg_id, - info.enriched.travels_to_vertex_index, - robust_policy, - current_output); + turn_operation_type const& op = turn.operations[i]; + if (op.visited.started()) + { + selected_op_index = i; + return true; + } + + signed_size_type const next_turn_index = get_next_turn_index(op); + + // In some cases there are two alternatives. + // For "ii", take the other one (alternate) + // UNLESS the other one is already visited + // For "uu", take the same one (see above); + // For "cc", take either one, but if there is a starting one, + // take that one. If next is dead end, skip that one. + if ( (op.operation == operation_continue + && traverse_possible(next_turn_index) + && ! result) + || (op.operation == OperationType + && ! op.visited.finished() + && (! result + || select_source(next_turn_index, op.seg_id, seg_id) + ) + ) + ) + { + if (op.operation == operation_continue) + { + max_remaining_distance = op.remaining_distance; + } + selected_op_index = i; + debug_traverse(turn, op, " Candidate"); + result = true; + } + + if (op.operation == operation_continue && result) + { + if (next_turn_index == start_turn_index) + { + selected_op_index = i; + debug_traverse(turn, op, " Candidate override (start)"); + } + else if (op.remaining_distance > max_remaining_distance) + { + max_remaining_distance = op.remaining_distance; + selected_op_index = i; + debug_traverse(turn, op, " Candidate override (remaining)"); + } + } } - else + + if (result) { - geometry::copy_segments<Reverse2>(g2, info.seg_id, - info.enriched.travels_to_vertex_index, - robust_policy, - current_output); + debug_traverse(turn, turn.operations[selected_op_index], " Accepted"); } - seg_id = info.seg_id; - ip = boost::begin(turns) + info.enriched.travels_to_ip_index; + + return result; } - else + + inline bool select_from_cluster(signed_size_type& turn_index, + int& op_index, signed_size_type start_turn_index, + sbs_type const& sbs, bool allow_pass_rank) { - ip = boost::begin(turns) + info.enriched.next_ip_index; - seg_id = info.seg_id; + bool const is_union = OperationType == operation_union; + bool const is_intersection = OperationType == operation_intersection; + + std::size_t selected_rank = 0; + std::size_t min_rank = 0; + bool result = false; + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[i]; + if (result && ranked_point.main_rank > selected_rank) + { + return result; + } + + turn_type const& ranked_turn = m_turns[ranked_point.turn_index]; + turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.op_index]; + + if (result && ranked_op.visited.finalized()) + { + // One of the arcs in the same direction as the selected result + // is already traversed. + return false; + } + + if (! allow_pass_rank && ranked_op.visited.finalized()) + { + // Skip this one, go to next + min_rank = ranked_point.main_rank; + continue; + } + + if (ranked_point.index == sort_by_side::index_to + && (ranked_point.main_rank > min_rank + || ranked_turn.both(operation_continue))) + { + if ((is_union + && ranked_op.enriched.count_left == 0 + && ranked_op.enriched.count_right > 0) + || (is_intersection + && ranked_op.enriched.count_right == 2)) + { + if (result && ranked_point.turn_index != start_turn_index) + { + // Don't override - only override if arrive at start + continue; + } + + turn_index = ranked_point.turn_index; + op_index = ranked_point.op_index; + + if (is_intersection + && ranked_turn.both(operation_intersection) + && ranked_op.visited.finalized()) + { + // Override: + // For a ii turn, even though one operation might be selected, + // it should take the other one if the first one is used in a completed ring + op_index = 1 - ranked_point.op_index; + } + + result = true; + selected_rank = ranked_point.main_rank; + } + else if (! allow_pass_rank) + { + return result; + } + } + } + return result; } - detail::overlay::append_no_dups_or_spikes(current_output, ip->point, - robust_policy); + inline bool select_turn_from_cluster(signed_size_type& turn_index, + int& op_index, signed_size_type start_turn_index, + point_type const& point) + { + bool const is_union = OperationType == operation_union; - return true; -} + turn_type const& turn = m_turns[turn_index]; + BOOST_ASSERT(turn.cluster_id >= 0); + typename Clusters::const_iterator mit = m_clusters.find(turn.cluster_id); + BOOST_ASSERT(mit != m_clusters.end()); -inline bool select_source(operation_type operation, - signed_size_type source1, - signed_size_type source2) -{ - return (operation == operation_intersection && source1 != source2) - || (operation == operation_union && source1 == source2) - ; -} + std::set<signed_size_type> const& ids = mit->second; + sbs_type sbs; + sbs.set_origin(point); -template -< - typename Turn, - typename Iterator -> -inline bool select_next_ip(operation_type operation, - Turn& turn, - segment_identifier const& seg_id, - Iterator& selected) -{ - if (turn.discarded) + for (typename std::set<signed_size_type>::const_iterator sit = ids.begin(); + sit != ids.end(); ++sit) + { + signed_size_type cluster_turn_index = *sit; + turn_type const& cluster_turn = m_turns[cluster_turn_index]; + if (cluster_turn.discarded) + { + // Defensive check, discarded turns should not be in cluster + continue; + } + + for (int i = 0; i < 2; i++) + { + sbs.add(cluster_turn.operations[i], cluster_turn_index, i, + m_geometry1, m_geometry2, false); + } + } + + sbs.apply(turn.point); + + int open_count = 0; + if (is_union) + { + // Check how many open spaces there are. + // TODO: might be moved to sbs itself, though it also uses turns + + std::size_t last_rank = 0; + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[i]; + + if (ranked_point.main_rank > last_rank + && ranked_point.index == sort_by_side::index_to) + { + turn_type const& ranked_turn = m_turns[ranked_point.turn_index]; + turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.op_index]; + if (ranked_op.enriched.count_left == 0 + && ranked_op.enriched.count_right > 0) + { + open_count++; + last_rank = ranked_point.main_rank; + } + } + } + } + + bool allow = false; + if (open_count > 1) + { + sbs.reverse(); + allow = true; + } + + return select_from_cluster(turn_index, op_index, start_turn_index, sbs, allow); + } + + inline void change_index_for_self_turn(signed_size_type& to_vertex_index, + turn_type const& start_turn, + turn_operation_type const& start_op, + int start_op_index) { - return false; + turn_operation_type const& other_op + = start_turn.operations[1 - start_op_index]; + if (start_op.seg_id.source_index != other_op.seg_id.source_index) + { + // Not a buffer/self-turn + return; + } + + // It travels to itself, can happen. If this is a buffer, it can + // sometimes travel to itself in the following configuration: + // + // +---->--+ + // | | + // | +---*----+ *: one turn, with segment index 2/7 + // | | | | + // | +---C | C: closing point (start/end) + // | | + // +------------+ + // + // If it starts on segment 2 and travels to itself on segment 2, that + // should be corrected to 7 because that is the shortest path + // + // Also a uu turn (touching with another buffered ring) might have this + // apparent configuration, but there it should + // always travel the whole ring + + bool const correct + = ! start_turn.both(operation_union) + && start_op.seg_id.segment_index == to_vertex_index; + +#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSE) + std::cout << " WARNING: self-buffer " + << " correct=" << correct + << " turn=" << operation_char(start_turn.operations[0].operation) + << operation_char(start_turn.operations[1].operation) + << " start=" << start_op.seg_id.segment_index + << " from=" << to_vertex_index + << " to=" << other_op.enriched.travels_to_vertex_index + << std::endl; +#endif + + if (correct) + { + to_vertex_index = other_op.enriched.travels_to_vertex_index; + } } - bool has_tp = false; + template <typename Ring> + inline traverse_error_type travel_to_next_turn(signed_size_type start_turn_index, + int start_op_index, + signed_size_type& turn_index, + int& op_index, + segment_identifier& seg_id, + Ring& current_ring, + bool is_start) + { + int const previous_op_index = op_index; + signed_size_type const previous_turn_index = turn_index; + turn_type& previous_turn = m_turns[turn_index]; + turn_operation_type& previous_op = previous_turn.operations[op_index]; - typedef typename std::iterator_traits - < - Iterator - >::value_type operation_type; + // If there is no next IP on this segment + if (previous_op.enriched.next_ip_index < 0) + { + if (previous_op.enriched.travels_to_vertex_index < 0 + || previous_op.enriched.travels_to_ip_index < 0) + { + return is_start + ? traverse_error_no_next_ip_at_start + : traverse_error_no_next_ip; + } + + signed_size_type to_vertex_index = previous_op.enriched.travels_to_vertex_index; + + if (is_start && + previous_op.enriched.travels_to_ip_index == start_turn_index) + { + change_index_for_self_turn(to_vertex_index, previous_turn, + previous_op, start_op_index); + } + + if (previous_op.seg_id.source_index == 0) + { + geometry::copy_segments<Reverse1>(m_geometry1, + previous_op.seg_id, to_vertex_index, + m_robust_policy, current_ring); + } + else + { + geometry::copy_segments<Reverse2>(m_geometry2, + previous_op.seg_id, to_vertex_index, + m_robust_policy, current_ring); + } + seg_id = previous_op.seg_id; + turn_index = previous_op.enriched.travels_to_ip_index; + } + else + { + turn_index = previous_op.enriched.next_ip_index; + seg_id = previous_op.seg_id; + } + + // turn_index is not yet finally selected, can change for clusters + bool const has_cluster = m_turns[turn_index].cluster_id >= 0; + if (has_cluster) + { + + if (! select_turn_from_cluster(turn_index, op_index, + start_turn_index, current_ring.back())) + { + return is_start + ? traverse_error_no_next_ip_at_start + : traverse_error_no_next_ip; + } + + if (is_start && turn_index == previous_turn_index) + { + op_index = previous_op_index; + } + } + + turn_type& current_turn = m_turns[turn_index]; + detail::overlay::append_no_dups_or_spikes(current_ring, current_turn.point, + m_robust_policy); + + if (is_start) + { + // Register the start + previous_op.visited.set_started(); + m_visitor.visit_traverse(m_turns, previous_turn, previous_op, "Start"); + } + + if (! has_cluster) + { + if (! select_operation(current_turn, + start_turn_index, + seg_id, + op_index)) + { + return is_start + ? traverse_error_dead_end_at_start + : traverse_error_dead_end; + } + } + + turn_operation_type& op = current_turn.operations[op_index]; + if (op.visited.finalized() || op.visited.visited()) + { + return traverse_error_visit_again; + } - typename operation_type::comparable_distance_type - max_remaining_distance = 0; + // Register the visit + set_visited(current_turn, op); + m_visitor.visit_traverse(m_turns, current_turn, op, "Visit"); - selected = boost::end(turn.operations); - for (Iterator it = boost::begin(turn.operations); - it != boost::end(turn.operations); - ++it) + return traverse_error_none; + } + + inline void finalize_visit_info() { - if (it->visited.started()) + for (typename boost::range_iterator<Turns>::type + it = boost::begin(m_turns); + it != boost::end(m_turns); + ++it) { - selected = it; - //std::cout << " RETURN"; - return true; + turn_type& turn = *it; + for (int i = 0; i < 2; i++) + { + turn_operation_type& op = turn.operations[i]; + op.visited.finalize(); + } } + } - // In some cases there are two alternatives. - // For "ii", take the other one (alternate) - // UNLESS the other one is already visited - // For "uu", take the same one (see above); - // For "cc", take either one, but if there is a starting one, - // take that one. - if ( (it->operation == operation_continue - && (! has_tp || it->visited.started() - ) - ) - || (it->operation == operation - && ! it->visited.finished() - && (! has_tp - || select_source(operation, - it->seg_id.source_index, seg_id.source_index) - ) - ) - ) + inline void set_visited(turn_type& turn, turn_operation_type& op) + { + // On "continue", set "visited" for ALL directions in this turn + if (op.operation == detail::overlay::operation_continue) { - if (it->operation == operation_continue) + for (int i = 0; i < 2; i++) { - max_remaining_distance = it->remaining_distance; + turn_operation_type& op = turn.operations[i]; + if (op.visited.none()) + { + op.visited.set_visited(); + } } - selected = it; - debug_traverse(turn, *it, " Candidate"); - has_tp = true; + } + else + { + op.visited.set_visited(); + } + } + + + template <typename Ring> + inline traverse_error_type traverse(Ring& ring, + signed_size_type start_turn_index, int start_op_index) + { + turn_type const& start_turn = m_turns[start_turn_index]; + turn_operation_type& start_op = m_turns[start_turn_index].operations[start_op_index]; + + detail::overlay::append_no_dups_or_spikes(ring, start_turn.point, + m_robust_policy); + + signed_size_type current_turn_index = start_turn_index; + int current_op_index = start_op_index; + segment_identifier current_seg_id; + + traverse_error_type error = travel_to_next_turn(start_turn_index, + start_op_index, + current_turn_index, current_op_index, current_seg_id, + ring, true); + + if (error != traverse_error_none) + { + // This is not necessarily a problem, it happens for clustered turns + // which are "build in" or otherwise point inwards + return error; + } + + if (current_turn_index == start_turn_index) + { + start_op.visited.set_finished(); + m_visitor.visit_traverse(m_turns, m_turns[current_turn_index], start_op, "Early finish"); + return traverse_error_none; } - if (it->operation == operation_continue && has_tp) + std::size_t const max_iterations = 2 + 2 * m_turns.size(); + for (std::size_t i = 0; i <= max_iterations; i++) { - if (it->remaining_distance > max_remaining_distance) + // We assume clockwise polygons only, non self-intersecting, closed. + // However, the input might be different, and checking validity + // is up to the library user. + + // Therefore we make here some sanity checks. If the input + // violates the assumptions, the output polygon will not be correct + // but the routine will stop and output the current polygon, and + // will continue with the next one. + + // Below three reasons to stop. + error = travel_to_next_turn(start_turn_index, start_op_index, + current_turn_index, current_op_index, current_seg_id, + ring, false); + + if (error != traverse_error_none) { - max_remaining_distance = it->remaining_distance; - selected = it; - debug_traverse(turn, *it, " Candidate override"); + return error; + } + + if (current_turn_index == start_turn_index + && current_op_index == start_op_index) + { + start_op.visited.set_finished(); + m_visitor.visit_traverse(m_turns, start_turn, start_op, "Finish"); + return traverse_error_none; } } + + return traverse_error_endless_loop; } - if (has_tp) + template <typename Rings> + void traverse_with_operation(turn_type const& start_turn, + std::size_t turn_index, int op_index, + Rings& rings, std::size_t& finalized_ring_size, + typename Backtrack::state_type& state) { - debug_traverse(turn, *selected, " Accepted"); + typedef typename boost::range_value<Rings>::type ring_type; + + turn_operation_type const& start_op = start_turn.operations[op_index]; + + if (! start_op.visited.none() + || ! start_op.enriched.startable + || start_op.visited.rejected() + || ! (start_op.operation == OperationType + || start_op.operation == detail::overlay::operation_continue)) + { + return; + } + + ring_type ring; + traverse_error_type traverse_error = traverse(ring, turn_index, op_index); + + if (traverse_error == traverse_error_none) + { + std::size_t const min_num_points + = core_detail::closure::minimum_ring_size + < + geometry::closure<ring_type>::value + >::value; + + if (geometry::num_points(ring) >= min_num_points) + { + clean_closing_dups_and_spikes(ring, m_robust_policy); + rings.push_back(ring); + + finalize_visit_info(); + finalized_ring_size++; + } + } + else + { + Backtrack::apply( + finalized_ring_size, + rings, ring, m_turns, start_turn, + m_turns[turn_index].operations[op_index], + traverse_error, + m_geometry1, m_geometry2, m_robust_policy, + state, m_visitor); + } } + template <typename Rings> + void iterate(Rings& rings, std::size_t& finalized_ring_size, + typename Backtrack::state_type& state, + int pass) + { + if (pass == 1) + { + if (OperationType == operation_intersection) + { + // Second pass currently only used for uu + return; + } + if (! m_has_uu) + { + // There is no uu found in first pass + return; + } + if (m_has_only_uu) + { + m_switch_at_uu = false; + } + } + + // Iterate through all unvisited points + for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index) + { + turn_type const& start_turn = m_turns[turn_index]; - return has_tp; -} + if (start_turn.discarded || start_turn.blocked()) + { + // Skip discarded and blocked turns + continue; + } + if (OperationType == operation_union) + { + if (start_turn.both(operation_union)) + { + // Start with a uu-turn only in the second pass + m_has_uu = true; + if (pass == 0) + { + continue; + } + } + else + { + m_has_only_uu = false; + } + } + + for (int op_index = 0; op_index < 2; op_index++) + { + traverse_with_operation(start_turn, turn_index, op_index, + rings, finalized_ring_size, state); + } + } + } +private : + Geometry1 const& m_geometry1; + Geometry2 const& m_geometry2; + Turns& m_turns; + Clusters const& m_clusters; + RobustPolicy const& m_robust_policy; + Visitor& m_visitor; + + // Next members are only used for operation union + bool m_has_uu; + bool m_has_only_uu; + bool m_switch_at_uu; +}; /*! @@ -256,185 +811,45 @@ template bool Reverse1, bool Reverse2, typename Geometry1, typename Geometry2, + operation_type OperationType, typename Backtrack = backtrack_check_self_intersections<Geometry1, Geometry2> > class traverse { public : - template <typename RobustPolicy, typename Turns, typename Rings> + template + < + typename RobustPolicy, + typename Turns, + typename Rings, + typename Visitor, + typename Clusters + > static inline void apply(Geometry1 const& geometry1, Geometry2 const& geometry2, - detail::overlay::operation_type operation, RobustPolicy const& robust_policy, - Turns& turns, Rings& rings) + Turns& turns, Rings& rings, + Clusters const& clusters, + Visitor& visitor) { - typedef typename boost::range_value<Rings>::type ring_type; - typedef typename boost::range_iterator<Turns>::type turn_iterator; - typedef typename boost::range_value<Turns>::type turn_type; - typedef typename boost::range_iterator + traversal < - typename turn_type::container_type - >::type turn_operation_iterator_type; + Reverse1, Reverse2, OperationType, + Geometry1, Geometry2, + Turns, Clusters, + RobustPolicy, Visitor, + Backtrack + > trav(geometry1, geometry2, turns, clusters, + robust_policy, visitor); - std::size_t const min_num_points - = core_detail::closure::minimum_ring_size - < - geometry::closure<ring_type>::value - >::value; - - std::size_t size_at_start = boost::size(rings); + std::size_t finalized_ring_size = boost::size(rings); typename Backtrack::state_type state; - do - { - state.reset(); - // Iterate through all unvisited points - for (turn_iterator it = boost::begin(turns); - state.good() && it != boost::end(turns); - ++it) - { - // Skip discarded ones - if (! (it->discarded || ! it->selectable_start || it->blocked())) - { - for (turn_operation_iterator_type iit = boost::begin(it->operations); - state.good() && iit != boost::end(it->operations); - ++iit) - { - if (iit->visited.none() - && ! iit->visited.rejected() - && (iit->operation == operation - || iit->operation == detail::overlay::operation_continue) - ) - { - set_visited_for_continue(*it, *iit); - - ring_type current_output; - detail::overlay::append_no_dups_or_spikes(current_output, - it->point, robust_policy); - - turn_iterator current = it; - turn_operation_iterator_type current_iit = iit; - segment_identifier current_seg_id; - - if (! detail::overlay::assign_next_ip<Reverse1, Reverse2>( - geometry1, geometry2, - turns, - current, current_output, - *iit, current_seg_id, - robust_policy)) - { - Backtrack::apply( - size_at_start, - rings, current_output, turns, *current_iit, - "No next IP", - geometry1, geometry2, robust_policy, state); - } - - if (! detail::overlay::select_next_ip( - operation, - *current, - current_seg_id, - current_iit)) - { - Backtrack::apply( - size_at_start, - rings, current_output, turns, *iit, - "Dead end at start", - geometry1, geometry2, robust_policy, state); - } - else - { - - iit->visited.set_started(); - detail::overlay::debug_traverse(*it, *iit, "-> Started"); - detail::overlay::debug_traverse(*current, *current_iit, "Selected "); - - - typename boost::range_size<Turns>::type i = 0; - - while (current_iit != iit && state.good()) - { - if (current_iit->visited.visited()) - { - // It visits a visited node again, without passing the start node. - // This makes it suspicious for endless loops - Backtrack::apply( - size_at_start, - rings, current_output, turns, *iit, - "Visit again", - geometry1, geometry2, robust_policy, state); - } - else - { - - - // We assume clockwise polygons only, non self-intersecting, closed. - // However, the input might be different, and checking validity - // is up to the library user. - - // Therefore we make here some sanity checks. If the input - // violates the assumptions, the output polygon will not be correct - // but the routine will stop and output the current polygon, and - // will continue with the next one. - - // Below three reasons to stop. - detail::overlay::assign_next_ip<Reverse1, Reverse2>( - geometry1, geometry2, - turns, current, current_output, - *current_iit, current_seg_id, - robust_policy); - - if (! detail::overlay::select_next_ip( - operation, - *current, - current_seg_id, - current_iit)) - { - // Should not occur in valid (non-self-intersecting) polygons - // Should not occur in self-intersecting polygons without spikes - // Might occur in polygons with spikes - Backtrack::apply( - size_at_start, - rings, current_output, turns, *iit, - "Dead end", - geometry1, geometry2, robust_policy, state); - } - else - { - detail::overlay::debug_traverse(*current, *current_iit, "Selected "); - } - - if (i++ > 2 + 2 * turns.size()) - { - // Sanity check: there may be never more loops - // than turn points. - // Turn points marked as "ii" can be visited twice. - Backtrack::apply( - size_at_start, - rings, current_output, turns, *iit, - "Endless loop", - geometry1, geometry2, robust_policy, state); - } - } - } - - if (state.good()) - { - iit->visited.set_finished(); - detail::overlay::debug_traverse(*current, *iit, "->Finished"); - if (geometry::num_points(current_output) >= min_num_points) - { - clean_closing_dups_and_spikes(current_output, robust_policy); - rings.push_back(current_output); - } - } - } - } - } - } - } - } while (! state.good()); + for (int pass = 0; pass < 2; pass++) + { + trav.iterate(rings, finalized_ring_size, state, pass); + } } }; diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp index 1ac77d5796..699997fc38 100644 --- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp @@ -58,6 +58,8 @@ enum method_type template <typename Point, typename SegmentRatio> struct turn_operation { + typedef SegmentRatio segment_ratio_type; + operation_type operation; segment_identifier seg_id; SegmentRatio fraction; @@ -91,23 +93,25 @@ template struct turn_info { typedef Point point_type; + typedef SegmentRatio segment_ratio_type; typedef Operation turn_operation_type; typedef Container container_type; Point point; method_type method; + int cluster_id; bool discarded; - bool selectable_start; // Can be used as starting-turn in traverse bool colocated; - + bool switch_source; // For u/u turns which can either switch or not Container operations; inline turn_info() : method(method_none) + , cluster_id(-1) , discarded(false) - , selectable_start(true) , colocated(false) + , switch_source(false) {} inline bool both(operation_type type) const diff --git a/boost/geometry/algorithms/detail/overlay/visit_info.hpp b/boost/geometry/algorithms/detail/overlay/visit_info.hpp index 4284a801a1..9e1e6b9056 100644 --- a/boost/geometry/algorithms/detail/overlay/visit_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/visit_info.hpp @@ -28,11 +28,13 @@ private : int m_visit_code; bool m_rejected; + bool m_final; public: inline visit_info() : m_visit_code(0) , m_rejected(false) + , m_final(false) {} inline void set_visited() { m_visit_code = VISITED; } @@ -49,15 +51,24 @@ public: inline bool started() const { return m_visit_code == STARTED; } inline bool finished() const { return m_visit_code == FINISHED; } inline bool rejected() const { return m_rejected; } + inline bool finalized() const { return m_final; } inline void clear() { - if (! rejected()) + if (! m_rejected && ! m_final) { m_visit_code = NONE; } } + inline void finalize() + { + if (visited() || started() || finished() ) + { + m_final = true; + } + } + #ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION friend std::ostream& operator<<(std::ostream &os, visit_info const& v) { diff --git a/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp b/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp index 399c8fbefc..607ba81531 100644 --- a/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp +++ b/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp @@ -17,6 +17,7 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_POINT_IS_EQUAL_OR_SPIKE_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_POINT_IS_EQUAL_OR_SPIKE_HPP +#include <boost/geometry/algorithms/detail/direction_code.hpp> #include <boost/geometry/algorithms/detail/recalculate.hpp> #include <boost/geometry/policies/robustness/robust_point_type.hpp> #include <boost/geometry/strategies/side.hpp> @@ -31,17 +32,6 @@ namespace boost { namespace geometry namespace detail { -template <std::size_t Index, typename Point1, typename Point2> -inline int sign_of_difference(Point1 const& point1, Point2 const& point2) -{ - return - math::equals(geometry::get<Index>(point1), geometry::get<Index>(point2)) - ? - 0 - : - (geometry::get<Index>(point1) > geometry::get<Index>(point2) ? 1 : -1); -} - // Checks if a point ("last_point") causes a spike w.r.t. // the specified two other points (segment_a, segment_b) // diff --git a/boost/geometry/algorithms/detail/recalculate.hpp b/boost/geometry/algorithms/detail/recalculate.hpp index 056f7c6e1d..fd3ed6b8da 100644 --- a/boost/geometry/algorithms/detail/recalculate.hpp +++ b/boost/geometry/algorithms/detail/recalculate.hpp @@ -21,8 +21,10 @@ #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/bounds.hpp> #include <boost/numeric/conversion/cast.hpp> -#include <boost/range.hpp> -#include <boost/type_traits.hpp> +#include <boost/range/begin.hpp> +#include <boost/range/end.hpp> +#include <boost/range/iterator.hpp> +#include <boost/range/size.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/algorithms/append.hpp> diff --git a/boost/geometry/algorithms/detail/relate/result.hpp b/boost/geometry/algorithms/detail/relate/result.hpp index b7b9264449..a92badf65b 100644 --- a/boost/geometry/algorithms/detail/relate/result.hpp +++ b/boost/geometry/algorithms/detail/relate/result.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2013, 2014, 2015. -// Modifications copyright (c) 2013-2015 Oracle and/or its affiliates. +// This file was modified by Oracle on 2013-2016. +// Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -133,9 +133,6 @@ public: matrix_handler() {} - matrix_handler(Matrix const&) - {} - result_type const& result() const { return m_matrix; @@ -604,10 +601,6 @@ public: bool interrupt; - inline mask_handler() - : interrupt(false) - {} - inline explicit mask_handler(Mask const& m) : interrupt(false) , m_mask(m) |